Merge pull request #800 from OPM/frankenstein

Frankenstein V2

CMakeLists.txt

@@ -81,6 +81,7 @@ include (CMakeLists_files.cmake)
 macro (config_hook)
 	opm_need_version_of ("dune-common")
 	opm_need_version_of ("dune-istl")
+	opm_need_version_of ("ewoms")
 endmacro (config_hook)
 
 macro (prereqs_hook)
@@ -148,10 +149,7 @@ if (HAVE_OPM_DATA)
  endif()
 
  # create a symbolic link from flow to flow_legacy
+message("create symlink")
 ADD_CUSTOM_TARGET(flow ALL
-	COMMAND ${CMAKE_COMMAND} -E create_symlink "flow_legacy" "${CMAKE_BINARY_DIR}/bin/flow")
+	                COMMAND ${CMAKE_COMMAND} -E create_symlink "flow_legacy" "${CMAKE_BINARY_DIR}/bin/flow")
-install(
-	FILES "${CMAKE_BINARY_DIR}/bin/flow"
-	DESTINATION "${CMAKE_INSTALL_PREFIX}/bin"
-	)
 

CMakeLists_files.cmake

@@ -36,6 +36,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/autodiff/ImpesTPFAAD.cpp
   opm/autodiff/moduleVersion.cpp
   opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
+  opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp
   opm/autodiff/SimulatorIncompTwophaseAd.cpp
   opm/autodiff/TransportSolverTwophaseAd.cpp
   opm/autodiff/BlackoilPropsAdFromDeck.cpp
@@ -49,7 +50,6 @@ list (APPEND MAIN_SOURCE_FILES
   opm/autodiff/WellMultiSegment.cpp
 	opm/autodiff/MultisegmentWells.cpp
   opm/autodiff/BlackoilSolventState.cpp
-  opm/autodiff/ThreadHandle.hpp
   opm/autodiff/MissingFeatures.cpp
   opm/polymer/PolymerState.cpp
   opm/polymer/PolymerBlackoilState.cpp
@@ -105,6 +105,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
   examples/find_zero.cpp
   examples/flow_legacy.cpp
   examples/flow_sequential.cpp
+  examples/flow_ebos.cpp
   examples/flow_multisegment.cpp
   examples/flow_solvent.cpp
   examples/sim_2p_incomp.cpp
@@ -124,6 +125,7 @@ list (APPEND PROGRAM_SOURCE_FILES
   examples/sim_2p_incomp.cpp
   examples/sim_2p_incomp_ad.cpp
   examples/sim_2p_comp_reorder.cpp
+  examples/flow_ebos.cpp
   examples/flow_legacy.cpp
   examples/flow_sequential.cpp
   examples/flow_solvent.cpp
@@ -143,6 +145,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/AutoDiffMatrix.hpp
   opm/autodiff/AutoDiff.hpp
   opm/autodiff/BackupRestore.hpp
+	opm/autodiff/BlackoilDetails.hpp
   opm/autodiff/BlackoilModel.hpp
   opm/autodiff/BlackoilModelBase.hpp
   opm/autodiff/BlackoilModelBase_impl.hpp
@@ -155,6 +158,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/Compat.hpp
   opm/autodiff/CPRPreconditioner.hpp
   opm/autodiff/createGlobalCellArray.hpp
+	opm/autodiff/DefaultBlackoilSolutionState.hpp
   opm/autodiff/BlackoilSequentialModel.hpp
   opm/autodiff/BlackoilSolventModel.hpp
   opm/autodiff/BlackoilSolventModel_impl.hpp
@@ -166,6 +170,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/DuneMatrix.hpp
   opm/autodiff/ExtractParallelGridInformationToISTL.hpp
   opm/autodiff/FlowMain.hpp
+  opm/autodiff/FlowMainEbos.hpp
   opm/autodiff/FlowMainPolymer.hpp
   opm/autodiff/FlowMainSequential.hpp
   opm/autodiff/FlowMainSolvent.hpp
@@ -173,6 +178,8 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/GridHelpers.hpp
   opm/autodiff/GridInit.hpp
   opm/autodiff/ImpesTPFAAD.hpp
+  opm/autodiff/ISTLSolver.hpp
+	opm/autodiff/IterationReport.hpp
   opm/autodiff/moduleVersion.hpp
   opm/autodiff/NewtonIterationBlackoilCPR.hpp
   opm/autodiff/NewtonIterationBlackoilInterface.hpp
@@ -187,18 +194,22 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp
   opm/autodiff/RateConverter.hpp
   opm/autodiff/RedistributeDataHandles.hpp
+	opm/autodiff/SimFIBODetails.hpp
   opm/autodiff/SimulatorBase.hpp
   opm/autodiff/SimulatorBase_impl.hpp
+	opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
+	opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp
   opm/autodiff/SimulatorIncompTwophaseAd.hpp
   opm/autodiff/SimulatorSequentialBlackoil.hpp
   opm/autodiff/TransportSolverTwophaseAd.hpp
   opm/autodiff/WellDensitySegmented.hpp
   opm/autodiff/WellStateFullyImplicitBlackoil.hpp
+  opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
   opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp
   opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
   opm/autodiff/VFPProperties.hpp
@@ -212,9 +223,11 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/WellHelpers.hpp
   opm/autodiff/StandardWells.hpp
   opm/autodiff/StandardWells_impl.hpp
+	opm/autodiff/StandardWellsDense.hpp
   opm/autodiff/StandardWellsSolvent.hpp
   opm/autodiff/StandardWellsSolvent_impl.hpp
   opm/autodiff/MissingFeatures.hpp
+  opm/autodiff/ThreadHandle.hpp
   opm/polymer/CompressibleTpfaPolymer.hpp
   opm/polymer/GravityColumnSolverPolymer.hpp
   opm/polymer/GravityColumnSolverPolymer_impl.hpp
@@ -243,6 +256,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/ParallelFileMerger.hpp
   opm/simulators/SimulatorCompressibleTwophase.hpp
   opm/simulators/SimulatorIncompTwophase.hpp
+	opm/simulators/thresholdPressures.hpp
 	opm/simulators/WellSwitchingLogger.hpp
   )
 

dune.module

@@ -10,4 +10,4 @@ Label: 2017.04-pre
 Maintainer: atgeirr@sintef.no
 MaintainerName: Atgeirr F. Rasmussen
 Url: http://opm-project.org
-Depends: opm-common opm-parser opm-output opm-material opm-core opm-grid dune-istl (>=2.2)
+Depends: opm-common opm-parser opm-output opm-material opm-core opm-grid dune-istl (>=2.2) ewoms

examples/flow_ebos.cpp

@@ -0,0 +1,39 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 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/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
+#include <dune/grid/CpGrid.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
+#include <opm/autodiff/FlowMainEbos.hpp>
+
+
+// ----------------- Main program -----------------
+int main(int argc, char** argv)
+{
+    Opm::FlowMainEbos mainfunc;
+    return mainfunc.execute(argc, argv);
+}

jenkins/build.sh

@@ -7,7 +7,8 @@ upstreams=(opm-common
            opm-output
            opm-material
            opm-core
-           opm-grid)
+           opm-grid
+           ewoms)
 
 declare -A upstreamRev
 upstreamRev[opm-common]=master
@@ -17,6 +18,7 @@ upstreamRev[opm-material]=master
 upstreamRev[opm-core]=master
 upstreamRev[opm-grid]=master
 upstreamRev[opm-output]=master
+upstreamRev[ewoms]=master
 
 if grep -q "opm-common=" <<< $ghprbCommentBody
 then

jenkins/run-norne.sh

@@ -2,10 +2,11 @@
 
 pushd .
 cd deps/opm-data
+test -z $SIM && SIM=flow
 
 # Run the norne case
 cd norne
-$WORKSPACE/$configuration/build-opm-simulators/bin/flow deck_filename=NORNE_ATW2013.DATA output_dir=OPM
+$WORKSPACE/$configuration/build-opm-simulators/bin/$SIM deck_filename=NORNE_ATW2013.DATA output_dir=OPM
 test $? -eq 0 || exit 1
 PATH=$WORKSPACE/$configuration/install/bin:$PATH ./plotwells.sh
 

opm/autodiff/BlackoilDetails.hpp

@@ -0,0 +1,363 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2015 NTNU
+  Copyright 2015 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/>.
+*/
+
+#ifndef OPM_BLACKOILDETAILS_HEADER_INCLUDED
+#define OPM_BLACKOILDETAILS_HEADER_INCLUDED
+
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
+
+#include <Eigen/Eigen>
+#include <Eigen/Sparse>
+
+namespace Opm {
+namespace detail {
+
+
+    inline
+    std::vector<int>
+    buildAllCells(const int nc)
+    {
+        std::vector<int> all_cells(nc);
+
+        for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
+
+        return all_cells;
+    }
+
+
+
+    template <class PU>
+    std::vector<bool>
+    activePhases(const PU& pu)
+    {
+        const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+        std::vector<bool> active(maxnp, false);
+
+        for (int p = 0; p < pu.MaxNumPhases; ++p) {
+            active[ p ] = pu.phase_used[ p ] != 0;
+        }
+
+        return active;
+    }
+
+
+
+    template <class PU>
+    std::vector<int>
+    active2Canonical(const PU& pu)
+    {
+        const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
+        std::vector<int> act2can(maxnp, -1);
+
+        for (int phase = 0; phase < maxnp; ++phase) {
+            if (pu.phase_used[ phase ]) {
+                act2can[ pu.phase_pos[ phase ] ] = phase;
+            }
+        }
+
+        return act2can;
+    }
+
+
+
+    inline
+    double getGravity(const double* g, const int dim) {
+        double grav = 0.0;
+        if (g) {
+            // Guard against gravity in anything but last dimension.
+            for (int dd = 0; dd < dim - 1; ++dd) {
+                assert(g[dd] == 0.0);
+            }
+            grav = g[dim - 1];
+        }
+        return grav;
+    }
+
+        /// \brief Compute the L-infinity norm of a vector
+        /// \warn This function is not suitable to compute on the well equations.
+        /// \param a The container to compute the infinity norm on.
+        ///          It has to have one entry for each cell.
+        /// \param info In a parallel this holds the information about the data distribution.
+        template <class ADB>
+        inline
+        double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() )
+        {
+            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& real_info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                double result=0;
+                real_info.computeReduction(a.value(), Reduction::makeLInfinityNormFunctor<double>(), result);
+                return result;
+            }
+            else
+#endif
+            {
+                if( a.value().size() > 0 ) {
+                    return a.value().matrix().template lpNorm<Eigen::Infinity> ();
+                }
+                else { // this situation can occur when no wells are present
+                    return 0.0;
+                }
+            }
+        }
+
+        /// \brief Compute the Euclidian norm of a vector
+        /// \warning In the case that num_components is greater than 1
+        ///          an interleaved ordering is assumed. E.g. for each cell
+        ///          all phases of that cell are stored consecutively. First
+        ///          the ones for cell 0, then the ones for cell 1, ... .
+        /// \param it              begin iterator for the given vector
+        /// \param end             end iterator for the given vector
+        /// \param num_components  number of components (i.e. phases) in the vector
+        /// \param pinfo           In a parallel this holds the information about the data distribution.
+        template <class Iterator>
+        inline
+        double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() )
+        {
+            static_cast<void>(num_components); // Suppress warning in the serial case.
+            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                typedef typename Iterator::value_type Scalar;
+                Scalar product = 0.0;
+                int size_per_component = (end - it);
+                size_per_component /= num_components; // two lines to supresse unused warning.
+                assert((end - it) == num_components * size_per_component);
+
+                if( num_components == 1 )
+                {
+                    auto component_container =
+                        boost::make_iterator_range(it, end);
+                    info.computeReduction(component_container,
+                                           Opm::Reduction::makeInnerProductFunctor<double>(),
+                                           product);
+                }
+                else
+                {
+                    auto& maskContainer = info.getOwnerMask();
+                    auto mask = maskContainer.begin();
+                    assert(static_cast<int>(maskContainer.size()) == size_per_component);
+
+                    for(int cell = 0; cell < size_per_component; ++cell, ++mask)
+                    {
+                        Scalar cell_product = (*it) * (*it);
+                        ++it;
+                        for(int component=1; component < num_components;
+                            ++component, ++it)
+                        {
+                            cell_product += (*it) * (*it);
+                        }
+                        product += cell_product * (*mask);
+                    }
+                }
+                return info.communicator().sum(product);
+            }
+            else
+#endif
+            {
+                double product = 0.0 ;
+                for( ; it != end; ++it ) {
+                    product += ( *it * *it );
+                }
+                return product;
+            }
+        }
+        /// \brief Compute the reduction within the convergence check.
+        /// \param[in] B     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. tempV.col(i) contains the
+        ///                   values
+        ///                  for phase i.
+        /// \param[in] R     A matrix with MaxNumPhases columns and the same number rows
+        ///                  as the number of cells of the grid. B.col(i) contains the values
+        ///                  for phase i.
+        /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
+        ///                   of R for the phase i.
+        /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
+        ///                   maximum of tempV for the phase i.
+        /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
+        ///                   of B for the phase i.
+        /// \param[out] maxNormWell The maximum of the well flux equations for each phase.
+        /// \param[in]  nc    The number of cells of the local grid.
+        /// \return The total pore volume over all cells.
+        inline
+        double
+        convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
+                             const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
+                             const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
+                             std::vector<double>& R_sum,
+                             std::vector<double>& maxCoeff,
+                             std::vector<double>& B_avg,
+                             std::vector<double>& maxNormWell,
+                             int nc,
+                             int np,
+                             const std::vector<double> pv,
+                             std::vector<double> residual_well)
+        {
+            const int nw = residual_well.size() / np;
+            assert(nw * np == int(residual_well.size()));
+
+            // Do the global reductions
+#if 0 // HAVE_MPI
+            if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
+
+                // Compute the global number of cells and porevolume
+                std::vector<int> v(nc, 1);
+                auto nc_and_pv = std::tuple<int, double>(0, 0.0);
+                auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
+                                                           Opm::Reduction::makeGlobalSumFunctor<double>());
+                auto nc_and_pv_containers  = std::make_tuple(v, pv);
+                info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
+
+                for ( int idx = 0; idx < np; ++idx )
+                {
+                    auto values     = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
+                    auto containers = std::make_tuple(B.col(idx),
+                                                      tempV.col(idx),
+                                                      R.col(idx));
+                    auto operators  = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalMaxFunctor<double>(),
+                                                      Opm::Reduction::makeGlobalSumFunctor<double>());
+                    info.computeReduction(containers, operators, values);
+                    B_avg[idx]       = std::get<0>(values)/std::get<0>(nc_and_pv);
+                    maxCoeff[idx]    = std::get<1>(values);
+                    R_sum[idx]       = std::get<2>(values);
+                    assert(np >= np);
+                    if (idx < np) {
+                        maxNormWell[idx] = 0.0;
+                        for ( int w = 0; w < nw; ++w ) {
+                            maxNormWell[idx]  = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
+                        }
+                    }
+                }
+                info.communicator().max(maxNormWell.data(), np);
+                // Compute pore volume
+                return std::get<1>(nc_and_pv);
+            }
+            else
+#endif
+            {
+                B_avg.resize(np);
+                maxCoeff.resize(np);
+                R_sum.resize(np);
+                maxNormWell.resize(np);
+                for ( int idx = 0; idx < np; ++idx )
+                {
+                    B_avg[idx] = B.col(idx).sum()/nc;
+                    maxCoeff[idx] = tempV.col(idx).maxCoeff();
+                    R_sum[idx] = R.col(idx).sum();
+
+                    assert(np >= np);
+                    if (idx < np) {
+                        maxNormWell[idx] = 0.0;
+                        for ( int w = 0; w < nw; ++w ) {
+                            maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
+                        }
+                    }
+                }
+                // Compute total pore volume
+                return std::accumulate(pv.begin(), pv.end(), 0.0);
+            }
+        }
+
+
+        template <class Scalar>
+        inline
+        double
+        convergenceReduction(const std::vector< std::vector< Scalar > >& B,
+                             const std::vector< std::vector< Scalar > >& tempV,
+                             const std::vector< std::vector< Scalar > >& R,
+                             std::vector< Scalar >& R_sum,
+                             std::vector< Scalar >& maxCoeff,
+                             std::vector< Scalar >& B_avg,
+                             std::vector< Scalar >& maxNormWell,
+                             const int nc,
+                             const int np,
+                             const std::vector< Scalar >& pv,
+                             const std::vector< Scalar >& residual_well)
+        {
+            const int nw = residual_well.size() / np;
+            assert(nw * np == int(residual_well.size()));
+
+            // Do the global reductions
+            {
+                B_avg.resize(np);
+                maxCoeff.resize(np);
+                R_sum.resize(np);
+                maxNormWell.resize(np);
+                for ( int idx = 0; idx < np; ++idx )
+                {
+                    B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(), 0.0 ) / nc;
+                    R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(), 0.0 );
+                    maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() ));
+
+                    assert(np >= np);
+                    if (idx < np) {
+                        maxNormWell[idx] = 0.0;
+                        for ( int w = 0; w < nw; ++w ) {
+                            maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_well[nw*idx + w]));
+                        }
+                    }
+                }
+                // Compute total pore volume
+                return std::accumulate(pv.begin(), pv.end(), 0.0);
+            }
+        }
+
+        /// \brief Compute the L-infinity norm of a vector representing a well equation.
+        /// \param a The container to compute the infinity norm on.
+        /// \param info In a parallel this holds the information about the data distribution.
+        template <class ADB>
+        inline
+        double infinityNormWell( const ADB& a, const boost::any& pinfo )
+        {
+            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
+            double result=0;
+            if( a.value().size() > 0 ) {
+                result = a.value().matrix().template lpNorm<Eigen::Infinity> ();
+            }
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& real_info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                result = real_info.communicator().max(result);
+            }
+#endif
+            return result;
+        }
+    } // namespace detail
+} // namespace Opm
+
+#endif // OPM_BLACKOILDETAILS_HEADER_INCLUDED

