Merge pull request #2051 from atgeirr/enable-onephase-revised-rebased

Enable onephase

CMakeLists.txt

@@ -19,6 +19,7 @@ cmake_minimum_required (VERSION 2.8)
 option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON)
 set( USE_OPENMP_DEFAULT OFF ) # Use of OpenMP is considered experimental
 option(BUILD_FLOW "Build the production oriented flow simulator?" ON)
+option(BUILD_FLOW_VARIANTS "Build the variants for flow by default?" OFF)
 option(BUILD_EBOS "Build the research oriented ebos simulator?" ON)
 option(BUILD_EBOS_EXTENSIONS "Build the variants for various extensions of ebos by default?" OFF)
 option(BUILD_EBOS_DEBUG_EXTENSIONS "Build the ebos variants which are purely for debugging by default?" OFF)
@@ -155,15 +156,36 @@ opm_add_test(flow
   flow/flow_ebos_oilwater_polymer.cpp
   flow/flow_ebos_oilwater_polymer_injectivity.cpp)
 
+if (NOT BUILD_FLOW_VARIANTS)
+  set(FLOW_VARIANTS_DEFAULT_ENABLE_IF "FALSE")
+else()
+  set(FLOW_VARIANTS_DEFAULT_ENABLE_IF "TRUE")
+endif()
+
+# Variant versions of Flow.
 opm_add_test(flow_blackoil_dunecpr
   ONLY_COMPILE
-  DEFAULT_ENABLE_IF ${FLOW_DEFAULT_ENABLE_IF}
+  DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}
   SOURCES flow/flow_blackoil_dunecpr.cpp
   EXE_NAME flow_blackoil_dunecpr
   DEPENDS "opmsimulators"
   LIBRARIES "opmsimulators")
 
+opm_add_test(flow_onephase
+  ONLY_COMPILE
+  DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}
+  SOURCES flow/flow_onephase.cpp
+  EXE_NAME flow_onephase
+  DEPENDS "opmsimulators"
+  LIBRARIES "opmsimulators")
 
+opm_add_test(flow_onephase_energy
+  ONLY_COMPILE
+  DEFAULT_ENABLE_IF ${FLOW_VARIANTS_DEFAULT_ENABLE_IF}
+  SOURCES flow/flow_onephase_energy.cpp
+  EXE_NAME flow_onephase_energy
+  DEPENDS "opmsimulators"
+  LIBRARIES "opmsimulators")
 
 if (BUILD_FLOW)
   install(TARGETS flow DESTINATION bin)

ebos/ecloutputblackoilmodule.hh

@@ -395,7 +395,16 @@ public:
             }
 
             if (oilPressure_.size() > 0) {
-                oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx));
+                if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                    oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(oilPhaseIdx));
+                }else{
+                    // put pressure in oil pressure for output
+                    if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
+                        oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(waterPhaseIdx));
+                    } else {
+                        oilPressure_[globalDofIdx] = Opm::getValue(fs.pressure(gasPhaseIdx));
+                    }
+                }
                 Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]);
             }
 

ebos/eclproblem.hh

@@ -665,17 +665,20 @@ public:
             minTimeStepSize_ = tuning.getTSMINZ(0);
         }
 
+        initFluidSystem_();
+
         // deal with DRSDT
         unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
-        maxDRs_.resize(ntpvt, 1e30);
-        dRsDtOnlyFreeGas_.resize(ntpvt, false);
         size_t numDof = this->model().numGridDof();
-        lastRs_.resize(numDof, 0.0);
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-        maxDRv_.resize(ntpvt, 1e30);
+            maxDRs_.resize(ntpvt, 1e30);
-        lastRv_.resize(numDof, 0.0);
+            dRsDtOnlyFreeGas_.resize(ntpvt, false);
-        maxOilSaturation_.resize(numDof, 0.0);
+            lastRs_.resize(numDof, 0.0);
+            maxDRv_.resize(ntpvt, 1e30);
+            lastRv_.resize(numDof, 0.0);
+            maxOilSaturation_.resize(numDof, 0.0);
+        }
 
-        initFluidSystem_();
         updateElementDepths_();
         readRockParameters_();
         readMaterialParameters_();

flow/flow_blackoil_dunecpr.cpp

@@ -43,7 +43,7 @@ SET_PROP(EclFlowProblemSimple, FluidState)
       typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
       typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
       static const bool compositionSwitchEnabled = Indices::gasEnabled;
- 
+
     public:
 //typedef Opm::BlackOilFluidSystemSimple<Scalar> type;
        typedef Opm::BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled,  Indices::numPhases > type;
@@ -89,7 +89,7 @@ namespace Opm {
 #endif
     SET_BOOL_PROP(EclFlowProblemSimple, EnableStorageCache, true);
     SET_BOOL_PROP(EclFlowProblemSimple, EnableIntensiveQuantityCache, true);
-    
+
     //SET_INT_PROP(EclFlowProblemSimple, NumWellAdjoint, 1);
     //SET_BOOL_PROP(EclFlowProblem, EnableStorageCache, true);
     //SET_BOOL_PROP(EclFlowProblem, EnableIntensiveQuantityCache, true);
@@ -98,8 +98,6 @@ namespace Opm {
 
 int main(int argc, char** argv)
 {
-  typedef TTAG(EclFlowProblemSimple) TypeTag;
+    typedef TTAG(EclFlowProblemSimple) TypeTag;
-  bool outputCout = true;
+    return mainFlow<TypeTag>(argc, argv);
-  bool outputFiles = true;
-  return mainFlow<TypeTag>(argc, argv, outputCout, outputFiles);
 }

flow/flow_onephase.cpp

@@ -0,0 +1,91 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015, 2017 IRIS AS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include "config.h"
+#include "flow/flow_tag.hpp"
+#include <opm/models/blackoil/blackoilonephaseindices.hh>
+
+
+BEGIN_PROPERTIES
+NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem));
+NEW_PROP_TAG(FluidState);
+NEW_PROP_TAG(FluidSystem);
+//! The indices required by the model
+SET_PROP(EclFlowProblemSimple, Indices)
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    typedef TTAG(EclFlowProblem) BaseTypeTag;
+    typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem;
+
+public:
+    typedef Ewoms::BlackOilOnePhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent),
+                                           GET_PROP_VALUE(TypeTag, EnablePolymer),
+                                           GET_PROP_VALUE(TypeTag, EnableEnergy),
+                                           GET_PROP_VALUE(TypeTag, EnableFoam),
+                                           /*PVOffset=*/0,
+                                           /*enebledCompIdx=*/FluidSystem::waterCompIdx>
+        type;
+};
+SET_PROP(EclFlowProblemSimple, FluidState)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    static const bool compositionSwitchEnabled = Indices::gasEnabled;
+
+public:
+    // typedef Opm::BlackOilFluidSystemSimple<Scalar> type;
+    typedef Opm::BlackOilFluidState<Evaluation,
+                                    FluidSystem,
+                                    enableTemperature,
+                                    enableEnergy,
+                                    compositionSwitchEnabled,
+                                    Indices::numPhases>
+        type;
+};
+
+// SET_PROP(EclFlowProblemSimple, FluidSystem)
+// {
+// private:
+//   //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//   typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//   typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+//   typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+// public:
+//   typedef Opm::BlackOilFluidSystem<Scalar,Indices> type;
+// };
+END_PROPERTIES
+
+int
+main(int argc, char** argv)
+{
+    typedef TTAG(EclFlowProblemSimple) TypeTag;
+    return mainFlow<TypeTag>(argc, argv);
+}

flow/flow_onephase_energy.cpp

@@ -0,0 +1,93 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015, 2017 IRIS AS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include "config.h"
+#include "flow/flow_tag.hpp"
+#include <opm/models/blackoil/blackoilonephaseindices.hh>
+
+
+BEGIN_PROPERTIES
+NEW_TYPE_TAG(EclFlowProblemSimple, INHERITS_FROM(EclFlowProblem));
+SET_BOOL_PROP(EclFlowProblemSimple, EnableEnergy, true);
+NEW_PROP_TAG(FluidState);
+NEW_PROP_TAG(FluidSystem);
+//! The indices required by the model
+SET_PROP(EclFlowProblemSimple, Indices)
+{
+private:
+    // it is unfortunately not possible to simply use 'TypeTag' here because this leads
+    // to cyclic definitions of some properties. if this happens the compiler error
+    // messages unfortunately are *really* confusing and not really helpful.
+    typedef TTAG(EclFlowProblem) BaseTypeTag;
+    typedef typename GET_PROP_TYPE(BaseTypeTag, FluidSystem) FluidSystem;
+
+public:
+    typedef Ewoms::BlackOilOnePhaseIndices<GET_PROP_VALUE(TypeTag, EnableSolvent),
+                                           GET_PROP_VALUE(TypeTag, EnablePolymer),
+                                           GET_PROP_VALUE(TypeTag, EnableEnergy),
+                                           GET_PROP_VALUE(TypeTag, EnableFoam),
+                                           /*PVOffset=*/0,
+                                           /*enebledCompIdx=*/FluidSystem::waterCompIdx>
+        type;
+};
+SET_PROP(EclFlowProblemSimple, FluidState)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
+    enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    static const bool compositionSwitchEnabled = Indices::gasEnabled;
+
+public:
+    // typedef Opm::BlackOilFluidSystemSimple<Scalar> type;
+    typedef Opm::BlackOilFluidState<Evaluation,
+                                    FluidSystem,
+                                    enableTemperature,
+                                    enableEnergy,
+                                    compositionSwitchEnabled,
+                                    Indices::numPhases>
+        type;
+};
+
+// SET_PROP(EclFlowProblemSimple, FluidSystem)
+// {
+// private:
+//   //typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//   typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+//   typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+//   typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+// public:
+//   typedef Opm::BlackOilFluidSystem<Scalar,Indices> type;
+// };
+END_PROPERTIES
+
+
+int
+main(int argc, char** argv)
+{
+    typedef TTAG(EclFlowProblemSimple) TypeTag;
+    return mainFlow<TypeTag>(argc, argv);
+}

flow/flow_tag.hpp

@@ -29,13 +29,18 @@
 #include <opm/common/utility/parameters/ParameterGroup.hpp>
 #include <opm/material/common/ResetLocale.hpp>
 
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/OpmLog/EclipsePRTLog.hpp>
+#include <opm/common/OpmLog/LogUtil.hpp>
+
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ArrayDimChecker.hpp>
 
 //#include <opm/material/fluidsystems/BlackOilFluidSystemSimple.hpp>
-//#include <opm/material/fluidsystems/BlackOilFluidSystemSimple.hpp> 
+//#include <opm/material/fluidsystems/BlackOilFluidSystemSimple.hpp>
 #include <opm/models/blackoil/blackoilintensivequantities.hh>
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 //#include <opm/material/fluidstates/BlackOilFluidStateSimple.hpp>
@@ -66,7 +71,7 @@ namespace Opm {
     Vanguard::setExternalDeck(&deck);
     Vanguard::setExternalEclState(&eclState);
   }
-  
+
 // ----------------- Main program -----------------
   template <class TypeTag>
   int flowEbosMain(int argc, char** argv, bool outputCout, bool outputFiles)