opm/autodiff/BlackoilModelBase.hpp

@@ -33,6 +33,8 @@
 #include <opm/autodiff/BlackoilModelEnums.hpp>
 #include <opm/autodiff/VFPProperties.hpp>
 #include <opm/autodiff/RateConverter.hpp>
+#include <opm/autodiff/IterationReport.hpp>
+#include <opm/autodiff/DefaultBlackoilSolutionState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
 
@@ -49,47 +51,6 @@ namespace Opm {
     class VFPProperties;
     class SimulationDataContainer;
 
-    /// Struct for containing iteration variables.
-    struct DefaultBlackoilSolutionState
-    {
-        typedef AutoDiffBlock<double> ADB;
-        explicit DefaultBlackoilSolutionState(const int np)
-            : pressure  (    ADB::null())
-            , temperature(   ADB::null())
-            , saturation(np, ADB::null())
-            , rs        (    ADB::null())
-            , rv        (    ADB::null())
-            , qs        (    ADB::null())
-            , bhp       (    ADB::null())
-            , canonical_phase_pressures(3, ADB::null())
-        {
-        }
-        ADB              pressure;
-        ADB              temperature;
-        std::vector<ADB> saturation;
-        ADB              rs;
-        ADB              rv;
-        ADB              qs;
-        ADB              bhp;
-        // Below are quantities stored in the state for optimization purposes.
-        std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
-    };
-
-
-
-
-    /// Class used for reporting the outcome of a nonlinearIteration() call.
-    struct IterationReport
-    {
-        bool failed;
-        bool converged;
-        int linear_iterations;
-        int well_iterations;
-    };
-
-
-
-
     /// Traits to encapsulate the types used by classes using or
     /// extending this model. Forward declared here, must be
     /// specialised for each concrete model class.

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -25,6 +25,7 @@
 #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
 
 #include <opm/autodiff/BlackoilModelBase.hpp>
+#include <opm/autodiff/BlackoilDetails.hpp>
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
@@ -93,73 +94,6 @@ typedef Eigen::Array<double,
                      Eigen::Dynamic,
                      Eigen::RowMajor> DataBlock;
 
-
-namespace detail {
-
-
-    inline
-    std::vector<int>
-    buildAllCells(const int nc)
-    {
-        std::vector<int> all_cells(nc);
-
-        for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
-
-        return all_cells;
-    }
-
-
-
-    template <class PU>
-    std::vector<bool>
-    activePhases(const PU& pu)
-    {
-        const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
-        std::vector<bool> active(maxnp, false);
-
-        for (int p = 0; p < pu.MaxNumPhases; ++p) {
-            active[ p ] = pu.phase_used[ p ] != 0;
-        }
-
-        return active;
-    }
-
-
-
-    template <class PU>
-    std::vector<int>
-    active2Canonical(const PU& pu)
-    {
-        const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
-        std::vector<int> act2can(maxnp, -1);
-
-        for (int phase = 0; phase < maxnp; ++phase) {
-            if (pu.phase_used[ phase ]) {
-                act2can[ pu.phase_pos[ phase ] ] = phase;
-            }
-        }
-
-        return act2can;
-    }
-
-
-
-    inline
-    double getGravity(const double* g, const int dim) {
-        double grav = 0.0;
-        if (g) {
-            // Guard against gravity in anything but last dimension.
-            for (int dd = 0; dd < dim - 1; ++dd) {
-                assert(g[dd] == 0.0);
-            }
-            grav = g[dim - 1];
-        }
-        return grav;
-    }
-
-} // namespace detail
-
-
     template <class Grid, class WellModel, class Implementation>
     BlackoilModelBase<Grid, WellModel, Implementation>::
     BlackoilModelBase(const ModelParameters&          param,
@@ -1155,135 +1089,6 @@ namespace detail {
         return linsolver_.computeNewtonIncrement(residual_);
     }
 
-
-
-
-
-    namespace detail
-    {
-        /// \brief Compute the L-infinity norm of a vector
-        /// \warn This function is not suitable to compute on the well equations.
-        /// \param a The container to compute the infinity norm on.
-        ///          It has to have one entry for each cell.
-        /// \param info In a parallel this holds the information about the data distribution.
-        inline
-        double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() )
-        {
-            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
-#if HAVE_MPI
-            if ( pinfo.type() == typeid(ParallelISTLInformation) )
-            {
-                const ParallelISTLInformation& real_info =
-                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
-                double result=0;
-                real_info.computeReduction(a.value(), Reduction::makeLInfinityNormFunctor<double>(), result);
-                return result;
-            }
-            else
-#endif
-            {
-                if( a.value().size() > 0 ) {
-                    return a.value().matrix().lpNorm<Eigen::Infinity> ();
-                }
-                else { // this situation can occur when no wells are present
-                    return 0.0;
-                }
-            }
-        }
-
-        /// \brief Compute the Euclidian norm of a vector
-        /// \warning In the case that num_components is greater than 1
-        ///          an interleaved ordering is assumed. E.g. for each cell
-        ///          all phases of that cell are stored consecutively. First
-        ///          the ones for cell 0, then the ones for cell 1, ... .
-        /// \param it              begin iterator for the given vector
-        /// \param end             end iterator for the given vector
-        /// \param num_components  number of components (i.e. phases) in the vector
-        /// \param pinfo           In a parallel this holds the information about the data distribution.
-        template <class Iterator>
-        inline
-        double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() )
-        {
-            static_cast<void>(num_components); // Suppress warning in the serial case.
-            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
-#if HAVE_MPI
-            if ( pinfo.type() == typeid(ParallelISTLInformation) )
-            {
-                const ParallelISTLInformation& info =
-                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
-                typedef typename Iterator::value_type Scalar;
-                Scalar product = 0.0;
-                int size_per_component = (end - it);
-                size_per_component /= num_components; // two lines to supresse unused warning.
-                assert((end - it) == num_components * size_per_component);
-
-                if( num_components == 1 )
-                {
-                    auto component_container =
-                        boost::make_iterator_range(it, end);
-                    info.computeReduction(component_container,
-                                           Opm::Reduction::makeInnerProductFunctor<double>(),
-                                           product);
-                }
-                else
-                {
-                    auto& maskContainer = info.getOwnerMask();
-                    auto mask = maskContainer.begin();
-                    assert(static_cast<int>(maskContainer.size()) == size_per_component);
-
-                    for(int cell = 0; cell < size_per_component; ++cell, ++mask)
-                    {
-                        Scalar cell_product = (*it) * (*it);
-                        ++it;
-                        for(int component=1; component < num_components;
-                            ++component, ++it)
-                        {
-                            cell_product += (*it) * (*it);
-                        }
-                        product += cell_product * (*mask);
-                    }
-                }
-                return info.communicator().sum(product);
-            }
-            else
-#endif
-            {
-                double product = 0.0 ;
-                for( ; it != end; ++it ) {
-                    product += ( *it * *it );
-                }
-                return product;
-            }
-        }
-
-        /// \brief Compute the L-infinity norm of a vector representing a well equation.
-        /// \param a The container to compute the infinity norm on.
-        /// \param info In a parallel this holds the information about the data distribution.
-        inline
-        double infinityNormWell( const ADB& a, const boost::any& pinfo )
-        {
-            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
-            double result=0;
-            if( a.value().size() > 0 ) {
-                result = a.value().matrix().lpNorm<Eigen::Infinity> ();
-            }
-#if HAVE_MPI
-            if ( pinfo.type() == typeid(ParallelISTLInformation) )
-            {
-                const ParallelISTLInformation& real_info =
-                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
-                result = real_info.communicator().max(result);
-            }
-#endif
-            return result;
-        }
-
-    } // namespace detail
-
-
-
-
-
     template <class Grid, class WellModel, class Implementation>
     void
     BlackoilModelBase<Grid, WellModel, Implementation>::

opm/autodiff/BlackoilModelParameters.cpp

@@ -51,6 +51,7 @@ namespace Opm
         update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
         compute_well_potentials_ = param.getDefault("compute_well_potentials", compute_well_potentials_);
         use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
+        deck_file_name_ = param.template get<std::string>("deck_filename");
     }
 
 

opm/autodiff/BlackoilModelParameters.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED
 #define OPM_BLACKOILMODELPARAMETERS_HEADER_INCLUDED
 
+#include <string>
+
 namespace Opm
 {
 
@@ -59,6 +61,9 @@ namespace Opm
         /// Try to detect oscillation or stagnation.
         bool use_update_stabilization_;
 
+        // The file name of the deck
+        std::string deck_file_name_;
+
         /// Construct from user parameters or defaults.
         explicit BlackoilModelParameters( const parameter::ParameterGroup& param );
 

opm/autodiff/DefaultBlackoilSolutionState.hpp

@@ -0,0 +1,57 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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/>.
+*/
+#ifndef OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED
+#define OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+namespace Opm {
+    /// Struct for containing iteration variables.
+    struct DefaultBlackoilSolutionState
+    {
+        typedef AutoDiffBlock<double> ADB;
+        explicit DefaultBlackoilSolutionState(const int np)
+            : pressure  (    ADB::null())
+            , temperature(   ADB::null())
+            , saturation(np, ADB::null())
+            , rs        (    ADB::null())
+            , rv        (    ADB::null())
+            , qs        (    ADB::null())
+            , bhp       (    ADB::null())
+            , wellVariables (    ADB::null())
+            , canonical_phase_pressures(3, ADB::null())
+        {
+        }
+        ADB              pressure;
+        ADB              temperature;
+        std::vector<ADB> saturation;
+        ADB              rs;
+        ADB              rv;
+        ADB              qs;
+        ADB              bhp;
+        ADB wellVariables;
+        // Below are quantities stored in the state for optimization purposes.
+        std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
+    };
+} // namespace Opm
+
+#endif // OPM_DEFAULTBLACKOILSOLUTIONSTATE_HEADER_INCLUDED

opm/autodiff/FlowMain.hpp

@@ -193,6 +193,7 @@ namespace Opm
         typedef BlackoilPropsAdFromDeck FluidProps;
         typedef FluidProps::MaterialLawManager MaterialLawManager;
         typedef typename Simulator::ReservoirState ReservoirState;
+        typedef typename Simulator::OutputWriter OutputWriter;
 
 
         // ------------   Data members   ------------
@@ -229,7 +230,7 @@ namespace Opm
         boost::any parallel_information_;
         // setupOutputWriter()
         std::unique_ptr<EclipseWriter> eclipse_writer_;
-        std::unique_ptr<BlackoilOutputWriter> output_writer_;
+        std::unique_ptr<OutputWriter> output_writer_;
         // setupLinearSolver
         std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
         // createSimulator()
@@ -770,12 +771,12 @@ namespace Opm
             // create output writer after grid is distributed, otherwise the parallel output
             // won't work correctly since we need to create a mapping from the distributed to
             // the global view
-            output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(),
+            output_writer_.reset(new OutputWriter(grid_init_->grid(),
-                                                          param_,
+                                                  param_,
-                                                          *eclipse_state_,
+                                                  *eclipse_state_,
-                                                          std::move(eclipse_writer_),
+                                                  std::move(eclipse_writer_),
-                                                          Opm::phaseUsageFromDeck(*deck_),
+                                                  Opm::phaseUsageFromDeck(*deck_),
-                                                          fluidprops_->permeability()));
+                                                  fluidprops_->permeability()));
         }
 
 

opm/autodiff/FlowMainEbos.hpp