@@ -74,7 +79,7 @@ namespace Opm {
     // we always want to use the default locale, and thus spare us the trouble
     // with incorrect locale settings.
     Opm::resetLocale();
-    
+
 #if HAVE_DUNE_FEM
     Dune::Fem::MPIManager::initialize(argc, argv);
 #else
@@ -118,13 +123,162 @@ namespace detail
 
         throw std::invalid_argument( "Cannot find input case " + casename );
     }
+
+
+    // This function is an extreme special case, if the program has been invoked
+    // *exactly* as:
+    //
+    //    flow   --version
+    //
+    // the call is intercepted by this function which will print "flow $version"
+    // on stdout and exit(0).
+    void handleVersionCmdLine(int argc, char** argv) {
+        for ( int i = 1; i < argc; ++i )
+        {
+            if (std::strcmp(argv[i], "--version") == 0) {
+                std::cout << "flow " << Opm::moduleVersionName() << std::endl;
+                std::exit(EXIT_SUCCESS);
+            }
+        }
+    }
+
+}
+
+
+
+enum class FileOutputMode {
+    //! \brief No output to files.
+    OUTPUT_NONE = 0,
+    //! \brief Output only to log files, no eclipse output.
+    OUTPUT_LOG_ONLY = 1,
+    //! \brief Output to all files.
+    OUTPUT_ALL = 3
+};
+
+
+
+void ensureOutputDirExists(const std::string& cmdline_output_dir)
+{
+    if (!boost::filesystem::is_directory(cmdline_output_dir)) {
+        try {
+            boost::filesystem::create_directories(cmdline_output_dir);
+        }
+        catch (...) {
+            throw std::runtime_error("Creation of output directory '" + cmdline_output_dir + "' failed\n");
+        }
+    }
+}
+
+
+// Setup the OpmLog backends
+FileOutputMode setupLogging(int mpi_rank_, const std::string& deck_filename, const std::string& cmdline_output_dir, const std::string& cmdline_output, bool output_cout_, const std::string& stdout_log_id) {
+
+    if (!cmdline_output_dir.empty()) {
+        ensureOutputDirExists(cmdline_output_dir);
+    }
+
+    // create logFile
+    using boost::filesystem::path;
+    path fpath(deck_filename);
+    std::string baseName;
+    std::ostringstream debugFileStream;
+    std::ostringstream logFileStream;
+
+    // Strip extension "." or ".DATA"
+    std::string extension = boost::to_upper_copy(fpath.extension().string());
+    if (extension == ".DATA" || extension == ".") {
+        baseName = boost::to_upper_copy(fpath.stem().string());
+    } else {
+        baseName = boost::to_upper_copy(fpath.filename().string());
+    }
+
+    std::string output_dir = cmdline_output_dir;
+    if (output_dir.empty()) {
+        output_dir = absolute(path(baseName).parent_path()).string();
+    }
+
+    logFileStream << output_dir << "/" << baseName;
+    debugFileStream << output_dir << "/" << baseName;
+
+    if (mpi_rank_ != 0) {
+        // Added rank to log file for non-zero ranks.
+        // This prevents message loss.
+        debugFileStream << "." << mpi_rank_;
+        // If the following file appears then there is a bug.
+        logFileStream << "." << mpi_rank_;
+    }
+    logFileStream << ".PRT";
+    debugFileStream << ".DBG";
+
+    FileOutputMode output;
+    {
+        static std::map<std::string, FileOutputMode> stringToOutputMode =
+            { {"none", FileOutputMode::OUTPUT_NONE },
+              {"false", FileOutputMode::OUTPUT_LOG_ONLY },
+              {"log", FileOutputMode::OUTPUT_LOG_ONLY },
+              {"all" , FileOutputMode::OUTPUT_ALL },
+              {"true" , FileOutputMode::OUTPUT_ALL }};
+        auto outputModeIt = stringToOutputMode.find(cmdline_output);
+        if (outputModeIt != stringToOutputMode.end()) {
+            output = outputModeIt->second;
+        }
+        else {
+            output = FileOutputMode::OUTPUT_ALL;
+            std::cerr << "Value " << cmdline_output <<
+                " is not a recognized output mode. Using \"all\" instead."
+                      << std::endl;
+        }
+    }
+
+    if (output > FileOutputMode::OUTPUT_NONE) {
+        std::shared_ptr<Opm::EclipsePRTLog> prtLog = std::make_shared<Opm::EclipsePRTLog>(logFileStream.str(), Opm::Log::NoDebugMessageTypes, false, output_cout_);
+        Opm::OpmLog::addBackend("ECLIPSEPRTLOG", prtLog);
+        prtLog->setMessageLimiter(std::make_shared<Opm::MessageLimiter>());
+        prtLog->setMessageFormatter(std::make_shared<Opm::SimpleMessageFormatter>(false));
+    }
+
+    if (output >= FileOutputMode::OUTPUT_LOG_ONLY) {
+        std::string debugFile = debugFileStream.str();
+        std::shared_ptr<Opm::StreamLog> debugLog = std::make_shared<Opm::EclipsePRTLog>(debugFileStream.str(), Opm::Log::DefaultMessageTypes, false, output_cout_);
+        Opm::OpmLog::addBackend("DEBUGLOG", debugLog);
+    }
+
+    std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(std::cout, Opm::Log::StdoutMessageTypes);
+    Opm::OpmLog::addBackend(stdout_log_id, streamLog);
+    streamLog->setMessageFormatter(std::make_shared<Opm::SimpleMessageFormatter>(true));
+
+    return output;
+}
+
+
+
+void setupMessageLimiter(const Opm::MessageLimits msgLimits,  const std::string& stdout_log_id) {
+    std::shared_ptr<Opm::StreamLog> stream_log = Opm::OpmLog::getBackend<Opm::StreamLog>(stdout_log_id);
+
+    const std::map<int64_t, int> limits = {{Opm::Log::MessageType::Note,
+                                            msgLimits.getCommentPrintLimit(0)},
+                                           {Opm::Log::MessageType::Info,
+                                            msgLimits.getMessagePrintLimit(0)},
+                                           {Opm::Log::MessageType::Warning,
+                                            msgLimits.getWarningPrintLimit(0)},
+                                           {Opm::Log::MessageType::Error,
+                                            msgLimits.getErrorPrintLimit(0)},
+                                           {Opm::Log::MessageType::Problem,
+                                            msgLimits.getProblemPrintLimit(0)},
+                                           {Opm::Log::MessageType::Bug,
+                                            msgLimits.getBugPrintLimit(0)}};
+    stream_log->setMessageLimiter(std::make_shared<Opm::MessageLimiter>(10, limits));
 }
 
 
 // ----------------- Main program -----------------
 template<class TypeTag>
-int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles)
+int mainFlow(int argc, char** argv)
 {
+    Dune::Timer externalSetupTimer;
+    externalSetupTimer.start();
+
+    detail::handleVersionCmdLine(argc, argv);
     // MPI setup.
 #if HAVE_DUNE_FEM
     Dune::Fem::MPIManager::initialize(argc, argv);
@@ -151,15 +305,21 @@ int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles)
     typedef TTAG(FlowEarlyBird) PreTypeTag;
     typedef GET_PROP_TYPE(PreTypeTag, Problem) PreProblem;
 
-    PreProblem::setBriefDescription("Simple Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project.");
+    PreProblem::setBriefDescription("Flow, an advanced reservoir simulator for ECL-decks provided by the Open Porous Media project.");
-    
+
     int status = Opm::FlowMainEbos<PreTypeTag>::setupParameters_(argc, argv);
-    if (status != 0)
+    if (status != 0) {
         // if setupParameters_ returns a value smaller than 0, there was no error, but
         // the program should abort. This is the case e.g. for the --help and the
         // --print-properties parameters.
+#if HAVE_MPI
+        MPI_Finalize();
+#endif
         return (status >= 0)?status:0;
+    }
 
+    FileOutputMode outputMode = FileOutputMode::OUTPUT_NONE;
+    bool outputCout = false;
     if (mpiRank == 0)
         outputCout = EWOMS_GET_PARAM(PreTypeTag, bool, EnableTerminalOutput);
 
@@ -176,38 +336,67 @@ int mainFlow(int argc, char** argv, bool outputCout, bool outputFiles)
 
     // Create Deck and EclipseState.
     try {
-      Opm::Parser parser;
+        if (outputCout) {
-      typedef std::pair<std::string, Opm::InputError::Action> ParseModePair;
+            std::cout << "Reading deck file '" << deckFilename << "'\n";
-      typedef std::vector<ParseModePair> ParseModePairs;
+            std::cout.flush();
-      ParseModePairs tmp;
+        }
-      tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE));
+        std::shared_ptr<Opm::Deck> deck;
-      tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN));
+        std::shared_ptr<Opm::EclipseState> eclipseState;
-      tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN));
+        std::shared_ptr<Opm::Schedule> schedule;
-      tmp.push_back(ParseModePair(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN));
+        std::shared_ptr<Opm::SummaryConfig> summaryConfig;
-      Opm::ParseContext parseContext(tmp);
+        {
-      Opm::ErrorGuard errorGuard;
+            Opm::Parser parser;
+            Opm::ParseContext parseContext;
+            Opm::ErrorGuard errorGuard;
+            outputMode = setupLogging(mpiRank,
+                                      deckFilename,
+                                      EWOMS_GET_PARAM(PreTypeTag, std::string, OutputDir),
+                                      EWOMS_GET_PARAM(PreTypeTag, std::string, OutputMode),
+                                      outputCout, "STDOUT_LOGGER");
 
-      std::shared_ptr<Opm::Deck> deck = std::make_shared< Opm::Deck >( parser.parseFile(deckFilename , parseContext, errorGuard) );
+            if (EWOMS_GET_PARAM(PreTypeTag, bool, EclStrictParsing))
-      if ( outputCout ) {
+                parseContext.update( Opm::InputError::DELAYED_EXIT1);
-          Opm::checkDeck(*deck, parser, parseContext, errorGuard);
+            else {
-	Opm::MissingFeatures::checkKeywords(*deck);
+                parseContext.update(Opm::ParseContext::PARSE_RANDOM_SLASH, Opm::InputError::IGNORE);
-      }
+                parseContext.update(Opm::ParseContext::PARSE_MISSING_DIMS_KEYWORD, Opm::InputError::WARN);
-      //Opm::Runspec runspec( *deck );
+                parseContext.update(Opm::ParseContext::SUMMARY_UNKNOWN_WELL, Opm::InputError::WARN);
-      //const auto& phases = runspec.phases();
+                parseContext.update(Opm::ParseContext::SUMMARY_UNKNOWN_GROUP, Opm::InputError::WARN);
+            }
 
-      std::shared_ptr<Opm::EclipseState> eclipseState = std::make_shared< Opm::EclipseState > ( *deck, parseContext, errorGuard );
+            Opm::FlowMainEbos<PreTypeTag>::printPRTHeader(outputCout);
-      Opm::flowEbosSetDeck<TypeTag>(*deck, *eclipseState);
+
-      return Opm::flowEbosMain<TypeTag>(argc, argv, outputCout, outputFiles);
+            deck.reset( new Opm::Deck( parser.parseFile(deckFilename , parseContext, errorGuard)));
+            Opm::MissingFeatures::checkKeywords(*deck, parseContext, errorGuard);
+            if ( outputCout )
+                Opm::checkDeck(*deck, parser, parseContext, errorGuard);
+
+            eclipseState.reset( new Opm::EclipseState(*deck, parseContext, errorGuard ));
+            schedule.reset(new Opm::Schedule(*deck, *eclipseState, parseContext, errorGuard));
+            summaryConfig.reset( new Opm::SummaryConfig(*deck, *schedule, eclipseState->getTableManager(), parseContext, errorGuard));
+            setupMessageLimiter(schedule->getMessageLimits(), "STDOUT_LOGGER");
+
+            Opm::checkConsistentArrayDimensions(*eclipseState, *schedule, parseContext, errorGuard);
+
+            if (errorGuard) {
+                errorGuard.dump();
+                errorGuard.clear();
+
+                throw std::runtime_error("Unrecoverable errors were encountered while loading input.");
+            }
+        }
+        bool outputFiles = (outputMode != FileOutputMode::OUTPUT_NONE);
+        Opm::flowEbosSetDeck<TypeTag>(*deck, *eclipseState);
+        return Opm::flowEbosMain<TypeTag>(argc, argv, outputCout, outputFiles);
     }
   catch (const std::invalid_argument& e)
     {
       if (outputCout) {
-	std::cerr << "Failed to create valid EclipseState object." << std::endl;
+        std::cerr << "Failed to create valid EclipseState object." << std::endl;
-	std::cerr << "Exception caught: " << e.what() << std::endl;
+        std::cerr << "Exception caught: " << e.what() << std::endl;
       }
       throw;
     }