@@ -0,0 +1,141 @@
+/*
+  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 IRIS AS
+  Copyright 2014 STATOIL ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
+#define OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
+
+#include <opm/autodiff/FlowMain.hpp>
+#include <opm/autodiff/BlackoilModelEbos.hpp>
+
+#include <ewoms/version.hh>
+
+namespace Opm
+{
+    // The FlowMain class is the ebos based black-oil simulator.
+    class FlowMainEbos : public FlowMainBase<FlowMainEbos, Dune::CpGrid, Opm::SimulatorFullyImplicitBlackoilEbos>
+    {
+    protected:
+        typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator;
+        typedef FlowMainBase<FlowMainEbos, Dune::CpGrid, Simulator> Base;
+        friend Base;
+
+        typedef typename TTAG(EclFlowProblem) TypeTag;
+        typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator;
+
+        // Print startup message if on output rank.
+        void printStartupMessage()
+        {
+            if (output_cout_) {
+                const int lineLen = 70;
+                const std::string version = moduleVersionName();
+                const std::string banner = "This is flow_ebos (version "+version+")";
+                const std::string ewomsVersion = "(eWoms version: " + Ewoms::versionString() + ")";
+                const int bannerPreLen = (lineLen - 2 - banner.size())/2;
+                const int bannerPostLen = bannerPreLen + (lineLen - 2 - banner.size())%2;
+                const int eVPreLen = (lineLen - 2 - ewomsVersion.size())/2;
+                const int eVPostLen = eVPreLen + (lineLen - 2 - ewomsVersion.size())%2;
+                std::cout << "**********************************************************************\n";
+                std::cout << "*                                                                    *\n";
+                std::cout << "*" << std::string(bannerPreLen, ' ') << banner << std::string(bannerPostLen, ' ') << "*\n";
+                std::cout << "*" << std::string(eVPreLen, ' ') << ewomsVersion << std::string(eVPostLen, ' ') << "*\n";
+                std::cout << "*                                                                    *\n";
+                std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
+                std::cout << "*            and is part of OPM. For more information see:           *\n";
+                std::cout << "*                       http://opm-project.org                       *\n";
+                std::cout << "*                                                                    *\n";
+                std::cout << "**********************************************************************\n\n";
+            }
+        }
+
+        // Parser the input and creates the Deck and EclipseState objects.
+        // Writes to:
+        //   deck_
+        //   eclipse_state_
+        // May throw if errors are encountered, here configured to be somewhat tolerant.
+        void readDeckInput()
+        {
+            std::string progName("flow_ebos");
+            std::string deckFile("--ecl-deck-file-name=");
+            deckFile += param_.get<std::string>("deck_filename");
+            char* ptr[2];
+            ptr[ 0 ] = const_cast< char * > (progName.c_str());
+            ptr[ 1 ] = const_cast< char * > (deckFile.c_str());
+            EbosSimulator::registerParameters();
+            Ewoms::setupParameters_< TypeTag > ( 2, ptr );
+            ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false));
+            ebosSimulator_->model().applyInitialSolution();
+
+            Base::deck_ = ebosSimulator_->gridManager().deck();
+            Base::eclipse_state_ = ebosSimulator_->gridManager().eclState();
+            IOConfig& ioConfig = Base::eclipse_state_->getIOConfig();
+            ioConfig.setOutputDir(Base::output_dir_);
+
+            // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
+            if (Base::param_.has("output_interval")) {
+                const int output_interval = Base::param_.get<int>("output_interval");
+                eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) );
+
+            }
+
+            // Possible to force initialization only behavior (NOSIM).
+            if (Base::param_.has("nosim")) {
+                const bool nosim = Base::param_.get<bool>("nosim");
+                ioConfig.overrideNOSIM( nosim );
+            }
+        }
+
+        // Setup linear solver.
+        // Writes to:
+        //   fis_solver_
+        void setupLinearSolver()
+        {
+            typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType;
+            Base::fis_solver_.reset( new ISTLSolverType( Base::param_, Base::parallel_information_ ) );
+        }
+
+        /// This is the main function of Flow.
+        // Create simulator instance.
+        // Writes to:
+        //   simulator_
+        void createSimulator()
+        {
+            // Create the simulator instance.
+            Base::simulator_.reset(new Simulator(*ebosSimulator_,
+                                                 Base::param_,
+                                                 *Base::geoprops_,
+                                                 *Base::fluidprops_,
+                                                 Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
+                                                 *Base::fis_solver_,
+                                                 Base::gravity_.data(),
+                                                 Base::deck_->hasKeyword("DISGAS"),
+                                                 Base::deck_->hasKeyword("VAPOIL"),
+                                                 Base::eclipse_state_,
+                                                 *Base::output_writer_,
+                                                 Base::threshold_pressures_));
+        }
+
+    private:
+        std::unique_ptr<EbosSimulator> ebosSimulator_;
+    };
+} // namespace Opm
+
+#endif // OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED

opm/autodiff/ISTLSolver.hpp

@@ -0,0 +1,413 @@
+/*
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_ISTLSOLVER_HEADER_INCLUDED
+#define OPM_ISTLSOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AdditionalObjectDeleter.hpp>
+#include <opm/autodiff/CPRPreconditioner.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
+#include <opm/autodiff/NewtonIterationUtilities.hpp>
+#include <opm/autodiff/ParallelRestrictedAdditiveSchwarz.hpp>
+#include <opm/autodiff/ParallelOverlappingILU0.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+
+#include <opm/common/Exceptions.hpp>
+#include <opm/core/linalg/ParallelIstlInformation.hpp>
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/paamg/amg.hh>
+
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+namespace Dune
+{
+
+namespace ISTLUtility {
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
+{
+    FieldMatrix<K,1,1> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
+{
+    FieldMatrix<K,2,2> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling FMatrixHelp::invert
+template <typename K>
+static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
+{
+    FieldMatrix<K,3,3> A ( matrix );
+    FMatrixHelp::invertMatrix(A, matrix );
+}
+
+//! invert matrix by calling matrix.invert
+template <typename K, int n>
+static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
+{
+    matrix.invert();
+}
+
+} // end ISTLUtility
+
+template <class Scalar, int n, int m>
+class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
+{
+public:
+    typedef Dune::FieldMatrix<Scalar, n, m>  BaseType;
+
+    using BaseType :: operator= ;
+    using BaseType :: rows;
+    using BaseType :: cols;
+    explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
+    void invert()
+    {
+        ISTLUtility::invertMatrix( *this );
+    }
+    const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
+    BaseType& asBase() { return static_cast< BaseType& > (*this); }
+};
+
+template<class K, int n, int m>
+void
+print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
+           typename FieldMatrix<K,n,m>::size_type I,
+           typename FieldMatrix<K,n,m>::size_type J,
+           typename FieldMatrix<K,n,m>::size_type therow, int width,
+           int precision)
+{
+    print_row(s, A.asBase(), I, J, therow, width, precision);
+}
+
+template<class K, int n, int m>
+K& firstmatrixelement (MatrixBlock<K,n,m>& A)
+{
+   return firstmatrixelement( A.asBase() );
+}
+
+
+
+template<typename Scalar, int n, int m>
+struct MatrixDimension< MatrixBlock< Scalar, n, m > >
+: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
+{
+};
+
+} // end namespace Dune
+
+namespace Opm
+{
+    /// This class solves the fully implicit black-oil system by
+    /// solving the reduced system (after eliminating well variables)
+    /// as a block-structured matrix (one block for all cell variables) for a fixed
+    /// number of cell variables np .
+    template < class MatrixBlockType, class VectorBlockType >
+    class ISTLSolver : public NewtonIterationBlackoilInterface
+    {
+        typedef typename MatrixBlockType :: field_type  Scalar;
+
+        typedef Dune::BCRSMatrix <MatrixBlockType>      Matrix;
+        typedef Dune::BlockVector<VectorBlockType>      Vector;
+
+    public:
+        typedef Dune::AssembledLinearOperator< Matrix, Vector, Vector > AssembledLinearOperatorType;
+
+        typedef NewtonIterationBlackoilInterface :: SolutionVector  SolutionVector;
+        /// Construct a system solver.
+        /// \param[in] param   parameters controlling the behaviour of the linear solvers
+        /// \param[in] parallelInformation In the case of a parallel run
+        ///                                with dune-istl the information about the parallelization.
+        ISTLSolver(const NewtonIterationBlackoilInterleavedParameters& param,
+                   const boost::any& parallelInformation_arg=boost::any())
+        : iterations_( 0 ),
+          parallelInformation_(parallelInformation_arg),
+          isIORank_(isIORank(parallelInformation_arg)),
+          parameters_( param )
+        {
+        }
+
+        /// Construct a system solver.
+        /// \param[in] param   ParameterGroup controlling the behaviour of the linear solvers
+        /// \param[in] parallelInformation In the case of a parallel run
+        ///                                with dune-istl the information about the parallelization.
+        ISTLSolver(const parameter::ParameterGroup& param,
+                   const boost::any& parallelInformation_arg=boost::any())
+        : iterations_( 0 ),
+          parallelInformation_(parallelInformation_arg),
+          isIORank_(isIORank(parallelInformation_arg)),
+          parameters_( param )
+        {
+        }
+
+        // dummy method that is not implemented for this class
+        SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
+        {
+            OPM_THROW(std::logic_error,"This method is not implemented");
+            return SolutionVector();
+        }
+
+        /// Solve the system of linear equations Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] residual   residual object containing A and b.
+        /// \return               the solution x
+
+        /// \copydoc NewtonIterationBlackoilInterface::iterations
+        int iterations () const { return iterations_; }
+
+        /// \copydoc NewtonIterationBlackoilInterface::parallelInformation
+        const boost::any& parallelInformation() const { return parallelInformation_; }
+
+    public:
+        /// \brief construct the CPR preconditioner and the solver.
+        /// \tparam P The type of the parallel information.
+        /// \param parallelInformation the information about the parallelization.
+        template<int category=Dune::SolverCategory::sequential, class LinearOperator, class POrComm>
+        void constructPreconditionerAndSolve(LinearOperator& linearOperator,
+                                             Vector& x, Vector& istlb,
+                                             const POrComm& parallelInformation_arg,
+                                             Dune::InverseOperatorResult& result) const
+        {
+            // Construct scalar product.
+            typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
+            typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
+            SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
+
+            // Communicate if parallel.
+            parallelInformation_arg.copyOwnerToAll(istlb, istlb);
+
+#if ! HAVE_UMFPACK
+            const bool useAmg = false ;
+            if( useAmg )
+            {
+                typedef ISTLUtility::CPRSelector< Matrix, Vector, Vector, POrComm>  CPRSelectorType;
+                typedef typename CPRSelectorType::AMG AMG;
+                typedef typename CPRSelectorType::Operator MatrixOperator;
+
+                std::unique_ptr< AMG > amg;
+                std::unique_ptr< MatrixOperator > opA;
+
+                if( ! std::is_same< LinearOperator, MatrixOperator > :: value )
+                {
+                    // create new operator in case linear operator and matrix operator differ
+                    opA.reset( CPRSelectorType::makeOperator( linearOperator.getmat(), parallelInformation_arg ) );
+                }
+
+                const double relax = 1.0;
+
+                // Construct preconditioner.
+                constructAMGPrecond( linearOperator, parallelInformation_arg, amg, opA, relax );
+
+                // Solve.
+                solve(linearOperator, x, istlb, *sp, *amg, result);
+            }
+            else
+#endif
+            {
+                // Construct preconditioner.
+                auto precond = constructPrecond(linearOperator, parallelInformation_arg);
+
+                // Solve.
+                solve(linearOperator, x, istlb, *sp, *precond, result);
+            }
+        }
+
+        typedef Dune::SeqILU0<Matrix, Vector, Vector> SeqPreconditioner;
+
+        template <class Operator>
+        std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
+        {
+            const double relax = 0.9;
+            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
+            return precond;
+        }
+
+#if HAVE_MPI
+        typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
+        typedef ParallelOverlappingILU0<Matrix,Vector,Vector,Comm> ParPreconditioner;
+        template <class Operator>
+        std::unique_ptr<ParPreconditioner>
+        constructPrecond(Operator& opA, const Comm& comm) const
+        {
+            typedef std::unique_ptr<ParPreconditioner> Pointer;
+            const double relax = 0.9;
+            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
+        }
+#endif
+
+        template <class LinearOperator, class MatrixOperator, class POrComm, class AMG >
+        void
+        constructAMGPrecond(LinearOperator& linearOperator, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >& opA, const double relax ) const
+        {
+            ISTLUtility::createAMGPreconditionerPointer( *opA, relax, comm, amg );
+        }
+
+
+        template <class MatrixOperator, class POrComm, class AMG >
+        void
+        constructAMGPrecond(MatrixOperator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg, std::unique_ptr< MatrixOperator >&, const double relax ) const
+        {
+            ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
+        }
+
+        /// \brief Solve the system using the given preconditioner and scalar product.
+        template <class Operator, class ScalarProd, class Precond>
+        void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
+        {
+            // TODO: Revise when linear solvers interface opm-core is done
+            // Construct linear solver.
+            // GMRes solver
+            if ( parameters_.newton_use_gmres_ ) {
+                Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
+                          parameters_.linear_solver_reduction_,
+                          parameters_.linear_solver_restart_,
+                          parameters_.linear_solver_maxiter_,
+                          parameters_.linear_solver_verbosity_);
+                // Solve system.
+                linsolve.apply(x, istlb, result);
+            }
+            else { // BiCGstab solver
+                Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
+                          parameters_.linear_solver_reduction_,
+                          parameters_.linear_solver_maxiter_,
+                          parameters_.linear_solver_verbosity_);
+                // Solve system.
+                linsolve.apply(x, istlb, result);
+            }
+        }
+
+
+        /// Solve the linear system Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] A   matrix A
+        /// \param[inout] x  solution to be computed x
+        /// \param[in] b   right hand side b
+        void solve(Matrix& A, Vector& x, Vector& b ) const
+        {
+            // Parallel version is deactivated until we figure out how to do it properly.
+#if HAVE_MPI
+            if (parallelInformation_.type() == typeid(ParallelISTLInformation))
+            {
+                typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
+                Comm istlComm(info.communicator());
+
+                // Construct operator, scalar product and vectors needed.
+                typedef Dune::OverlappingSchwarzOperator<Matrix, Vector, Vector,Comm> Operator;
+                Operator opA(A, istlComm);
+                solve( opA, x, b, istlComm  );
+            }
+            else
+#endif
+            {
+                // Construct operator, scalar product and vectors needed.
+                Dune::MatrixAdapter< Matrix, Vector, Vector> opA( A );
+                solve( opA, x, b );
+            }
+        }
+
+        /// Solve the linear system Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] A   matrix A
+        /// \param[inout] x  solution to be computed x
+        /// \param[in] b   right hand side b
+        template <class Operator, class Comm >
+        void solve(Operator& opA, Vector& x, Vector& b, Comm& comm) const
+        {
+            Dune::InverseOperatorResult result;
+            // Parallel version is deactivated until we figure out how to do it properly.
+#if HAVE_MPI
+            if (parallelInformation_.type() == typeid(ParallelISTLInformation))
+            {
+                const size_t size = opA.getmat().N();
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
+
+                // As we use a dune-istl with block size np the number of components
+                // per parallel is only one.
+                info.copyValuesTo(comm.indexSet(), comm.remoteIndices(),
+                                  size, 1);
+                // Construct operator, scalar product and vectors needed.
+                constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, b, comm, result);
+            }
+            else
+#endif
+            {
+                OPM_THROW(std::logic_error,"this method if for parallel solve only");
+            }
+
+            checkConvergence( result );
+        }
+
+        /// Solve the linear system Ax = b, with A being the
+        /// combined derivative matrix of the residual and b
+        /// being the residual itself.
+        /// \param[in] A   matrix A
+        /// \param[inout] x  solution to be computed x
+        /// \param[in] b   right hand side b
+        template <class Operator>
+        void solve(Operator& opA, Vector& x, Vector& b ) const
+        {
+            Dune::InverseOperatorResult result;
+            // Construct operator, scalar product and vectors needed.
+            Dune::Amg::SequentialInformation info;
+            constructPreconditionerAndSolve(opA, x, b, info, result);
+            checkConvergence( result );
+        }
+
+        void checkConvergence( const Dune::InverseOperatorResult& result ) const
+        {
+            // store number of iterations
+            iterations_ = result.iterations;
+
+            // Check for failure of linear solver.
+            if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
+                const std::string msg("Convergence failure for linear solver.");
+                if (isIORank_) {
+                    OpmLog::problem(msg);
+                }
+                OPM_THROW_NOLOG(LinearSolverProblem, msg);
+            }
+        }
+    protected:
+        mutable int iterations_;
+        boost::any parallelInformation_;
+        bool isIORank_;
+
+        NewtonIterationBlackoilInterleavedParameters parameters_;
+    }; // end ISTLSolver
+
+} // namespace Opm
+#endif

opm/autodiff/IterationReport.hpp

@@ -0,0 +1,37 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2014, 2015 Statoil ASA.
+  Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2015 NTNU
+
+  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/>.
+*/
+
+#ifndef OPM_ITERATION_REPORT_HEADER_INCLUDED
+#define OPM_ITERATION_REPORT_HEADER_INCLUDED
+
+namespace Opm {
+    /// Class used for reporting the outcome of a nonlinearIteration() call.
+    struct IterationReport
+    {
+        bool failed;
+        bool converged;
+        int linear_iterations;
+        int well_iterations;
+    };
+} // namespace Opm
+
+#endif // OPM_ITERATION_REPORT_HEADER_INCLUDED

opm/autodiff/NewtonIterationBlackoilInterleaved.cpp

@@ -34,13 +34,9 @@
 #include <opm/common/Exceptions.hpp>
 #include <opm/core/linalg/ParallelIstlInformation.hpp>
 
+#include <opm/autodiff/ISTLSolver.hpp>
+
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
-#include <dune/istl/scalarproducts.hh>
-#include <dune/istl/operators.hh>
-#include <dune/istl/preconditioners.hh>
-#include <dune/istl/solvers.hh>
-#include <dune/istl/owneroverlapcopy.hh>
-#include <dune/istl/paamg/amg.hh>
 
 #if HAVE_UMFPACK
 #include <Eigen/UmfPackSupport>
@@ -49,90 +45,6 @@
 #endif
 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
-
-namespace Dune
-{
-
-namespace ISTLUtility {
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,1,1> &matrix)
-{
-    FieldMatrix<K,1,1> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,2,2> &matrix)
-{
-    FieldMatrix<K,2,2> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling FMatrixHelp::invert
-template <typename K>
-static inline void invertMatrix (FieldMatrix<K,3,3> &matrix)
-{
-    FieldMatrix<K,3,3> A ( matrix );
-    FMatrixHelp::invertMatrix(A, matrix );
-}
-
-//! invert matrix by calling matrix.invert
-template <typename K, int n>
-static inline void invertMatrix (FieldMatrix<K,n,n> &matrix)
-{
-    matrix.invert();
-}
-
-} // end ISTLUtility
-
-template <class Scalar, int n, int m>
-class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
-{
-public:
-    typedef Dune::FieldMatrix<Scalar, n, m>  BaseType;
-
-    using BaseType :: operator= ;
-    using BaseType :: rows;
-    using BaseType :: cols;
-    explicit MatrixBlock( const Scalar scalar = 0 ) : BaseType( scalar ) {}
-    void invert()
-    {
-        ISTLUtility::invertMatrix( *this );
-    }
-    const BaseType& asBase() const { return static_cast< const BaseType& > (*this); }
-    BaseType& asBase() { return static_cast< BaseType& > (*this); }
-};
-
-template<class K, int n, int m>
-void
-print_row (std::ostream& s, const MatrixBlock<K,n,m>& A,
-           typename FieldMatrix<K,n,m>::size_type I,
-           typename FieldMatrix<K,n,m>::size_type J,
-           typename FieldMatrix<K,n,m>::size_type therow, int width,
-           int precision)
-{
-    print_row(s, A.asBase(), I, J, therow, width, precision);
-}
-
-template<class K, int n, int m>
-K& firstmatrixelement (MatrixBlock<K,n,m>& A)
-{
-   return firstmatrixelement( A.asBase() );
-}
-
-
-
-template<typename Scalar, int n, int m>
-struct MatrixDimension< MatrixBlock< Scalar, n, m > >
-: public MatrixDimension< typename MatrixBlock< Scalar, n, m >::BaseType >
-{
-};
-
-} // end namespace Dune
-
 namespace Opm
 {
 
@@ -161,9 +73,12 @@ namespace Opm
         typedef Dune::FieldVector<Scalar, np    >       VectorBlockType;
 
         typedef Dune::MatrixBlock<Scalar, np, np >      MatrixBlockType;
+
         typedef Dune::BCRSMatrix <MatrixBlockType>      Mat;
         typedef Dune::BlockVector<VectorBlockType>      Vector;
 
+        typedef Opm::ISTLSolver< MatrixBlockType, VectorBlockType > ISTLSolverType;
+
     public:
         typedef NewtonIterationBlackoilInterface :: SolutionVector  SolutionVector;
         /// Construct a system solver.
@@ -172,9 +87,7 @@ namespace Opm
          ///                               with dune-istl the information about the parallelization.
         NewtonIterationBlackoilInterleavedImpl(const NewtonIterationBlackoilInterleavedParameters& param,
                                                const boost::any& parallelInformation_arg=boost::any())
-        : iterations_( 0 ),
+        : istlSolver_( param, parallelInformation_arg ),
-          parallelInformation_(parallelInformation_arg),
-          isIORank_(isIORank(parallelInformation_arg)),
           parameters_( param )
         {
         }
@@ -186,110 +99,12 @@ namespace Opm
         /// \return               the solution x
 
         /// \copydoc NewtonIterationBlackoilInterface::iterations
-        int iterations () const { return iterations_; }
+        int iterations () const { return istlSolver_.iterations(); }
 
         /// \copydoc NewtonIterationBlackoilInterface::parallelInformation
-        const boost::any& parallelInformation() const { return parallelInformation_; }
+        const boost::any& parallelInformation() const { return istlSolver_.parallelInformation(); }
 
     public:
-        /// \brief construct the CPR preconditioner and the solver.
-        /// \tparam P The type of the parallel information.
-        /// \param parallelInformation the information about the parallelization.
-        template<int category=Dune::SolverCategory::sequential, class O, class POrComm>
-        void constructPreconditionerAndSolve(O& opA,
-                                             Vector& x, Vector& istlb,
-                                             const POrComm& parallelInformation_arg,
-                                             Dune::InverseOperatorResult& result) const
-        {
-            // Construct scalar product.
-            typedef Dune::ScalarProductChooser<Vector, POrComm, category> ScalarProductChooser;
-            typedef std::unique_ptr<typename ScalarProductChooser::ScalarProduct> SPPointer;
-            SPPointer sp(ScalarProductChooser::construct(parallelInformation_arg));
-
-            // Communicate if parallel.
-            parallelInformation_arg.copyOwnerToAll(istlb, istlb);
-
-#if ! HAVE_UMFPACK
-            const bool useAmg = false ;
-            if( useAmg )
-            {
-                typedef ISTLUtility::CPRSelector< Mat, Vector, Vector, POrComm>  CPRSelectorType;
-                typedef typename CPRSelectorType::AMG AMG;
-                std::unique_ptr< AMG > amg;
-                // Construct preconditioner.
-                constructAMGPrecond(opA, parallelInformation_arg, amg);
-
-                // Solve.
-                solve(opA, x, istlb, *sp, *amg, result);
-            }
-            else
-#endif
-            {
-                // Construct preconditioner.
-                auto precond = constructPrecond(opA, parallelInformation_arg);
-
-                // Solve.
-                solve(opA, x, istlb, *sp, *precond, result);
-            }
-        }
-
-        typedef Dune::SeqILU0<Mat, Vector, Vector> SeqPreconditioner;
-
-        template <class Operator>
-        std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
-        {
-            const double relax = 0.9;
-            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
-            return precond;
-        }
-
-#if HAVE_MPI
-        typedef Dune::OwnerOverlapCopyCommunication<int, int> Comm;
-        typedef ParallelOverlappingILU0<Mat,Vector,Vector,Comm> ParPreconditioner;
-        template <class Operator>
-        std::unique_ptr<ParPreconditioner>
-        constructPrecond(Operator& opA, const Comm& comm) const
-        {
-            typedef std::unique_ptr<ParPreconditioner> Pointer;
-            const double relax = 0.9;
-            return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
-        }
-#endif
-
-        template <class Operator, class POrComm, class AMG >
-        void
-        constructAMGPrecond(Operator& opA, const POrComm& comm, std::unique_ptr< AMG >& amg ) const
-        {
-            const double relax = 1.0;
-            ISTLUtility::createAMGPreconditionerPointer( opA, relax, comm, amg );
-        }
-
-        /// \brief Solve the system using the given preconditioner and scalar product.
-        template <class Operator, class ScalarProd, class Precond>
-        void solve(Operator& opA, Vector& x, Vector& istlb, ScalarProd& sp, Precond& precond, Dune::InverseOperatorResult& result) const
-        {
-            // TODO: Revise when linear solvers interface opm-core is done
-            // Construct linear solver.
-            // GMRes solver
-            if ( parameters_.newton_use_gmres_ ) {
-                Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond,
-                          parameters_.linear_solver_reduction_,
-                          parameters_.linear_solver_restart_,
-                          parameters_.linear_solver_maxiter_,
-                          parameters_.linear_solver_verbosity_);
-                // Solve system.
-                linsolve.apply(x, istlb, result);
-            }
-            else { // BiCGstab solver
-                Dune::BiCGSTABSolver<Vector> linsolve(opA, sp, precond,
-                          parameters_.linear_solver_reduction_,
-                          parameters_.linear_solver_maxiter_,
-                          parameters_.linear_solver_verbosity_);
-                // Solve system.
-                linsolve.apply(x, istlb, result);
-            }
-        }
-
         void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
                                    Mat& istlA) const
         {
@@ -455,45 +270,8 @@ namespace Opm
             Vector x(istlA.M());
             x = 0.0;
 
-            Dune::InverseOperatorResult result;
+            // solve linear system using ISTL methods
-            // Parallel version is deactivated until we figure out how to do it properly.
+            istlSolver_.solve( istlA, x, istlb );
-#if HAVE_MPI
-            if (parallelInformation_.type() == typeid(ParallelISTLInformation))
-            {
-                typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
-                const ParallelISTLInformation& info =
-                    boost::any_cast<const ParallelISTLInformation&>( parallelInformation_);
-                Comm istlComm(info.communicator());
-                // As we use a dune-istl with block size np the number of components
-                // per parallel is only one.
-                info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices(),
-                                  size, 1);
-                // Construct operator, scalar product and vectors needed.
-                typedef Dune::OverlappingSchwarzOperator<Mat,Vector,Vector,Comm> Operator;
-                Operator opA(istlA, istlComm);
-                constructPreconditionerAndSolve<Dune::SolverCategory::overlapping>(opA, x, istlb, istlComm, result);
-            }
-            else
-#endif
-            {
-                // Construct operator, scalar product and vectors needed.
-                typedef Dune::MatrixAdapter<Mat,Vector,Vector> Operator;
-                Operator opA(istlA);
-                Dune::Amg::SequentialInformation info;
-                constructPreconditionerAndSolve(opA, x, istlb, info, result);
-            }
-
-            // store number of iterations
-            iterations_ = result.iterations;
-
-            // Check for failure of linear solver.
-            if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
-                const std::string msg("Convergence failure for linear solver.");
-                if (isIORank_) {
-                    OpmLog::problem(msg);
-                }
-                OPM_THROW_NOLOG(LinearSolverProblem, msg);
-            }
 
             // Copy solver output to dx.
             for (int i = 0; i < size; ++i) {
@@ -512,10 +290,7 @@ namespace Opm
         }
 
     protected:
-        mutable int iterations_;
+        ISTLSolverType istlSolver_;
-        boost::any parallelInformation_;
-        bool isIORank_;
-
         NewtonIterationBlackoilInterleavedParameters parameters_;
     }; // end NewtonIterationBlackoilInterleavedImpl
 

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -66,7 +66,7 @@ namespace Opm
         {
             newton_use_gmres_        = false;
             linear_solver_reduction_ = 1e-2;
-            linear_solver_maxiter_   = 75;
+            linear_solver_maxiter_   = 150;
             linear_solver_restart_   = 40;
             linear_solver_verbosity_ = 0;
             require_full_sparsity_pattern_ = false;

opm/autodiff/NonlinearSolver.hpp

@@ -24,6 +24,9 @@
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <memory>
 
 namespace Opm {
@@ -40,6 +43,10 @@ namespace Opm {
         typedef ADB::V V;
         typedef ADB::M M;
 
+        typedef double Scalar;
+        typedef Dune::FieldVector<Scalar, 3    >       VectorBlockType;
+        typedef Dune::BlockVector<VectorBlockType>      BVector;
+
         // Available relaxation scheme types.
         enum RelaxType { DAMPEN, SOR };
 
@@ -131,7 +138,17 @@ namespace Opm {
         /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
         std::vector<V>
         computeFluidInPlace(const ReservoirState& x,
-                            const std::vector<int>& fipnum) const;
+                            const std::vector<int>& fipnum) const
+        {
+            return model_->computeFluidInPlace(x, fipnum);
+        }
+
+        std::vector<std::vector<double>>
+        computeFluidInPlace(const std::vector<int>& fipnum) const
+        {
+            return model_->computeFluidInPlace(fipnum);
+        }
+
 
         /// Reference to physical model.
         const PhysicalModel& model() const;
@@ -146,6 +163,8 @@ namespace Opm {
         /// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
         void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const;
 
+        void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const;
+
         /// The greatest relaxation factor (i.e. smallest factor) allowed.
         double relaxMax() const          { return param_.relax_max_; }
 

opm/autodiff/NonlinearSolver_impl.hpp

@@ -24,6 +24,8 @@
 #define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED
 
 #include <opm/autodiff/NonlinearSolver.hpp>
+#include <opm/common/Exceptions.hpp>
+#include <opm/common/ErrorMacros.hpp>
 
 namespace Opm
 {
@@ -99,15 +101,6 @@ namespace Opm
         return wellIterationsLast_;
     }
 
-    template <class PhysicalModel>
-    std::vector<V> 
-    NonlinearSolver<PhysicalModel>::computeFluidInPlace(const ReservoirState& x,
-                                                        const std::vector<int>& fipnum) const
-    {
-        return model_->computeFluidInPlace(x, fipnum);
-    }
-
-
     template <class PhysicalModel>
     int
     NonlinearSolver<PhysicalModel>::
@@ -289,6 +282,38 @@ namespace Opm
         return;
     }
 
+    template <class PhysicalModel>
+    void
+    NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
+    {
+        // The dxOld is updated with dx.
+        // If omega is equal to 1., no relaxtion will be appiled.
+
+        BVector tempDxOld = dxOld;
+        dxOld = dx;
+
+        switch (relaxType()) {
+            case DAMPEN:
+                if (omega == 1.) {
+                    return;
+                }
+                dx *= omega;
+                return;
+            case SOR:
+                if (omega == 1.) {
+                    return;
+                }
+                dx *= omega;
+                tempDxOld *= (1.-omega);
+                dx += tempDxOld;
+                return;
+            default:
+                OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
+        }
+
+        return;
+    }
+
 
 } // namespace Opm
 

opm/autodiff/SimFIBODetails.hpp

@@ -0,0 +1,151 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014-2016 IRIS AS
+  Copyright 2015 Andreas Lauser
+
+  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/>.
+*/
+#ifndef OPM_SIM_FIBO_DETAILS_HPP
+#define OPM_SIM_FIBO_DETAILS_HPP
+
+#include <utility>
+#include <algorithm>
+#include <locale>
+#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
+#include <opm/core/utility/initHydroCarbonState.hpp>
+#include <opm/core/well_controls.h>
+#include <opm/core/wells/DynamicListEconLimited.hpp>
+
+namespace Opm
+{
+    namespace SimFIBODetails {
+        typedef std::unordered_map<std::string, const Well* > WellMap;
+
+        inline WellMap
+        mapWells(const std::vector< const Well* >& wells)
+        {
+            WellMap wmap;
+
+            for (std::vector< const Well* >::const_iterator
+                     w = wells.begin(), e = wells.end();
+                 w != e; ++w)
+            {
+                wmap.insert(std::make_pair((*w)->name(), *w));
+            }
+
+            return wmap;
+        }
+
+        inline int
+        resv_control(const WellControls* ctrl)
+        {
+            int i, n = well_controls_get_num(ctrl);
+
+            bool match = false;
+            for (i = 0; (! match) && (i < n); ++i) {
+                match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
+            }
+
+            if (! match) { i = 0; }
+
+            return i - 1; // -1 if no match, undo final "++" otherwise
+        }
+
+        inline bool
+        is_resv(const Wells& wells,
+                const int    w)
+        {
+            return (0 <= resv_control(wells.ctrls[w]));
+        }
+
+        inline bool
+        is_resv(const WellMap&     wmap,
+                const std::string& name,
+                const std::size_t  step)
+        {
+            bool match = false;
+
+            WellMap::const_iterator i = wmap.find(name);
+
+            if (i != wmap.end()) {
+                const Well* wp = i->second;
+
+                match = (wp->isProducer(step) &&
+                         wp->getProductionProperties(step)
+                         .hasProductionControl(WellProducer::RESV))
+                    ||  (wp->isInjector(step) &&
+                         wp->getInjectionProperties(step)
+                         .hasInjectionControl(WellInjector::RESV));
+            }
+
+            return match;
+        }
+
+        inline std::vector<int>
+        resvWells(const Wells*      wells,
+                  const std::size_t step,
+                  const WellMap&    wmap)
+        {
+            std::vector<int> resv_wells;
+            if( wells )
+            {
+                for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
+                    if (is_resv(*wells, w) ||
+                        ((wells->name[w] != 0) &&
+                         is_resv(wmap, wells->name[w], step)))
+                    {
+                        resv_wells.push_back(w);
+                    }
+                }
+            }
+
+            return resv_wells;
+        }
+
+        inline void
+        historyRates(const PhaseUsage&               pu,
+                     const WellProductionProperties& p,
+                     std::vector<double>&            rates)
+        {
+            assert (! p.predictionMode);
+            assert (rates.size() ==
+                    std::vector<double>::size_type(pu.num_phases));
+
+            if (pu.phase_used[ BlackoilPhases::Aqua ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Aqua ];
+
+                rates[i] = p.WaterRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Liquid ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Liquid ];
+
+                rates[i] = p.OilRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Vapour ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Vapour ];
+
+                rates[i] = p.GasRate;
+            }
+        }
+    } // namespace SimFIBODetails
+} // namespace Opm
+
+#endif // OPM_SIM_FIBO_DETAILS_HPP

opm/autodiff/SimulatorBase_impl.hpp

@@ -27,6 +27,7 @@
 #include <opm/core/utility/initHydroCarbonState.hpp>
 #include <opm/core/well_controls.h>
 #include <opm/core/wells/DynamicListEconLimited.hpp>