-    
+
     return EXIT_SUCCESS;
 }
 #endif

opm/core/props/phaseUsageFromDeck.hpp

@@ -54,14 +54,9 @@ namespace Opm
             }
         }
 
-        // Only 2 or 3 phase systems handled.
-        if (pu.num_phases < 2 || pu.num_phases > 3) {
-            OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
-        }
-
         // We need oil systems, since we do not support the keywords needed for
         // water-gas systems.
-        if (!pu.phase_used[BlackoilPhases::Liquid]) {
+        if (!pu.phase_used[BlackoilPhases::Liquid] && !(pu.num_phases == 1)) {
             OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
         }
 

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -422,28 +422,35 @@ namespace Opm {
 
                 Scalar saturationsOld[FluidSystem::numPhases] = { 0.0 };
                 Scalar oilSaturationOld = 1.0;
-                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
-                    oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
-                }
-
-                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
-                    saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
-                    oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
-                }
-
-                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-                    saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
-                }
 
+                // NB fix me! adding pressures changes to satutation changes does not make sense
                 Scalar tmp = pressureNew - pressureOld;
                 resultDelta += tmp*tmp;
                 resultDenom += pressureNew*pressureNew;
 
-                for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
+                if (FluidSystem::numActivePhases() > 1) {
-                    const Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
+                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    resultDelta += tmpSat * tmpSat;
+                        saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
-                    resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
+                        oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
+                    }
+
+                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
+                        priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
+                    {
+                        saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
+                        oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
+                    }
+
+                    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                        saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
+                    }
+                    for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) {
+                        Scalar tmpSat = saturationsNew[phaseIdx] - saturationsOld[phaseIdx];
+                        resultDelta += tmpSat*tmpSat;
+                        resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx];
+                        assert(std::isfinite(resultDelta));
+                        assert(std::isfinite(resultDenom));
+                    }
                 }
             }
 

opm/simulators/flow/FlowMainEbos.hpp

@@ -432,9 +432,11 @@ namespace Opm
                 return;
             }
 
-            // Run relperm diagnostics
+            // Run relperm diagnostics if we have more than one phase.
-            RelpermDiagnostics diagnostic;
+            if (FluidSystem::numActivePhases() > 1) {
-            diagnostic.diagnosis(eclState(), deck(), this->grid());
+                RelpermDiagnostics diagnostic;
+                diagnostic.diagnosis(eclState(), deck(), this->grid());
+            }
         }
 
         // Run the simulator.

opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp

@@ -60,6 +60,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
     {
         ++current_step_;
         current_time_ += dt_;
+        assert(dt_ > 0);
         // store used time step sizes
         steps_.push_back( dt_ );
         return *this;
@@ -71,7 +72,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
         double remaining = (total_time_ - current_time_);
         // apply max time step if it was set
         dt_ = std::min( dt_estimate, max_time_step_ );
-
+        assert(dt_ > 0);
         if( remaining > 0 ) {
 
             // set new time step (depending on remaining time)
@@ -81,6 +82,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
                 if( dt_ > max_time_step_ ) {
                     dt_ = 0.5 * remaining;
                 }
+                assert(dt_ > 0);
                 return;
             }
 
@@ -89,6 +91,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
 
             if( 1.5 * dt_ > remaining ) {
                 dt_ = 0.5 * remaining;
+                assert(dt_ > 0);
                 return;
             }
         }
@@ -102,6 +105,7 @@ AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ ()
 
     double AdaptiveSimulatorTimer::currentStepLength () const
     {
+      assert(dt_ > 0);
         return dt_;
     }
 

opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp

@@ -293,9 +293,10 @@ namespace Opm {
                     double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
                                                                                substepTimer.simulationTimeElapsed());
 
+                    assert(dtEstimate > 0);
                     // limit the growth of the timestep size by the growth factor
                     dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