+#include <opm/autodiff/BlackoilModel.hpp>
 
 namespace Opm
 {

opm/autodiff/SimulatorFullyImplicitBlackoil.hpp

@@ -22,6 +22,7 @@
 #define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
 
 #include <opm/autodiff/SimulatorBase.hpp>
+#include <opm/autodiff/BlackoilModel.hpp>
 #include <opm/autodiff/NonlinearSolver.hpp>
 
 namespace Opm {

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -0,0 +1,755 @@
+/*
+  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
+  Copyright 2015 Andreas Lauser
+
+  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/>.
+*/
+
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
+
+//#include <opm/autodiff/SimulatorBase.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp>
+#include <opm/autodiff/IterationReport.hpp>
+#include <opm/autodiff/NonlinearSolver.hpp>
+#include <opm/autodiff/BlackoilModelEbos.hpp>
+#include <opm/autodiff/BlackoilModelParameters.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
+#include <opm/autodiff/StandardWellsDense.hpp>
+#include <opm/autodiff/RateConverter.hpp>
+#include <opm/autodiff/SimFIBODetails.hpp>
+
+#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
+#include <opm/core/utility/initHydroCarbonState.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+
+#include <opm/common/Exceptions.hpp>
+#include <opm/common/ErrorMacros.hpp>
+
+namespace Opm {
+
+class SimulatorFullyImplicitBlackoilEbos;
+//class StandardWellsDense<FluidSystem>;
+
+/// a simulator for the blackoil model
+class SimulatorFullyImplicitBlackoilEbos
+{
+public:
+    typedef typename TTAG(EclFlowProblem) TypeTag;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
+
+    typedef WellStateFullyImplicitBlackoilDense WellState;
+    typedef BlackoilState ReservoirState;
+    typedef BlackoilOutputWriterEbos OutputWriter;
+    typedef BlackoilModelEbos Model;
+    typedef BlackoilModelParameters ModelParameters;
+    typedef NonlinearSolver<Model> Solver;
+    typedef StandardWellsDense<FluidSystem, BlackoilIndices> WellModel;
+
+
+    /// Initialise from parameters and objects to observe.
+    /// \param[in] param       parameters, this class accepts the following:
+    ///     parameter (default)            effect
+    ///     -----------------------------------------------------------
+    ///     output (true)                  write output to files?
+    ///     output_dir ("output")          output directoty
+    ///     output_interval (1)            output every nth step
+    ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+    ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+    ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+    ///     nl_maxiter (30)                max nonlinear iterations in transport
+    ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+    ///     num_transport_substeps (1)     number of transport steps per pressure step
+    ///     use_segregation_split (false)  solve for gravity segregation (if false,
+    ///                                    segregation is ignored).
+    ///
+    /// \param[in] geo           derived geological properties
+    /// \param[in] props         fluid and rock properties
+    /// \param[in] rock_comp_props     if non-null, rock compressibility properties
+    /// \param[in] linsolver     linear solver
+    /// \param[in] gravity       if non-null, gravity vector
+    /// \param[in] has_disgas    true for dissolved gas option
+    /// \param[in] has_vapoil    true for vaporized oil option
+    /// \param[in] eclipse_state the object which represents an internalized ECL deck
+    /// \param[in] output_writer
+    /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
+    SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
+                                       const parameter::ParameterGroup& param,
+                                       DerivedGeology& geo,
+                                       BlackoilPropsAdInterface& props,
+                                       const RockCompressibility* rock_comp_props,
+                                       NewtonIterationBlackoilInterface& linsolver,
+                                       const double* gravity,
+                                       const bool has_disgas,
+                                       const bool has_vapoil,
+                                       std::shared_ptr<EclipseState> eclipse_state,
+                                       BlackoilOutputWriterEbos& output_writer,
+                                       const std::vector<double>& threshold_pressures_by_face)
+        : ebosSimulator_(ebosSimulator),
+          param_(param),
+          model_param_(param),
+          solver_param_(param),
+          props_(props),
+          rock_comp_props_(rock_comp_props),
+          gravity_(gravity),
+          geo_(geo),
+          solver_(linsolver),
+          has_disgas_(has_disgas),
+          has_vapoil_(has_vapoil),
+          terminal_output_(param.getDefault("output_terminal", true)),
+          output_writer_(output_writer),
+          threshold_pressures_by_face_(threshold_pressures_by_face),
+          is_parallel_run_( false )
+    {
+
+        // Misc init.
+        const int num_cells = AutoDiffGrid::numCells(grid());
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+
+        rateConverter_.reset(new RateConverterType(props_, std::vector<int>(AutoDiffGrid::numCells(grid()), 0)));
+
+#if HAVE_MPI
+        if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+        {
+            const ParallelISTLInformation& info =
+                boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
+            // Only rank 0 does print to std::cout
+            terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
+            is_parallel_run_ = ( info.communicator().size() > 1 );
+        }
+#endif
+    }
+
+    /// Run the simulation.
+    /// This will run succesive timesteps until timer.done() is true. It will
+    /// modify the reservoir and well states.
+    /// \param[in,out] timer       governs the requested reporting timesteps
+    /// \param[in,out] state       state of reservoir: pressure, fluxes
+    /// \return                    simulation report, with timing data
+    SimulatorReport run(SimulatorTimer& timer,
+                        ReservoirState& state)
+    {
+        WellState prev_well_state;
+
+        if (output_writer_.isRestart()) {
+            // This is a restart, populate WellState and ReservoirState state objects from restart file
+            output_writer_.initFromRestartFile(props_.phaseUsage(), props_.permeability(), grid(), state, prev_well_state);
+            initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
+        }
+
+        // Create timers and file for writing timing info.
+        Opm::time::StopWatch solver_timer;
+        double stime = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
+        std::ofstream tstep_os(tstep_filename.c_str());
+
+        // adaptive time stepping
+        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
+        if( param_.getDefault("timestep.adaptive", true ) )
+        {
+            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
+        }
+
+        std::string restorefilename = param_.getDefault("restorefile", std::string("") );
+        if( ! restorefilename.empty() )
+        {
+            // -1 means that we'll take the last report step that was written
+            //const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
+
+            //            output_writer_.restore( timer,
+            //                                    state,
+            //                                    prev_well_state,
+            //                                    restorefilename,
+            //                                    desiredRestoreStep );
+        }
+
+        unsigned int totalLinearizations = 0;
+        unsigned int totalNonlinearIterations = 0;
+        unsigned int totalLinearIterations = 0;
+        bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
+        std::vector<double> well_potentials;
+        DynamicListEconLimited dynamic_list_econ_limited;
+
+        bool ooip_computed = false;
+        std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
+        //Get compressed cell fipnum.
+        std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
+        if (fipnum_global.empty()) {
+            std::fill(fipnum.begin(), fipnum.end(), 0);
+        } else {
+            for (size_t c = 0; c < fipnum.size(); ++c) {
+                fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
+            }
+        }
+        std::vector<std::vector<double>> OOIP;
+
+        // Main simulation loop.
+        while (!timer.done()) {
+            // Report timestep.
+            step_timer.start();
+            if ( terminal_output_ )
+            {
+                std::ostringstream ss;
+                timer.report(ss);
+                OpmLog::note(ss.str());
+            }
+
+            // Create wells and well state.
+            WellsManager wells_manager(eclState(),
+                                       timer.currentStepNum(),
+                                       Opm::UgGridHelpers::numCells(grid()),
+                                       Opm::UgGridHelpers::globalCell(grid()),
+                                       Opm::UgGridHelpers::cartDims(grid()),
+                                       Opm::UgGridHelpers::dimensions(grid()),
+                                       Opm::UgGridHelpers::cell2Faces(grid()),
+                                       Opm::UgGridHelpers::beginFaceCentroids(grid()),
+                                       props_.permeability(),
+                                       dynamic_list_econ_limited,
+                                       is_parallel_run_,
+                                       well_potentials );
+
+            const Wells* wells = wells_manager.c_wells();
+            WellState well_state;
+            well_state.init(wells, state, prev_well_state);
+
+            // give the polymer and surfactant simulators the chance to do their stuff
+            handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
+
+            // Compute reservoir volumes for RESV controls.
+            computeRESV(timer.currentStepNum(), wells, state, well_state);
+
+            // Run a multiple steps of the solver depending on the time step control.
+            solver_timer.start();
+
+            const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
+            const WellModel well_model(wells, model_param_, terminal_output_, pv);
+
+            auto solver = createSolver(well_model);
+
+            // write the inital state at the report stage
+            if (timer.initialStep()) {
+
+                // calculate Intensive Quantities
+                const auto& gridManager = ebosSimulator_.gridManager();
+                const auto& gridView = gridManager.gridView();
+                auto elemIt = gridView.template begin<0>();
+                auto elemEndIt = gridView.template end<0>();
+                ElementContext elemCtx(ebosSimulator_);
+                for (; elemIt != elemEndIt; ++ elemIt) {
+                    // this is convenient, but slightly inefficient: one only needs to update
+                    // the primary intensive quantities
+                    elemCtx.updateAll(*elemIt);
+                }
+
+                // No per cell data is written for initial step, but will be
+                // for subsequent steps, when we have started simulating
+                output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
+            }
+
+            // Compute orignal FIP;
+            if (!ooip_computed) {
+                OOIP = solver->computeFluidInPlace(fipnum);
+                FIPUnitConvert(eclState().getUnits(), OOIP);
+                ooip_computed = true;
+            }
+
+            if( terminal_output_ )
+            {
+                std::ostringstream step_msg;
+                boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
+                step_msg.imbue(std::locale(std::locale::classic(), facet));
+                step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
+                         << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
+                         << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
+                         << ", date = " << timer.currentDateTime()
+                         << "\n";
+                OpmLog::info(step_msg.str());
+            }
+
+            solver->model().beginReportStep();
+
+            // If sub stepping is enabled allow the solver to sub cycle
+            // in case the report steps are too large for the solver to converge
+            //
+            // \Note: The report steps are met in any case
+            // \Note: The sub stepping will require a copy of the state variables
+            if( adaptiveTimeStepping ) {
+                adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
+            }
+            else {
+                // solve for complete report step
+                solver->step(timer, state, well_state);
+
+                if( terminal_output_ )
+                {
+                    std::ostringstream iter_msg;
+                    iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
+                    if (solver->wellIterations() != 0) {
+                        iter_msg << " days well iterations = " << solver->wellIterations() << ", ";
+                    }
+                    iter_msg << "non-linear iterations = " << solver->nonlinearIterations()
+                             << ", total linear iterations = " << solver->linearIterations()
+                             << "\n";
+                    OpmLog::info(iter_msg.str());
+                }
+            }
+
+            solver->model().endReportStep();
+
+            // take time that was used to solve system for this reportStep
+            solver_timer.stop();
+
+            // accumulate the number of nonlinear and linear Iterations
+            totalLinearizations += solver->linearizations();
+            totalNonlinearIterations += solver->nonlinearIterations();
+            totalLinearIterations += solver->linearIterations();
+
+            // Report timing.
+            const double st = solver_timer.secsSinceStart();
+
+            // Compute current FIP.
+            std::vector<std::vector<double>> COIP;
+            COIP = solver->computeFluidInPlace(fipnum);
+            FIPUnitConvert(eclState().getUnits(), COIP);
+            std::vector<double> OOIP_totals = FIPTotals(OOIP, state);
+            std::vector<double> COIP_totals = FIPTotals(COIP, state);
+
+            if ( terminal_output_ )
+            {
+                outputFluidInPlace(OOIP_totals, COIP_totals,eclState().getUnits(), 0);
+                for (size_t reg = 0; reg < OOIP.size(); ++reg) {
+                    outputFluidInPlace(OOIP[reg], COIP[reg], eclState().getUnits(), reg+1);
+                }
+            }
+
+            // accumulate total time
+            stime += st;
+
+            if ( terminal_output_ )
+            {
+                std::string msg;
+                msg = "Fully implicit solver took: " + std::to_string(st) + " seconds. Total solver time taken: " + std::to_string(stime) + " seconds.";
+                OpmLog::note(msg);
+            }
+
+            if ( output_writer_.output() ) {
+                SimulatorReport step_report;
+                step_report.pressure_time = st;
+                step_report.total_time =  step_timer.secsSinceStart();
+                step_report.reportParam(tstep_os);
+            }
+
+            // Increment timer, remember well state.
+            ++timer;
+
+            // write simulation state at the report stage
+            output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
+
+            prev_well_state = well_state;
+            // The well potentials are only computed if they are needed
+            // For now thay are only used to determine default guide rates for group controlled wells
+            if ( is_well_potentials_computed ) {
+                computeWellPotentials(wells, well_state, well_potentials);
+            }
+
+            updateListEconLimited(solver, eclState().getSchedule(), timer.currentStepNum(), wells,
+                                  well_state, dynamic_list_econ_limited);
+        }
+
+        // Stop timer and create timing report
+        total_timer.stop();
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        report.total_linearizations = totalLinearizations;
+        report.total_newton_iterations = totalNonlinearIterations;
+        report.total_linear_iterations = totalLinearIterations;
+        return report;
+    }
+
+    const Grid& grid() const
+    { return ebosSimulator_.gridManager().grid(); }
+
+protected:
+    void handleAdditionalWellInflow(SimulatorTimer& timer,
+                                    WellsManager& wells_manager,
+                                    WellState& well_state,
+                                    const Wells* wells)
+    { }
+
+    std::unique_ptr<Solver> createSolver(const WellModel& well_model)
+    {
+        auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
+                                                      model_param_,
+                                                      props_,
+                                                      geo_,
+                                                      rock_comp_props_,
+                                                      well_model,
+                                                      solver_,
+                                                      terminal_output_));
+
+        return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
+    }
+
+    void computeRESV(const std::size_t step,
+                     const Wells* wells,
+                     const BlackoilState& x,
+                     WellState& xw)
+    {
+        typedef SimFIBODetails::WellMap WellMap;
+
+        const auto w_ecl = eclState().getSchedule().getWells(step);
+        const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
+
+        const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
+
+        const std::size_t number_resv_wells        = resv_wells.size();
+        std::size_t       global_number_resv_wells = number_resv_wells;
+#if HAVE_MPI
+        if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
+        {
+            const auto& info =
+                boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
+            global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
+            if ( global_number_resv_wells )
+            {
+                // At least one process has resv wells. Therefore rate converter needs
+                // to calculate averages over regions that might cross process
+                // borders. This needs to be done by all processes and therefore
+                // outside of the next if statement.
+                rateConverter_->defineState(x, boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()));
+            }
+        }
+        else
+#endif
+        {
+            if ( global_number_resv_wells )
+            {
+                rateConverter_->defineState(x);
+            }
+        }
+
+        if (! resv_wells.empty()) {
+            const PhaseUsage&                    pu = props_.phaseUsage();
+            const std::vector<double>::size_type np = props_.numPhases();
+
+            std::vector<double> distr (np);
+            std::vector<double> hrates(np);
+            std::vector<double> prates(np);
+
+            for (std::vector<int>::const_iterator
+                     rp = resv_wells.begin(), e = resv_wells.end();
+                 rp != e; ++rp)
+            {
+                WellControls* ctrl = wells->ctrls[*rp];
+                const bool is_producer = wells->type[*rp] == PRODUCER;
+
+                // RESV control mode, all wells
+                {
+                    const int rctrl = SimFIBODetails::resv_control(ctrl);
+
+                    if (0 <= rctrl) {
+                        const std::vector<double>::size_type off = (*rp) * np;
+
+                        if (is_producer) {
+                            // Convert to positive rates to avoid issues
+                            // in coefficient calculations.
+                            std::transform(xw.wellRates().begin() + (off + 0*np),
+                                           xw.wellRates().begin() + (off + 1*np),
+                                           prates.begin(), std::negate<double>());
+                        } else {
+                            std::copy(xw.wellRates().begin() + (off + 0*np),
+                                      xw.wellRates().begin() + (off + 1*np),
+                                      prates.begin());
+                        }
+
+                        const int fipreg = 0; // Hack.  Ignore FIP regions.
+                        rateConverter_->calcCoeff(prates, fipreg, distr);
+
+                        well_controls_iset_distr(ctrl, rctrl, & distr[0]);
+                    }
+                }
+
+                // RESV control, WCONHIST wells.  A bit of duplicate
+                // work, regrettably.
+                if (is_producer && wells->name[*rp] != 0) {
+                    WellMap::const_iterator i = wmap.find(wells->name[*rp]);
+
+                    if (i != wmap.end()) {
+                        const auto* wp = i->second;
+
+                        const WellProductionProperties& p =
+                            wp->getProductionProperties(step);
+
+                        if (! p.predictionMode) {
+                            // History matching (WCONHIST/RESV)
+                            SimFIBODetails::historyRates(pu, p, hrates);
+
+                            const int fipreg = 0; // Hack.  Ignore FIP regions.
+                            rateConverter_->calcCoeff(hrates, fipreg, distr);
+
+                            // WCONHIST/RESV target is sum of all
+                            // observed phase rates translated to
+                            // reservoir conditions.  Recall sign
+                            // convention: Negative for producers.
+                            const double target =
+                                - std::inner_product(distr.begin(), distr.end(),
+                                                     hrates.begin(), 0.0);
+
+                            well_controls_clear(ctrl);
+                            well_controls_assert_number_of_phases(ctrl, int(np));
+
+                            static const double invalid_alq = -std::numeric_limits<double>::max();
+                            static const int invalid_vfp = -std::numeric_limits<int>::max();
+
+                            const int ok_resv =
+                                well_controls_add_new(RESERVOIR_RATE, target,
+                                                      invalid_alq, invalid_vfp,
+                                                      & distr[0], ctrl);
+
+                            // For WCONHIST the BHP limit is set to 1 atm.
+                            // or a value specified using WELTARG
+                            double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
+                            const int ok_bhp =
+                                well_controls_add_new(BHP, bhp_limit,
+                                                      invalid_alq, invalid_vfp,
+                                                      NULL, ctrl);
+
+                            if (ok_resv != 0 && ok_bhp != 0) {
+                                xw.currentControls()[*rp] = 0;
+                                well_controls_set_current(ctrl, 0);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        if( wells )
+        {
+            for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
+                WellControls* ctrl = wells->ctrls[w];
+                const bool is_producer = wells->type[w] == PRODUCER;
+                if (!is_producer && wells->name[w] != 0) {
+                    WellMap::const_iterator i = wmap.find(wells->name[w]);
+                    if (i != wmap.end()) {
+                        const auto* wp = i->second;
+                        const WellInjectionProperties& injector = wp->getInjectionProperties(step);
+                        if (!injector.predictionMode) {
+                            //History matching WCONINJEH
+                            static const double invalid_alq = -std::numeric_limits<double>::max();
+                            static const int invalid_vfp = -std::numeric_limits<int>::max();
+                            // For WCONINJEH the BHP limit is set to a large number
+                            // or a value specified using WELTARG
+                            double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
+                            const int ok_bhp =
+                                well_controls_add_new(BHP, bhp_limit,
+                                                      invalid_alq, invalid_vfp,
+                                                      NULL, ctrl);
+                            if (!ok_bhp) {
+                                OPM_THROW(std::runtime_error, "Failed to add well control.");
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    void computeWellPotentials(const Wells*                    wells,
+                               const WellState& xw,
+                               std::vector<double>& well_potentials)
+    {
+        const int nw = wells->number_of_wells;
+        const int np = wells->number_of_phases;
+        well_potentials.clear();
+        well_potentials.resize(nw*np,0.0);
+        for (int w = 0; w < nw; ++w) {
+            for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                for (int phase = 0; phase < np; ++phase) {
+                    well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
+                }
+            }
+        }
+    }
+
+
+    void updateListEconLimited(const std::unique_ptr<Solver>& solver,
+                               const Schedule& schedule,
+                               const int current_step,
+                               const Wells* wells,
+                               const WellState& well_state,
+                               DynamicListEconLimited& list_econ_limited) const
+    {
+        solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
+                                                          well_state, list_econ_limited);
+    }
+
+    void FIPUnitConvert(const UnitSystem& units,
+                        std::vector<std::vector<double>>& fip)
+    {
+        if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
+            for (size_t i = 0; i < fip.size(); ++i) {
+                fip[i][0] = unit::convert::to(fip[i][0], unit::stb);
+                fip[i][1] = unit::convert::to(fip[i][1], unit::stb);
+                fip[i][2] = unit::convert::to(fip[i][2], 1000*unit::cubic(unit::feet));
+                fip[i][3] = unit::convert::to(fip[i][3], 1000*unit::cubic(unit::feet));
+                fip[i][4] = unit::convert::to(fip[i][4], unit::stb);
+                fip[i][5] = unit::convert::to(fip[i][5], unit::stb);
+                fip[i][6] = unit::convert::to(fip[i][6], unit::psia);
+            }
+        }
+        if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
+            for (size_t i = 0; i < fip.size(); ++i) {
+                fip[i][6] = unit::convert::to(fip[i][6], unit::barsa);
+            }
+        }
+    }
+
+
+    std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip, const ReservoirState& state)
+    {
+        std::vector<double> totals(7,0.0);
+        for (int i = 0; i < 5; ++i) {
+            for (size_t reg = 0; reg < fip.size(); ++reg) {
+                totals[i] += fip[reg][i];
+            }
+        }
+        const int numCells = Opm::AutoDiffGrid::numCells(grid());
+        const auto& pv = geo_.poreVolume();
+        double pv_hydrocarbon_sum = 0.0;
+        double p_pv_hydrocarbon_sum = 0.0;
+
+        for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+            const auto& intQuants = *ebosSimulator_.model().cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
+            const auto& fs = intQuants.fluidState();
+
+            const double& p = fs.pressure(FluidSystem::oilPhaseIdx).value();
+            const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
+            if ( ! is_parallel_run_ )
+            {
+                totals[5] += pv[cellIdx];
+                pv_hydrocarbon_sum += pv[cellIdx] * hydrocarbon;
+                p_pv_hydrocarbon_sum += p * pv[cellIdx] * hydrocarbon;
+            }
+            else {
+                OPM_THROW(std::logic_error, "FIP not yet implemented for MPI");
+            }
+        }
+        totals[6] = unit::convert::to( (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum), unit::barsa);
+        return totals;
+    }
+
+
+
+    void outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
+    {
+        std::ostringstream ss;
+        if (!reg) {
+            ss << "                                                  ===================================================\n"
+               << "                                                  :                   Field Totals                  :\n";
+        } else {
+            ss << "                                                  ===================================================\n"
+               << "                                                  :        FIPNUM report region  "
+               << std::setw(2) << reg << "                 :\n";
+        }
+        if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
+            ss << "                                                  :      PAV  =" << std::setw(14) << cip[6] << " BARSA                 :\n"
+               << std::fixed << std::setprecision(0)
+               << "                                                  :      PORV =" << std::setw(14) << cip[5] << "   RM3                 :\n";
+            if (!reg) {
+                ss << "                                                  : Pressure is weighted by hydrocarbon pore volume :\n"
+                   << "                                                  : Porv volumes are taken at reference conditions  :\n";
+            }
+            ss << "                         :--------------- Oil    SM3 ---------------:-- Wat    SM3 --:--------------- Gas    SM3 ---------------:\n";
+        }
+        if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
+            ss << "                                                  :      PAV  =" << std::setw(14) << cip[6] << "  PSIA                 :\n"
+               << std::fixed << std::setprecision(0)
+               << "                                                  :      PORV =" << std::setw(14) << cip[5] << "   RB                  :\n";
+            if (!reg) {
+                ss << "                                                  : Pressure is weighted by hydrocarbon pore voulme :\n"
+                   << "                                                  : Pore volumes are taken at reference conditions  :\n";
+            }
+            ss << "                         :--------------- Oil    STB ---------------:-- Wat    STB --:--------------- Gas   MSCF ---------------:\n";
+        }
+        ss << "                         :      Liquid        Vapour        Total   :      Total     :      Free        Dissolved       Total   :" << "\n"
+           << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
+           << ":Currently   in place    :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
+           << std::setw(13) << cip[0] << "   :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
+           << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
+           << ":Originally  in place    :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
+           << std::setw(13) << oip[0] << "   :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
+           << ":========================:==========================================:================:==========================================:\n";
+        OpmLog::note(ss.str());
+    }
+
+
+    const EclipseState& eclState() const
+    { return *ebosSimulator_.gridManager().eclState(); }
+
+    EclipseState& eclState()
+    { return *ebosSimulator_.gridManager().eclState(); }
+
+    // Data.
+    Simulator& ebosSimulator_;
+
+    typedef RateConverter::
+    SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
+                               std::vector<int> > RateConverterType;
+    typedef typename Solver::SolverParameters SolverParameters;
+
+    const parameter::ParameterGroup param_;
+    ModelParameters model_param_;
+    SolverParameters solver_param_;
+
+    // Observed objects.
+    BlackoilPropsAdInterface& props_;
+    const RockCompressibility* rock_comp_props_;
+    const double* gravity_;
+    // Solvers
+    DerivedGeology& geo_;
+    NewtonIterationBlackoilInterface& solver_;
+    // Misc. data
+    std::vector<int> allcells_;
+    const bool has_disgas_;
+    const bool has_vapoil_;
+    bool       terminal_output_;
+    // output_writer
+    OutputWriter& output_writer_;
+    std::unique_ptr<RateConverterType> rateConverter_;
+    // Threshold pressures.
+    std::vector<double> threshold_pressures_by_face_;
+    // Whether this a parallel simulation or not
+    bool is_parallel_run_;
+
+};
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.cpp

@@ -0,0 +1,282 @@
+/*
+  Copyright (c) 2014 SINTEF ICT, Applied Mathematics.
+  Copyright (c) 2015-2016 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 "SimulatorFullyImplicitBlackoilOutputEbos.hpp"
+
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/autodiff/Compat.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/BackupRestore.hpp>
+
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+//For OutputWriterHelper
+#include <map>
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
+
+#ifdef HAVE_OPM_GRID
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <dune/common/version.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+#endif
+namespace Opm
+{
+    namespace detail {
+
+        struct WriterCallEbos : public ThreadHandle :: ObjectInterface
+        {
+            BlackoilOutputWriterEbos& writer_;
+            std::unique_ptr< SimulatorTimerInterface > timer_;
+            const SimulationDataContainer state_;
+            const WellState wellState_;
+            data::Solution simProps_;
+            const bool substep_;
+
+            explicit WriterCallEbos( BlackoilOutputWriterEbos& writer,
+                                 const SimulatorTimerInterface& timer,
+                                 const SimulationDataContainer& state,
+                                 const WellState& wellState,
+                                 const data::Solution& simProps,
+                                 bool substep )
+                : writer_( writer ),
+                  timer_( timer.clone() ),
+                  state_( state ),
+                  wellState_( wellState ),
+                  simProps_( simProps ),
+                  substep_( substep )
+            {
+            }
+
+            // callback to writer's serial writeTimeStep method
+            void run ()
+            {
+                // write data
+                writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ );
+            }
+        };
+    }
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepWithoutCellProperties(
+                  const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  bool substep)
+    {
+        data::Solution noCellProperties;
+        writeTimeStepWithCellProperties(timer, localState, localWellState, noCellProperties, substep);
+    }
+
+
+
+
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepWithCellProperties(
+                  const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  const data::Solution& sol,
+                  bool substep)
+    {
+        bool isIORank = output_ ;
+        if( parallelOutput_ && parallelOutput_->isParallel() )
+        {
+            // If this is not the initial write and no substep, then the well
+            // state used in the computation is actually the one of the last
+            // step. We need that well state for the gathering. Otherwise
+            // It an exception with a message like "global state does not
+            // contain well ..." might be thrown.
+            int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
+                (timer.reportStepNum() - 1) : timer.reportStepNum();
+            // collect all solutions to I/O rank
+            isIORank = parallelOutput_->collectToIORank( localState, localWellState, sol, wellStateStepNumber );
+        }
+
+        const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
+        const WellState& wellState  = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
+
+        // serial output is only done on I/O rank
+        if( isIORank )
+        {
+            if( asyncOutput_ ) {
+                // dispatch the write call to the extra thread
+                asyncOutput_->dispatch( detail::WriterCallEbos( *this, timer, state, wellState, sol, substep ) );
+            }
+            else {
+                // just write the data to disk
+                writeTimeStepSerial( timer, state, wellState, sol, substep );
+            }
+        }
+    }
+
+
+
+    void
+    BlackoilOutputWriterEbos::
+    writeTimeStepSerial(const SimulatorTimerInterface& timer,
+                        const SimulationDataContainer& state,
+                        const WellState& wellState,
+                        const data::Solution& sol,
+                        bool substep)
+    {
+        // ECL output
+        if ( eclWriter_ )
+        {
+            const auto& initConfig = eclipseState_.getInitConfig();
+            if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
+                std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
+            } else {
+                data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...)
+                combined_sol.insert(sol.begin(), sol.end());           // ... insert "extra" data (KR, VISC, ...)
+                eclWriter_->writeTimeStep(timer.reportStepNum(),
+                                          substep,
+                                          timer.simulationTimeElapsed(),
+                                          combined_sol,
+                                          wellState.report(phaseUsage_));
+            }
+        }
+
+        // write backup file
+        if( backupfile_.is_open() )
+        {
+            int reportStep      = timer.reportStepNum();
+            int currentTimeStep = timer.currentStepNum();
+            if( (reportStep == currentTimeStep || // true for SimulatorTimer
+                 currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep
+                 timer.done() ) // true for AdaptiveSimulatorTimer at reportStep
+               && lastBackupReportStep_ != reportStep ) // only backup report step once
+            {
+                // store report step
+                lastBackupReportStep_ = reportStep;
+                // write resport step number
+                backupfile_.write( (const char *) &reportStep, sizeof(int) );
+
+                try {
+                    backupfile_ << state;
+
+                    const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
+                    backupfile_ << boWellState;
+                }
+                catch ( const std::bad_cast& e )
+                {
+                }
+
+                backupfile_ << std::flush;
+            }
+        } // end backup
+    }
+
+    void
+    BlackoilOutputWriterEbos::
+    restore(SimulatorTimerInterface& timer,
+            BlackoilState& state,
+            WellStateFullyImplicitBlackoilDense& wellState,
+            const std::string& filename,
+            const int desiredResportStep )
+    {
+        std::ifstream restorefile( filename.c_str() );
+        if( restorefile )
+        {
+            std::cout << "============================================================================"<<std::endl;
+            std::cout << "Restoring from ";
+            if( desiredResportStep < 0 ) {
+                std::cout << "last";
+            }
+            else {
+                std::cout << desiredResportStep;
+            }
+            std::cout << " report step! filename = " << filename << std::endl << std::endl;
+
+            int reportStep;
+            restorefile.read( (char *) &reportStep, sizeof(int) );
+
+            const int readReportStep = (desiredResportStep < 0) ?
+                std::numeric_limits<int>::max() : desiredResportStep;
+
+            while( reportStep <= readReportStep && ! timer.done() && restorefile )
+            {
+                restorefile >> state;
+                restorefile >> wellState;
+
+                // No per cell data is written for restore steps, but will be
+                // for subsequent steps, when we have started simulating
+                writeTimeStepWithoutCellProperties( timer, state, wellState );
+
+                // some output
+                std::cout << "Restored step " << timer.reportStepNum() << " at day "
+                          <<  unit::convert::to(timer.simulationTimeElapsed(),unit::day) << std::endl;
+
+                if( readReportStep == reportStep ) {
+                    break;
+                }
+
+                // if the stream is not valid anymore we just use the last state read
+                if( ! restorefile ) {
+                    std::cerr << "Reached EOF, using last state read!" << std::endl;
+                    break;
+                }
+
+                // try to read next report step
+                restorefile.read( (char *) &reportStep, sizeof(int) );
+
+                // if read failed, exit loop
+                if( ! restorefile ) {
+                    break;
+                }
+
+                // next step
+                timer.advance();
+
+                if( timer.reportStepNum() != reportStep ) {
+                    break;
+                }
+            }
+        }
+        else
+        {
+            std::cerr << "Warning: Couldn't open restore file '" << filename << "'" << std::endl;
+        }
+    }
+
+
+    bool BlackoilOutputWriterEbos::isRestart() const {
+        const auto& initconfig = eclipseState_.getInitConfig();
+        return initconfig.restartRequested();
+    }
+}

opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp

@@ -0,0 +1,744 @@
+/*
+  Copyright (c) 2014 SINTEF ICT, Applied Mathematics.
+  Copyright (c) 2015 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/>.
+*/
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUTEBOS_HEADER_INCLUDED
+#include <opm/core/grid.h>
+#include <opm/core/simulator/SimulatorTimerInterface.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/output/eclipse/EclipseReader.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/wells/DynamicListEconLimited.hpp>
+
+#include <opm/output/eclipse/EclipseWriter.hpp>
+
+#include <opm/autodiff/Compat.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/ParallelDebugOutput.hpp>
+
+#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
+#include <opm/autodiff/ThreadHandle.hpp>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+
+
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <thread>
+#include <memory>
+
+#include <boost/filesystem.hpp>
+
+#ifdef HAVE_OPM_GRID
+#include <dune/grid/CpGrid.hpp>
+#endif
+namespace Opm
+{
+    class BlackoilState;
+
+    /** \brief Wrapper class for VTK, Matlab, and ECL output. */
+    class BlackoilOutputWriterEbos
+    {
+
+    public:
+        // constructor creating different sub writers
+        template <class Grid>
+        BlackoilOutputWriterEbos(const Grid& grid,
+                                 const parameter::ParameterGroup& param,
+                                 const EclipseState& eclipseState,
+                                 std::unique_ptr<EclipseWriter>&& eclWriter,
+                                 const Opm::PhaseUsage &phaseUsage,
+                                 const double* permeability );
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will extract the
+         *        requested output cell properties specified by the RPTRST keyword
+         *        and write these to file.
+         */
+        template<class Model>
+        void writeTimeStep(const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           const Model& physicalModel,
+                           bool substep = false);
+
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will write all
+         *        CellData in simProps to the file as well.
+         */
+        void writeTimeStepWithCellProperties(
+                           const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           const data::Solution& sol,
+                           bool substep = false);
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection with
+         *        visualization tools like ResInsight. This function will not write
+         *        any cell properties (e.g., those requested by RPTRST keyword)
+         */
+        void writeTimeStepWithoutCellProperties(
+                           const SimulatorTimerInterface& timer,
+                           const SimulationDataContainer& reservoirState,
+                           const Opm::WellState& wellState,
+                           bool substep = false);
+
+        /*!
+         * \brief Write a blackoil reservoir state to disk for later inspection withS
+         *        visualization tools like ResInsight. This is the function which does
+         *        the actual write to file.
+         */
+        void writeTimeStepSerial(const SimulatorTimerInterface& timer,
+                                 const SimulationDataContainer& reservoirState,
+                                 const Opm::WellState& wellState,
+                                 const data::Solution& simProps,
+                                 bool substep);
+
+        /** \brief return output directory */
+        const std::string& outputDirectory() const { return outputDir_; }
+
+        /** \brief return true if output is enabled */
+        bool output () const { return output_; }
+
+        void restore(SimulatorTimerInterface& timer,
+                     BlackoilState& state,
+                     WellStateFullyImplicitBlackoilDense& wellState,
+                     const std::string& filename,
+                     const int desiredReportStep);
+
+
+        template <class Grid>
+        void initFromRestartFile(const PhaseUsage& phaseusage,
+                                 const double* permeability,
+                                 const Grid& grid,
+                                 SimulationDataContainer& simulatorstate,
+                                 WellStateFullyImplicitBlackoilDense& wellstate);
+
+        bool isRestart() const;
+
+    protected:
+        const bool output_;
+        std::unique_ptr< ParallelDebugOutputInterface > parallelOutput_;
+
+        // Parameters for output.
+        const std::string outputDir_;
+        const int output_interval_;
+
+        int lastBackupReportStep_;
+
+        std::ofstream backupfile_;
+        Opm::PhaseUsage phaseUsage_;
+        std::unique_ptr<EclipseWriter> eclWriter_;
+        const EclipseState& eclipseState_;
+
+        std::unique_ptr< ThreadHandle > asyncOutput_;
+    };
+
+
+    //////////////////////////////////////////////////////////////
+    //
+    //  Implementation
+    //
+    //////////////////////////////////////////////////////////////
+    template <class Grid>
+    inline
+    BlackoilOutputWriterEbos::
+    BlackoilOutputWriterEbos(const Grid& grid,
+                             const parameter::ParameterGroup& param,
+                             const Opm::EclipseState& eclipseState,
+                             std::unique_ptr<EclipseWriter>&& eclWriter,
+                             const Opm::PhaseUsage &phaseUsage,
+                             const double* permeability )
+      : output_( param.getDefault("output", true) ),
+        parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability, phaseUsage ) : 0 ),
+        outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ),
+        output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ),
+        lastBackupReportStep_( -1 ),
+        phaseUsage_( phaseUsage ),
+        eclipseState_(eclipseState),
+        asyncOutput_()
+    {
+        // For output.
+        if (output_ && parallelOutput_->isIORank() ) {
+            eclWriter_ = std::move(eclWriter);
+
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(outputDir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+
+            // create output thread if enabled and rank is I/O rank
+            // async output is enabled by default if pthread are enabled
+#if HAVE_PTHREAD
+            const bool asyncOutputDefault = false;
+#else
+            const bool asyncOutputDefault = false;
+#endif
+            if( param.getDefault("async_output", asyncOutputDefault ) )
+            {
+#if HAVE_PTHREAD
+                asyncOutput_.reset( new ThreadHandle() );
+#else
+                OPM_THROW(std::runtime_error,"Pthreads were not found, cannot enable async_output");
+#endif
+            }
+
+            std::string backupfilename = param.getDefault("backupfile", std::string("") );
+            if( ! backupfilename.empty() )
+            {
+                backupfile_.open( backupfilename.c_str() );
+            }
+        }
+    }
+
+
+    template <class Grid>
+    inline void
+    BlackoilOutputWriterEbos::
+    initFromRestartFile( const PhaseUsage& phaseusage,
+                         const double* permeability,
+                         const Grid& grid,
+                         SimulationDataContainer& simulatorstate,
+                         WellStateFullyImplicitBlackoilDense& wellstate)
+    {
+        // gives a dummy dynamic_list_econ_limited
+        DynamicListEconLimited dummy_list_econ_limited;
+        WellsManager wellsmanager(eclipseState_,
+                                  eclipseState_.getInitConfig().getRestartStep(),
+                                  Opm::UgGridHelpers::numCells(grid),
+                                  Opm::UgGridHelpers::globalCell(grid),
+                                  Opm::UgGridHelpers::cartDims(grid),
+                                  Opm::UgGridHelpers::dimensions(grid),
+                                  Opm::UgGridHelpers::cell2Faces(grid),
+                                  Opm::UgGridHelpers::beginFaceCentroids(grid),
+                                  permeability,
+                                  dummy_list_econ_limited
+                                  // We need to pass the optionaly arguments
+                                  // as we get the following error otherwise
+                                  // with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
+                                  // converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructo
+                                  , false,
+                                  std::vector<double>(),
+                                  std::unordered_set<std::string>());
+
+        const Wells* wells = wellsmanager.c_wells();
+        wellstate.resize(wells, simulatorstate); //Resize for restart step
+        auto restarted = Opm::init_from_restart_file(
+                                eclipseState_,
+                                Opm::UgGridHelpers::numCells(grid) );
+
+        solutionToSim( restarted.first, phaseusage, simulatorstate );
+        wellsToState( restarted.second, phaseusage, wellstate );
+    }
+
+
+
+
+
+    namespace detail {
+        template<class Model>
+    void getOutputDataEbos(data::Solution& output,
+                           const Opm::PhaseUsage& phaseUsage,
+                           const Model& model,
+                           const RestartConfig& restartConfig,
+                           const int reportStepNum,
+                           const bool log)
+    {
+            typedef typename Model::FluidSystem FluidSystem;
+
+            //Get the value of each of the keys
+            std::map<std::string, int> outKeywords = restartConfig.getRestartKeywords(reportStepNum);
+            for (auto& keyValue : outKeywords) {
+                keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
+            }
+
+            const auto& ebosModel = model.ebosSimulator().model();
+
+            //Get shorthands for water, oil, gas
+            const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
+            const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
+            const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
+
+            // extract everything which can possibly be written to disk
+            int numCells = ebosModel.numGridDof();
+
+            std::vector<double> pressureOil(numCells);
+            std::vector<double> temperature(numCells);
+
+            std::vector<double> satWater(numCells);
+            std::vector<double> satGas(numCells);
+
+            std::vector<double> bWater(numCells);
+            std::vector<double> bOil(numCells);
+            std::vector<double> bGas(numCells);
+
+            std::vector<double> rhoWater(numCells);
+            std::vector<double> rhoOil(numCells);
+            std::vector<double> rhoGas(numCells);
+
+            std::vector<double> muWater(numCells);
+            std::vector<double> muOil(numCells);
+            std::vector<double> muGas(numCells);
+
+            std::vector<double> krWater(numCells);
+            std::vector<double> krOil(numCells);
+            std::vector<double> krGas(numCells);
+
+            std::vector<double> Rs(numCells);
+            std::vector<double> Rv(numCells);
+            std::vector<double> RsSat(numCells);
+            std::vector<double> RvSat(numCells);
+
+            for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+                const auto& intQuants = *ebosModel.cachedIntensiveQuantities(cellIdx, /*timeIdx=*/0);
+                const auto& fs = intQuants.fluidState();
+
+                pressureOil[cellIdx] = fs.pressure(FluidSystem::oilPhaseIdx).value();
+
+                temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value();
+
+                if (aqua_active) {
+                    satWater[cellIdx] = fs.saturation(FluidSystem::waterPhaseIdx).value();
+                    bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();
+                    rhoWater[cellIdx] = fs.density(FluidSystem::waterPhaseIdx).value();
+                    muWater[cellIdx] = fs.viscosity(FluidSystem::waterPhaseIdx).value();
+                    krWater[cellIdx] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
+                }
+                if (vapour_active) {
+                    satGas[cellIdx] = fs.saturation(FluidSystem::gasPhaseIdx).value();
+                    bGas[cellIdx] = fs.invB(FluidSystem::gasPhaseIdx).value();
+                    rhoGas[cellIdx] = fs.density(FluidSystem::gasPhaseIdx).value();
+                    muGas[cellIdx] = fs.viscosity(FluidSystem::gasPhaseIdx).value();
+                    krGas[cellIdx] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
+                    Rs[cellIdx] = fs.Rs().value();
+                    Rv[cellIdx] = fs.Rv().value();
+                    RsSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
+                                                                             FluidSystem::oilPhaseIdx,
+                                                                             intQuants.pvtRegionIndex(),
+                                                                             /*maxOilSaturation=*/1.0).value();
+                    RvSat[cellIdx] = FluidSystem::saturatedDissolutionFactor(fs,
+                                                                             FluidSystem::gasPhaseIdx,
+                                                                             intQuants.pvtRegionIndex(),
+                                                                             /*maxOilSaturation=*/1.0).value();
+                }
+
+                // oil is always active
+                bOil[cellIdx] = fs.invB(FluidSystem::oilPhaseIdx).value();              
+                rhoOil[cellIdx] = fs.density(FluidSystem::oilPhaseIdx).value();               
+                muOil[cellIdx] = fs.viscosity(FluidSystem::oilPhaseIdx).value();
+                krOil[cellIdx] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
+            }
+
+            /**
+             * Oil Pressures
+             */
+            outKeywords["PRESSURE"] = 0;
+            output.insert("PRESSURE",
+                       UnitSystem::measure::pressure,
+                       std::move(pressureOil),
+                       data::TargetType::RESTART_SOLUTION);
+
+            /**
+             * Temperatures
+             */
+            outKeywords["TEMP"] = 0;
+            output.insert("TEMP",
+                       UnitSystem::measure::temperature,
+                       std::move(temperature),
+                       data::TargetType::RESTART_SOLUTION);
+
+            /**
+             * Water and gas saturation.
+             */
+            outKeywords["SWAT"] = 0;
+            outKeywords["SGAS"] = 0;
+            output.insert("SWAT",
+                       UnitSystem::measure::identity,
+                       std::move(satWater),
+                       data::TargetType::RESTART_SOLUTION);
+            output.insert("SGAS",
+                       UnitSystem::measure::identity,
+                       std::move(satGas),
+                       data::TargetType::RESTART_SOLUTION);
+
+            /**
+             * the dissolution factors
+             */
+            outKeywords["RS"] = 0;
+            outKeywords["RV"] = 0;
+            output.insert("RS",
+                       UnitSystem::measure::gas_oil_ratio,
+                       std::move(Rs),
+                       data::TargetType::RESTART_SOLUTION);
+            output.insert("RV",
+                       UnitSystem::measure::oil_gas_ratio,
+                       std::move(Rv),
+                       data::TargetType::RESTART_SOLUTION);
+
+            /**
+             * Formation volume factors for water, oil, gas
+             */
+            if (aqua_active && outKeywords["BW"] > 0) {
+                outKeywords["BW"] = 0;
+                output.insert("BW",
+                              Opm::UnitSystem::measure::water_inverse_formation_volume_factor,
+                              std::move(bWater),
+                              data::TargetType::RESTART_AUXILLARY);
+
+            }
+            if (liquid_active && outKeywords["BO"]  > 0) {
+                outKeywords["BO"] = 0;
+                output.insert("BO",
+                           Opm::UnitSystem::measure::oil_inverse_formation_volume_factor,
+                           std::move(bOil),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+            if (vapour_active && outKeywords["BG"] > 0) {
+                outKeywords["BG"] = 0;
+                output.insert("BG",
+                           Opm::UnitSystem::measure::gas_inverse_formation_volume_factor,
+                           std::move(bGas),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+
+            /**
+             * Densities for water, oil gas
+             */
+            if (outKeywords["DEN"] > 0) {
+                outKeywords["DEN"] = 0;
+                if (aqua_active) {
+                    output.insert("WAT_DEN",
+                               Opm::UnitSystem::measure::density,
+                               std::move(rhoWater),
+                               data::TargetType::RESTART_AUXILLARY);
+                }
+                if (liquid_active) {
+                    output.insert("OIL_DEN",
+                               Opm::UnitSystem::measure::density,
+                               std::move(rhoOil),
+                               data::TargetType::RESTART_AUXILLARY);
+                }
+                if (vapour_active) {
+                output.insert("GAS_DEN",
+                           Opm::UnitSystem::measure::density,
+                           std::move(rhoGas),
+                           data::TargetType::RESTART_AUXILLARY);
+                }
+            }
+
+            /**
+             * Viscosities for water, oil gas
+             */
+            if (outKeywords["VISC"] > 0) {
+                outKeywords["VISC"] = 0;
+                if (aqua_active) {
+                output.insert("WAT_VISC",
+                           Opm::UnitSystem::measure::viscosity,
+                           std::move(muWater),
+                           data::TargetType::RESTART_AUXILLARY);
+                }
+                if (liquid_active) {
+                output.insert("OIL_VISC",
+                           Opm::UnitSystem::measure::viscosity,
+                           std::move(muOil),
+                           data::TargetType::RESTART_AUXILLARY);
+                }
+                if (vapour_active) {
+                output.insert("GAS_VISC",
+                           Opm::UnitSystem::measure::viscosity,
+                           std::move(muGas),
+                           data::TargetType::RESTART_AUXILLARY);
+                }
+            }
+
+            /**
+             * Relative permeabilities for water, oil, gas
+             */
+            if (aqua_active && outKeywords["KRW"] > 0) {
+                outKeywords["KRW"] = 0;
+                output.insert("WATKR",
+                           Opm::UnitSystem::measure::permeability,
+                           std::move(krWater),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+            if (liquid_active && outKeywords["KRO"] > 0) {
+                outKeywords["KRO"] = 0;
+                output.insert("OILKR",
+                           Opm::UnitSystem::measure::permeability,
+                           std::move(krOil),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+            if (vapour_active && outKeywords["KRG"] > 0) {
+                outKeywords["KRG"] = 0;
+                output.insert("GASKR",
+                           Opm::UnitSystem::measure::permeability,
+                           std::move(krGas),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+
+            /**
+             * Vaporized and dissolved gas/oil ratio
+             */
+            if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
+                outKeywords["RSSAT"] = 0;
+                output.insert("RSSAT",
+                           Opm::UnitSystem::measure::gas_oil_ratio,
+                           std::move(RsSat),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+            if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) {
+                outKeywords["RVSAT"] = 0;
+                output.insert("RVSAT",
+                           Opm::UnitSystem::measure::oil_gas_ratio,
+                           std::move(RvSat),
+                           data::TargetType::RESTART_AUXILLARY);
+            }
+
+
+            /**
+             * Bubble point and dew point pressures
+             */
+            if (log && vapour_active &&
+                    liquid_active && outKeywords["PBPD"] > 0) {
+                Opm::OpmLog::warning("Bubble/dew point pressure output unsupported",
+                                     "Writing bubble points and dew points (PBPD) to file is unsupported, "
+                                     "as the simulator does not use these internally.");
+            }
+
+            //Warn for any unhandled keyword
+            if (log) {
+                for (auto& keyValue : outKeywords) {
+                    if (keyValue.second > 0) {
+                        std::string logstring = "Keyword '";
+                        logstring.append(keyValue.first);
+                        logstring.append("' is unhandled for output to file.");
+                        Opm::OpmLog::warning("Unhandled output keyword", logstring);
+                    }
+                }
+            }
+
+        }
+
+        /**
+         * Checks if the summaryConfig has a keyword with the standardized field, region, or block prefixes.
+         */
+        inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) {
+            std::string field_kw = "F" + keyword;
+            std::string region_kw = "R" + keyword;
+            std::string block_kw = "B" + keyword;
+            return summaryConfig.hasKeyword(field_kw)
+                    || summaryConfig.hasKeyword(region_kw)
+                    || summaryConfig.hasKeyword(block_kw);
+        }
+
+
+        /**
+         * Returns the data as asked for in the summaryConfig
+         */
+        template<class Model>
+        void getSummaryData(data::Solution& output,
+                            const Opm::PhaseUsage& phaseUsage,
+                            const Model& physicalModel,
+                            const SummaryConfig& summaryConfig) {
+
+            const typename Model::FIPData& fip = physicalModel.getFIPData();
+
+            //Get shorthands for water, oil, gas
+            const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
+            const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
+            const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
+
+            /**
+             * Now process all of the summary config files
+             */
+            // Water in place
+            if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) {
+                output.insert("WIP",
+                              Opm::UnitSystem::measure::volume,
+                              fip.fip[Model::FIPData::FIP_AQUA],
+                              data::TargetType::SUMMARY );
+            }
+            if (liquid_active) {
+                const std::vector<double>& oipl = fip.fip[Model::FIPData::FIP_LIQUID];
+                const int size = oipl.size();
+
+                const std::vector<double>& oipg = vapour_active ? fip.fip[Model::FIPData::FIP_VAPORIZED_OIL] : std::vector<double>(size,0.0);
+                std::vector<double> oip = oipl;
+                if (vapour_active) {
+                    oip.insert(oip.end(), oipg.begin(), oipg.end());
+                }
+
+                //Oil in place (liquid phase only)
+                if (hasFRBKeyword(summaryConfig, "OIPL")) {
+                    output.insert("OIPL",
+                                  Opm::UnitSystem::measure::volume,
+                                  oipl,
+                                  data::TargetType::SUMMARY );
+                }
+
+                //Oil in place (gas phase only)
+                if (hasFRBKeyword(summaryConfig, "OIPG")) {
+                    output.insert("OIPG",
+                                  Opm::UnitSystem::measure::volume,
+                                  oipg,
+                                  data::TargetType::SUMMARY );
+                }
+                // Oil in place (in liquid and gas phases)
+                if (hasFRBKeyword(summaryConfig, "OIP")) {
+                    output.insert("OIP",
+                                  Opm::UnitSystem::measure::volume,
+                                  oip,
+                                  data::TargetType::SUMMARY );
+                }
+
+
+
+            }
+            if (vapour_active) {
+                const std::vector<double>& gipg = fip.fip[Model::FIPData::FIP_VAPOUR];
+                const int size = gipg.size();
+
+                const std::vector<double>& gipl= liquid_active ? fip.fip[Model::FIPData::FIP_DISSOLVED_GAS] : std::vector<double>(size,0.0);
+                std::vector<double> gip = gipg;
+                if (liquid_active) {
+                    gip.insert(gip.end(), gipl.begin(), gipl.end());
+                }
+
+                // Gas in place (gas phase only)
+                if (hasFRBKeyword(summaryConfig, "GIPG")) {
+                    output.insert("GIPG",
+                                  Opm::UnitSystem::measure::volume,
+                                  gipg,
+                                  data::TargetType::SUMMARY );
+                }
+                // Gas in place (liquid phase only)
+                if (hasFRBKeyword(summaryConfig, "GIPL")) {
+                    output.insert("GIPL",
+                                  Opm::UnitSystem::measure::volume,
+                                  gipl,
+                                  data::TargetType::SUMMARY );
+                }
+                // Gas in place (in both liquid and gas phases)
+                if (hasFRBKeyword(summaryConfig, "GIP")) {
+                    output.insert("GIP",
+                                  Opm::UnitSystem::measure::volume,
+                                  gip,
+                                  data::TargetType::SUMMARY );
+                }
+            }
+            // Cell pore volume in reservoir conditions
+            if (hasFRBKeyword(summaryConfig, "RPV")) {
+                output.insert("RPV",
+                              Opm::UnitSystem::measure::volume,
+                              fip.fip[Model::FIPData::FIP_PV],
+                              data::TargetType::SUMMARY );
+            }
+            // Pressure averaged value (hydrocarbon pore volume weighted)
+            if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) {
+                output.insert("PRH",
+                              Opm::UnitSystem::measure::pressure,
+                              fip.fip[Model::FIPData::FIP_WEIGHTED_PRESSURE],
+                              data::TargetType::SUMMARY );
+            }
+        }
+
+    }
+
+
+
+
+    template<class Model>
+    inline void
+    BlackoilOutputWriterEbos::
+    writeTimeStep(const SimulatorTimerInterface& timer,
+                  const SimulationDataContainer& localState,
+                  const WellState& localWellState,
+                  const Model& physicalModel,
+                  bool substep)
+    {
+        data::Solution cellData{};
+        const RestartConfig& restartConfig = eclipseState_.getRestartConfig();
+        const SummaryConfig& summaryConfig = eclipseState_.getSummaryConfig();
+        const int reportStepNum = timer.reportStepNum();
+        bool logMessages = output_ && parallelOutput_->isIORank();
+
+        if( output_ && !parallelOutput_->isParallel() )
+        {
+
+            detail::getOutputDataEbos( cellData, phaseUsage_, physicalModel,
+                                    restartConfig, reportStepNum, logMessages );
+            detail::getSummaryData( cellData, phaseUsage_, physicalModel, summaryConfig );
+        }
+        else
+        {
+            if ( logMessages )
+            {
+                std::map<std::string, int> rstKeywords = restartConfig.getRestartKeywords(reportStepNum);
+                std::vector<const char*> keywords =
+                    { "WIP", "OIPL", "OIPG", "OIP", "GIPG", "GIPL", "GIP",
+                      "RPV", "FRPH", "RPRH"};
+
+                std::ostringstream str;
+                str << "Output of restart/summary config not supported in parallel. Requested keywords were ";
+                std::size_t no_kw = 0;
+
+                auto func = [&] (const char* kw)
+                    {
+                        if ( detail::hasFRBKeyword(summaryConfig, kw) )
+                        {
+                            str << kw << " ";
+                            ++ no_kw;
+                        }
+                    };
+
+                std::for_each(keywords.begin(), keywords.end(), func);
+
+                for (auto& keyValue : rstKeywords)
+                {
+                        str << keyValue.first << " ";
+                        ++ no_kw;
+                }
+
+                if ( no_kw )
+                {
+                    Opm::OpmLog::warning("Unhandled ouput request", str.str());
+                }
+            }
+        }
+
+        writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep);
+    }
+}
+#endif

opm/autodiff/StandardWells_impl.hpp

@@ -675,17 +675,18 @@ namespace Opm
                         int table_id = well_controls_iget_vfp(wc, ctrl_index);
 
                         const WellType& well_type = wells().type[w];
+                        const int perf = wells().well_connpos[w]; //first perforation.
                         if (well_type == INJECTOR) {
                             double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                             well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
                         }
                         else if (well_type == PRODUCER) {
                             double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                             well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
                         }
@@ -790,17 +791,18 @@ namespace Opm
                 //Set *BHP* target by calculating bhp from THP
                 const WellType& well_type = wells().type[w];
 
+                const int perf = wells().well_connpos[w]; // first perforation.
                 if (well_type == INJECTOR) {
                     double dp = wellhelpers::computeHydrostaticCorrection(
                             wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
-                            wellPerforationDensities(), gravity_);
+                            wellPerforationDensities()[perf], gravity_);
 
                     xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
                 }
                 else if (well_type == PRODUCER) {
                     double dp = wellhelpers::computeHydrostaticCorrection(
                             wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
-                            wellPerforationDensities(), gravity_);
+                            wellPerforationDensities()[perf], gravity_);
 
                     xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
                 }
@@ -1116,11 +1118,12 @@ namespace Opm
 
                     //Set *BHP* target by calculating bhp from THP
                     const WellType& well_type = wells().type[w];
+                    const int perf = wells().well_connpos[w]; //first perforation
 
                     if (well_type == INJECTOR) {
                         double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
                         const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
                         // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
                         if ( bhp < bhps[w]) {
@@ -1130,7 +1133,7 @@ namespace Opm
                     else if (well_type == PRODUCER) {
                         double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                         const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
                         // apply the strictest of the bhp controlls i.e. largest bhp for producers

opm/autodiff/VFPHelpers.hpp

@@ -25,7 +25,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
-
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 /**
  * This file contains a set of helper functions used by VFPProd / VFPInj.
@@ -35,14 +36,7 @@ namespace detail {
 
 
 typedef AutoDiffBlock<double> ADB;
-
+typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
-
-
-
-
-
-
-
 
 
 /**
@@ -52,7 +46,12 @@ inline double zeroIfNan(const double& value) {
     return (std::isnan(value)) ? 0.0 : value;
 }
 
-
+/**
+ * Returns zero if input value is NaN
+ */
+inline double zeroIfNan(const EvalWell& value) {
+    return (std::isnan(value.value())) ? 0.0 : value.value();
+}
 
 
 

opm/autodiff/VFPInjProperties.cpp

@@ -89,7 +89,36 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
 
 
 
+VFPInjProperties::EvalWell VFPInjProperties::bhp(const int table_id,
+                                                   const EvalWell& aqua,
+                                                   const EvalWell& liquid,
+                                                   const EvalWell& vapour,
+                                                   const double& thp) const {
 
+    //Get the table
+    const VFPInjTable* table = detail::getTable(m_tables, table_id);
+    EvalWell bhp = 0.0;
+
+    //Find interpolation variables
+    EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
+
+    //Compute the BHP for each well independently
+    if (table != nullptr) {
+        //First, find the values to interpolate between
+        //Value of FLO is negative in OPM for producers, but positive in VFP table
+        auto flo_i = detail::findInterpData(flo.value(), table->getFloAxis());
+        auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant
+
+        detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
+
+        bhp = bhp_val.dflo * flo;
+        bhp.setValue(bhp_val.value); // thp is assumed constant i.e.
+    }
+    else {
+        bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table
+    }
+    return bhp;
+}
 
 
 

opm/autodiff/VFPInjProperties.hpp

@@ -23,6 +23,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
 #include <opm/core/wells.h>
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 #include <vector>
 #include <map>
@@ -91,6 +93,13 @@ public:
             const ADB& vapour,
             const ADB& thp) const;
 
+    typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
+    EvalWell bhp(const int table_id,
+            const EvalWell& aqua,
+            const EvalWell& liquid,
+            const EvalWell& vapour,
+            const double& thp) const;
+
     /**
      * Linear interpolation of bhp as a function of the input parameters
      * @param table_id Table number to use

opm/autodiff/VFPProdProperties.cpp

@@ -24,7 +24,8 @@
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
-
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 #include <opm/autodiff/VFPHelpers.hpp>
 
 
@@ -74,7 +75,42 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
 }
 
 
+VFPProdProperties::EvalWell VFPProdProperties::bhp(const int table_id,
+                                                   const EvalWell& aqua,
+                                                   const EvalWell& liquid,
+                                                   const EvalWell& vapour,
+                                                   const double& thp,
+                                                   const double& alq) const {
 
+    //Get the table
+    const VFPProdTable* table = detail::getTable(m_tables, table_id);
+    EvalWell bhp = 0.0;
+
+    //Find interpolation variables
+    EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
+    EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
+    EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
+
+    //Compute the BHP for each well independently
+    if (table != nullptr) {
+        //First, find the values to interpolate between
+        //Value of FLO is negative in OPM for producers, but positive in VFP table
+        auto flo_i = detail::findInterpData(-flo.value(), table->getFloAxis());
+        auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant
+        auto wfr_i = detail::findInterpData( wfr.value(), table->getWFRAxis());
+        auto gfr_i = detail::findInterpData( gfr.value(), table->getGFRAxis());
+        auto alq_i = detail::findInterpData( alq, table->getALQAxis()); //assume constant
+
+        detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
+
+        bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo);
+        bhp.setValue(bhp_val.value);
+    }
+    else {
+        bhp.setValue(-1e100); //Signal that this value has not been calculated properly, due to "missing" table
+    }
+    return bhp;
+}
 
 
 

opm/autodiff/VFPProdProperties.hpp

@@ -23,6 +23,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
 #include <opm/core/wells.h>
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 #include <vector>
 #include <map>
@@ -99,6 +101,15 @@ public:
             const ADB& thp,
             const ADB& alq) const;
 
+
+    typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
+    EvalWell bhp(const int table_id,
+            const EvalWell& aqua,
+            const EvalWell& liquid,
+            const EvalWell& vapour,
+            const double& thp,
+            const double& alq) const;
+
     /**
      * Linear interpolation of bhp as a function of the input parameters
      * @param table_id Table number to use

opm/autodiff/WellHelpers.hpp

@@ -133,7 +133,7 @@ namespace Opm {
          */
         inline
         double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
-                                            const Vector& well_perforation_densities, const double gravity) {
+                                            const double& rho, const double gravity) {
             if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
             {
                 // This is a well with no perforations.
@@ -144,13 +144,12 @@ namespace Opm {
             }
             const double well_ref_depth = wells.depth_ref[w];
             const double dh = vfp_ref_depth - well_ref_depth;
-            const int perf = wells.well_connpos[w];
-            const double rho = well_perforation_densities[perf];
             const double dp = rho*gravity*dh;
 
             return dp;
         }
 
+
         inline
         Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
                 const Vector& well_perforation_densities, const double gravity) {
@@ -159,7 +158,23 @@ namespace Opm {
 
 #pragma omp parallel for schedule(static)
             for (int i=0; i<nw; ++i) {
-                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity);
+                const int perf = wells.well_connpos[i];
+                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
+            }
+
+            return retval;
+        }
+
+        inline
+        std::vector<double> computeHydrostaticCorrection(const Wells& wells, const std::vector<double>& vfp_ref_depth,
+                                                         const std::vector<double>& well_perforation_densities, const double gravity) {
+            const int nw = wells.number_of_wells;
+            std::vector<double> retval(nw,0.0);
+
+#pragma omp parallel for schedule(static)
+            for (int i=0; i<nw; ++i) {
+                const int perf = wells.well_connpos[i];
+                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
             }
 
             return retval;

opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp

@@ -0,0 +1,183 @@
+/*
+  Copyright 2016 IRIS
+
+  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/>.
+*/
+
+#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED
+#define OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED
+
+
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/autodiff/BlackoilModelEnums.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <vector>
+#include <cassert>
+#include <string>
+#include <utility>
+#include <map>
+#include <algorithm>
+#include <array>
+
+namespace Opm
+{
+
+    /// The state of a set of wells, tailored for use by the fully
+    /// implicit blackoil simulator.
+    class WellStateFullyImplicitBlackoilDense
+        : public WellStateFullyImplicitBlackoil
+    {
+        typedef WellStateFullyImplicitBlackoil  BaseType;
+    public:
+        typedef BaseType :: WellMapType WellMapType;
+
+        using BaseType :: wellRates;
+        using BaseType :: bhp;
+        using BaseType :: perfPress;
+        using BaseType :: wellMap;
+        using BaseType :: numWells;
+        using BaseType :: numPhases;
+        using BaseType :: perfPhaseRates;
+        using BaseType :: currentControls;
+        using BaseType :: wellPotentials;
+
+        /// Allocate and initialize if wells is non-null.  Also tries
+        /// to give useful initial values to the bhp(), wellRates()
+        /// and perfPhaseRates() fields, depending on controls
+        template <class State, class PrevState>
+        void init(const Wells* wells, const State& state, const PrevState& prevState)
+        {
+            // call init on base class
+            BaseType :: init(wells, state, prevState);
+
+            // if there are no well, do nothing in init
+            if (wells == 0) {
+                return;
+            }
+
+            const int nw = wells->number_of_wells;
+            if( nw == 0 ) return ;
+
+            const int np = wells->number_of_phases;
+
+            well_solutions_.clear();
+            well_solutions_.resize(nw * np, 0.0);
+            std::vector<double> g = {1.0,1.0,0.01};
+            for (int w = 0; w < nw; ++w) {
+                WellControls* wc = wells->ctrls[w];
+
+                // The current control in the well state overrides
+                // the current control set in the Wells struct, which
+                // is instead treated as a default.
+                const int current = currentControls()[w];
+                well_controls_set_current( wc, current);
+                const WellType& well_type = wells->type[w];
+
+                switch (well_controls_iget_type(wc, current)) {
+                case THP: // Intentional fall-through
+                case BHP:
+                {
+                    if (well_type == INJECTOR) {
+                        for (int p = 0; p < np; ++p)  {
+                            well_solutions_[w] += wellRates()[np*w + p] * wells->comp_frac[np*w + p];
+                        }
+                    } else {
+                        for (int p = 0; p < np; ++p) {
+                            well_solutions_[w] += g[p] * wellRates()[np*w + p];
+                        }
+                    }
+                }
+                    break;
+
+
+                case RESERVOIR_RATE: // Intentional fall-through
+                case SURFACE_RATE:
+                {
+                    wellSolutions()[w] = bhp()[w];
+
+                }
+                    break;
+                }
+                assert(np == 3);
+                double total_rates = 0.0;
+                for (int p = 0; p < np; ++p)  {
+                    total_rates += g[p] * wellRates()[np*w + p];
+                }
+
+                if(std::abs(total_rates) > 0) {
+                    wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water;
+                    wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas
+                } else {
+                    wellSolutions()[nw + w] = wells->comp_frac[np*w + Water];
+                    wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas];
+                }
+            }
+
+
+            // intialize wells that have been there before
+            // order may change so the mapping is based on the well name
+            if( ! prevState.wellMap().empty() )
+            {
+                typedef typename WellMapType :: const_iterator const_iterator;
+                const_iterator end = prevState.wellMap().end();
+                int nw_old = prevState.bhp().size();
+                for (int w = 0; w < nw; ++w) {
+                    std::string name( wells->name[ w ] );
+                    const_iterator it = prevState.wellMap().find( name );
+                    if( it != end )
+                    {
+                        const int oldIndex = (*it).second[ 0 ];
+                        const int newIndex = w;
+
+                        // wellSolutions
+                        for( int i = 0;  i < np; ++i)
+                        {
+                            wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ];
+                        }
+
+                    }
+                }
+            }
+        }
+
+        template <class State>
+        void resize(const Wells* wells, const State& state) {
+            const WellStateFullyImplicitBlackoilDense dummy_state{}; // Init with an empty previous state only resizes
+            init(wells, state, dummy_state) ;
+        }
+
+
+        /// One rate per phase and well connection.
+        std::vector<double>& wellSolutions() { return well_solutions_; }
+        const std::vector<double>& wellSolutions() const { return well_solutions_; }
+
+        data::Wells report(const PhaseUsage& pu) const override {
+            data::Wells res = BaseType::report(pu);
+            return res;
+        }
+
+    private:
+        std::vector<double> well_solutions_;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILDENSE_HEADER_INCLUDED