-
+                    assert(dtEstimate > 0);
                     // further restrict time step size growth after convergence problems
                     if (restarts > 0) {
                         dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
@@ -309,7 +310,7 @@ namespace Opm {
                         const double maxPredictionTHPTimestep = 16.0 * unit::day;
                         dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
                     }
-
+                    assert(dtEstimate > 0);
                     if (timestepVerbose_) {
                         std::ostringstream ss;
                         substepReport.reportStep(ss);

opm/simulators/timestepping/TimeStepControl.cpp

@@ -128,12 +128,16 @@ namespace Opm
     {
         // shift errors
         for( int i=0; i<2; ++i ) {
-          errors_[ i ] = errors_[i+1];
+            errors_[ i ] = errors_[i+1];
         }
 
         // store new error
         const double error = relChange.relativeChange();
         errors_[ 2 ] = error;
+        for( int i=0; i<2; ++i ) {
+            assert(std::isfinite(errors_[i]));
+            assert(errors_[i]>0);
+        }
 
         if( error > tol_ )
         {

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -251,9 +251,19 @@ namespace Opm {
             const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
             const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
             const auto& fs = intQuants.fluidState();
-
+            // copy of get perfpressure in Standard well
-            const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
+            // exept for value
-            cellPressures[cellIdx] = p;
+            double perf_pressure = 0.0;
+            if (Indices::oilEnabled) {
+                perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
+            } else {
+                if (Indices::waterEnabled) {
+                    perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
+                } else {
+                    perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
+                }
+            }
+            cellPressures[cellIdx] = perf_pressure;
         }
         well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_);
 
@@ -1525,8 +1535,8 @@ namespace Opm {
     int
     BlackoilWellModel<TypeTag>::numComponents() const
     {
-        if (numPhases() == 2) {
+        if (numWells() > 0 && numPhases() < 3) {
-            return 2;
+            return numPhases();
         }
         int numComp = FluidSystem::numComponents;
         if (has_solvent_) {

opm/simulators/wells/RateConverter.hpp

@@ -493,11 +493,13 @@ namespace Opm {
 
                     // sum p, rs, rv, and T.
                     double hydrocarbonPV = pv_cell*hydrocarbon;
-                    pv += hydrocarbonPV;
+                    if (hydrocarbonPV > 0) {
-                    p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
+                        pv += hydrocarbonPV;
-                    rs += fs.Rs().value()*hydrocarbonPV;
+                        p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
-                    rv += fs.Rv().value()*hydrocarbonPV;
+                        rs += fs.Rs().value()*hydrocarbonPV;
-                    T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
+                        rv += fs.Rv().value()*hydrocarbonPV;
+                        T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
+                    }
                 }
 
                 for (const auto& reg : rmap_.activeRegions()) {
@@ -521,7 +523,7 @@ namespace Opm {
                 }
             }
 
-	    /**
+            /**
              * Region identifier.
              *
              * Integral type.

opm/simulators/wells/StandardWell.hpp

@@ -59,6 +59,7 @@ namespace Opm
         using typename Base::Indices;
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
+        using typename Base::FluidState;
 
         using Base::numEq;
 
@@ -293,6 +294,8 @@ namespace Opm
 
         EvalWell extendEval(const Eval& in) const;
 
+        Eval getPerfCellPressure(const FluidState& fs) const;
+
         // xw = inv(D)*(rw - C*x)
         void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -238,6 +238,10 @@ namespace Opm
     StandardWell<TypeTag>::
     wellVolumeFraction(const unsigned compIdx) const
     {
+        if (FluidSystem::numActivePhases() == 1) {
+            return EvalWell(numWellEq_ + numEq, 1.0);
+        }
+
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
             return primary_variables_evaluation_[WFrac];
         }
@@ -305,6 +309,24 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    typename StandardWell<TypeTag>::Eval
+    StandardWell<TypeTag>::getPerfCellPressure(const typename StandardWell<TypeTag>::FluidState& fs) const
+    {
+        Eval pressure;
+        if (Indices::oilEnabled) {
+            pressure = fs.pressure(FluidSystem::oilPhaseIdx);
+        } else {
+            if (Indices::waterEnabled) {
+                pressure = fs.pressure(FluidSystem::waterPhaseIdx);
+            } else {
+                pressure = fs.pressure(FluidSystem::gasPhaseIdx);
+            }
+        }
+        return pressure;
+    }
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
@@ -321,7 +343,7 @@ namespace Opm
     {
 
         const auto& fs = intQuants.fluidState();
-        const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
+        const EvalWell pressure = extendEval(getPerfCellPressure(fs));
         const EvalWell rs = extendEval(fs.Rs());
         const EvalWell rv = extendEval(fs.Rv());
         std::vector<EvalWell> b_perfcells_dense(num_components_, EvalWell{numWellEq_ + numEq, 0.0});
@@ -669,7 +691,9 @@ namespace Opm
                         || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) {
 
                     const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
-                    const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value();
+                    const auto& fs = intQuants.fluidState();
+                    Eval perf_pressure = getPerfCellPressure(fs);
+                    const double drawdown  = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value();
                     const bool new_well = schedule.hasWellEvent(name(), ScheduleEvents::NEW_WELL, current_step_);
                     double productivity_index = cq_s[compIdx].value() / drawdown;
                     scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger);
@@ -679,12 +703,15 @@ namespace Opm
 
         }
 
-
         // add vol * dF/dt + Q to the well equations;
         for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
             // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
             // since all the rates are under surface condition
-            EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
+            EvalWell resWell_loc(numWellEq_ + numEq, 0.0);
+            if (FluidSystem::numActivePhases() > 1) {
+                assert(dt > 0);
+                resWell_loc += (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
+            }
             resWell_loc -= getQs(componentIdx) * well_efficiency_factor_;
             for (int pvIdx = 0; pvIdx < numWellEq_; ++pvIdx) {
                 invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
@@ -919,29 +946,30 @@ namespace Opm
         const double relaxation_factor_fractions = (well_type_ == PRODUCER) ?
                                          relaxationFactorFractionsProducer(old_primary_variables, dwells)
                                        : 1.0;
+        if (FluidSystem::numActivePhases() > 1) {
+            // update the second and third well variable (The flux fractions)
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
+                const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
+                // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+                primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+            }
 
-        // update the second and third well variable (The flux fractions)
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
-            const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
+                const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
-            const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
+                primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
-            // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+            }
-            primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+
+            if (has_solvent) {
+                const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
+                const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
+                primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
+            }
+
+            processFractions();
         }
 
-        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
-            const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
-            primary_variables_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
-        }
-
-        if (has_solvent) {
-            const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
-            const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
-            primary_variables_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
-        }
-
-        processFractions();
-
         // updating the total rates Q_t
         const double relaxation_factor_rate = relaxationFactorRate(old_primary_variables, dwells);
         primary_variables_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
@@ -956,6 +984,10 @@ namespace Opm
         }
 
         updateExtraPrimaryVariables(dwells);
+
+        for (double v : primary_variables_) {
+            assert(Opm::isfinite(v));
+        }
     }
 
 
@@ -1073,49 +1105,53 @@ namespace Opm
     updateWellStateFromPrimaryVariables(WellState& well_state, Opm::DeferredLogger& deferred_logger) const
     {
         const PhaseUsage& pu = phaseUsage();
-        assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
-        const int oil_pos = pu.phase_pos[Oil];
 
         std::vector<double> F(number_of_phases_, 0.0);
-        F[oil_pos] = 1.0;
+        if (FluidSystem::numActivePhases() == 1) {
+            F[0] = 1.0;
+        } else {
+            assert( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) );
+            const int oil_pos = pu.phase_pos[Oil];
+            F[oil_pos] = 1.0;
 
-        if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
+            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) {
-            const int water_pos = pu.phase_pos[Water];
+                const int water_pos = pu.phase_pos[Water];
-            F[water_pos] = primary_variables_[WFrac];
+                F[water_pos] = primary_variables_[WFrac];
-            F[oil_pos] -= F[water_pos];
+                F[oil_pos] -= F[water_pos];
-        }
-
-        if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
-            const int gas_pos = pu.phase_pos[Gas];
-            F[gas_pos] = primary_variables_[GFrac];
-            F[oil_pos] -= F[gas_pos];
-        }
-
-        double F_solvent = 0.0;
-        if (has_solvent) {
-            F_solvent = primary_variables_[SFrac];
-            F[oil_pos] -= F_solvent;
-        }
-
-        // convert the fractions to be Q_p / G_total to calculate the phase rates
-        for (int p = 0; p < number_of_phases_; ++p) {
-            const double scal = scalingFactor(p);
-            // for injection wells, there should only one non-zero scaling factor
-            if (scal > 0) {
-                F[p] /= scal ;
-            } else {
-                // this should only happens to injection wells
-                F[p] = 0.;
             }
-        }
 
-        // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
+            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) {
-        // More testing is needed to make sure this is correct for well groups and THP.
+                const int gas_pos = pu.phase_pos[Gas];
-        if (has_solvent){
+                F[gas_pos] = primary_variables_[GFrac];
-            F_solvent /= scalingFactor(contiSolventEqIdx);
+                F[oil_pos] -= F[gas_pos];
-            F[pu.phase_pos[Gas]] += F_solvent;
+            }
-        }
 
+            double F_solvent = 0.0;
+            if (has_solvent) {
+                F_solvent = primary_variables_[SFrac];
+                F[oil_pos] -= F_solvent;
+            }
+
+            // convert the fractions to be Q_p / G_total to calculate the phase rates
+            for (int p = 0; p < number_of_phases_; ++p) {
+                const double scal = scalingFactor(p);
+                // for injection wells, there should only one non-zero scaling factor
+                if (scal > 0) {
+                    F[p] /= scal ;
+                } else {
+                    // this should only happens to injection wells
+                    F[p] = 0.;
+                }
+            }
+
+            // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
+            // More testing is needed to make sure this is correct for well groups and THP.
+            if (has_solvent){
+                F_solvent /= scalingFactor(contiSolventEqIdx);
+                F[pu.phase_pos[Gas]] += F_solvent;
+            }
+
+        }
         well_state.bhp()[index_of_well_] = primary_variables_[Bhp];
 
         // calculate the phase rates based on the primary variables
@@ -1292,6 +1328,10 @@ namespace Opm
 
             break;
         } // end of switch
+
+        for (double v : primary_variables_) {
+            assert(Opm::isfinite(v));
+        }
     }
 
 
@@ -1324,7 +1364,8 @@ namespace Opm
             const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
             const auto& fs = int_quantities.fluidState();
             // the pressure of the reservoir grid block the well connection is in
-            const double p_r = fs.pressure(FluidSystem::oilPhaseIdx).value();
+            Eval perf_pressure = getPerfCellPressure(fs);
+            double p_r = perf_pressure.value();
 
             // calculating the b for the connection
             std::vector<double> b_perf(num_components_);
@@ -1550,8 +1591,9 @@ namespace Opm
             const int cell_idx = well_cells_[perf];
             const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
             const auto& fs = intQuants.fluidState();
+            Eval perf_pressure = getPerfCellPressure(fs);
+            const double pressure = perf_pressure.value();
 
-            const double pressure = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
             const double bhp = getBhp().value();
 
             // Pressure drawdown (also used to determine direction of flow)
@@ -1755,6 +1797,7 @@ namespace Opm
             // TODO: to check why should be perf - 1
             const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1];
             const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
+            // strange choice
             const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
 
             if (waterPresent) {
@@ -2518,57 +2561,57 @@ namespace Opm
             }
         }
 
-        if (std::abs(total_well_rate) > 0.) {
+        if (FluidSystem::numActivePhases() > 1) {
-            const auto pu = phaseUsage();
+            if (std::abs(total_well_rate) > 0.) {
-            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                const auto pu = phaseUsage();
-                primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
-            }
-            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
-            }
-            if (has_solvent) {
-                primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
-            }
-        } else { // total_well_rate == 0
-            if (well_type_ == INJECTOR) {
-                auto phase = well_ecl_.getInjectionProperties().injectorType;
-                // only single phase injection handled
                 if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    if (phase == Well2::InjectorType::WATER) {
+                    primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
-                        primary_variables_[WFrac] = 1.0;
-                    } else {
-                        primary_variables_[WFrac] = 0.0;
-                    }
                 }
-
                 if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    if (phase == Well2::InjectorType::GAS) {
+                    primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / total_well_rate ;
-                        primary_variables_[GFrac] = 1.0 - wsolvent();
+                }
-                        if (has_solvent) {
+                if (has_solvent) {
-                            primary_variables_[SFrac] = wsolvent();
+                    primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
+                }
+            } else { // total_well_rate == 0
+                if (well_type_ == INJECTOR) {
+                    auto phase = well_ecl_.getInjectionProperties().injectorType;
+                    // only single phase injection handled
+                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                        if (phase == Well2::InjectorType::WATER) {
+                            primary_variables_[WFrac] = 1.0;
+                        } else {
+                            primary_variables_[WFrac] = 0.0;
+                        }
+                    }
+                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                        if (phase == Well2::InjectorType::GAS) {
+                            primary_variables_[GFrac] = 1.0 - wsolvent();
+                            if (has_solvent) {
+                                primary_variables_[SFrac] = wsolvent();
+                            }
+                        } else {
+                            primary_variables_[GFrac] = 0.0;
                         }
-                    } else {
-                        primary_variables_[GFrac] = 0.0;
                     }
-                }
 
-                // TODO: it is possible to leave injector as a oil well,
+                    // TODO: it is possible to leave injector as a oil well,
-                // when F_w and F_g both equals to zero, not sure under what kind of circumstance
+                    // when F_w and F_g both equals to zero, not sure under what kind of circumstance
-                // this will happen.
+                    // this will happen.
-            } else if (well_type_ == PRODUCER) { // producers
+                } else if (well_type_ == PRODUCER) { // producers
-                // TODO: the following are not addressed for the solvent case yet
+                    // TODO: the following are not addressed for the solvent case yet
-                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    primary_variables_[WFrac] = 1.0 / np;
+                        primary_variables_[WFrac] = 1.0 / np;
+                    }
+                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                        primary_variables_[GFrac] = 1.0 / np;
+                    }
+                } else {
+                    OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
                 }
-                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    primary_variables_[GFrac] = 1.0 / np;
-                }
-            } else {
-                OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
             }
         }
 
-
         // BHP
         primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
 
@@ -2579,6 +2622,10 @@ namespace Opm
                 primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
             }
         }
+
+        for (double v : primary_variables_) {
+            assert(Opm::isfinite(v));
+        }
     }
 
 
@@ -2824,28 +2871,31 @@ namespace Opm
         // 0.95 is a experimental value, which remains to be optimized
         double relaxation_factor = 1.0;
 
-        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        if (FluidSystem::numActivePhases() > 1) {
-            const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-            relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
+                const double relaxation_factor_w = relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
-        }
+                relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
+            }
 
-        if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
+                const double relaxation_factor_g = relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
-            relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
+                relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
-        }
+            }
 
-        if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
-            // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
+                && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-            // not be negative oil fraction later
+                // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so
-            const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
+                // there will not be negative oil fraction later
-            const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
+                const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
-            const double possible_updated_sum = original_sum - relaxed_update;
+                const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
+                const double possible_updated_sum = original_sum - relaxed_update;
 
-            if (possible_updated_sum > 1.0) {
+                if (possible_updated_sum > 1.0) {
-                assert(relaxed_update != 0.);
+                    assert(relaxed_update != 0.);
 
-                const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
+                    const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
-                relaxation_factor *= further_relaxation_factor;
+                    relaxation_factor *= further_relaxation_factor;
+                }
             }
         }
 

opm/simulators/wells/WellInterface.hpp

@@ -94,6 +94,7 @@ namespace Opm
         static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
         static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
         static const bool has_energy = GET_PROP_VALUE(TypeTag, EnableEnergy);
+        static const bool has_temperature = GET_PROP_VALUE(TypeTag, EnableTemperature);
         // flag for polymer molecular weight related
         static const bool has_polymermw = GET_PROP_VALUE(TypeTag, EnablePolymerMW);
         static const bool has_foam = GET_PROP_VALUE(TypeTag, EnableFoam);
@@ -106,7 +107,13 @@ namespace Opm
         // For the conversion between the surface volume rate and reservoir voidage rate
         using RateConverterType = RateConverter::
         SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
-
+        static const bool compositionSwitchEnabled = Indices::gasEnabled;
+        using FluidState = Opm::BlackOilFluidState<Eval,
+                                                   FluidSystem,
+                                                   has_temperature,
+                                                   has_energy,
+                                                   compositionSwitchEnabled,
+                                                   Indices::numPhases >;
         /// Constructor
         WellInterface(const Well2& well, const int time_step, const Wells* wells,
                       const ModelParameters& param,