Merge pull request #262 from totto82/useOPM-output

Add ecl output from opm-output

ebos/collecttoiorank.hh

@@ -23,6 +23,9 @@
 #ifndef EWOMS_PARALLELSERIALOUTPUT_HH
 #define EWOMS_PARALLELSERIALOUTPUT_HH
 
+#include <opm/output/data/Cells.hpp>
+#include <opm/output/data/Solution.hpp>
+
 //#if HAVE_OPM_GRID
 #include <dune/grid/common/p2pcommunicator.hh>
 #include <dune/grid/utility/persistentcontainer.hh>
@@ -90,9 +93,6 @@ namespace Ewoms
             IndexMapType& localIndexMap_;
             IndexMapStorageType& indexMaps_;
             std::map< const int, const int > globalPosition_;
-#ifndef NDEBUG
-            std::set< int > checkPosition_;
-#endif
 
         public:
             DistributeIndexMapping( const std::vector<int>& globalIndex,
@@ -111,7 +111,7 @@ namespace Ewoms
                     globalPosition_.insert( std::make_pair( globalIndex[ index ], index ) );
                 }
 
-                // on I/O rank we need to create a mapping from local to global
+                // we need to create a mapping from local to global
                 if( ! indexMaps_.empty() )
                 {
                     // for the ioRank create a localIndex to index in global state map
@@ -122,10 +122,6 @@ namespace Ewoms
                     {
                         const int id = distributedGlobalIndex_[ localIndexMap_[ i ] ];
                         indexMap[ i ] = globalPosition_[ id ] ;
-#ifndef NDEBUG
-                        assert( checkPosition_.find( id ) == checkPosition_.end() );
-                        checkPosition_.insert( id );
-#endif
                     }
                 }
             }
@@ -164,10 +160,6 @@ namespace Ewoms
                     buffer.read( globalId );
                     assert( globalPosition_.find( globalId ) != globalPosition_.end() );
                     indexMap[ index ] = globalPosition_[ globalId ];
-#ifndef NDEBUG
-                    assert( checkPosition_.find( globalId ) == checkPosition_.end() );
-                    checkPosition_.insert( globalId );
-#endif
                 }
             }
         };
@@ -178,12 +170,10 @@ namespace Ewoms
             typename GridManager::Grid, typename GridManager::EquilGrid > :: value ;
 
         CollectDataToIORank( const GridManager& gridManager )
-            : toIORankComm_( ),
+            : toIORankComm_( )
-              isIORank_( gridManager.grid().comm().rank() == ioRank ),
-              isParallel_( gridManager.grid().comm().size() > 1 )
         {
             // index maps only have to be build when reordering is needed
-            if( ! needsReordering && ! isParallel_ )
+            if( ! needsReordering && ! isParallel() )
             {
                 return ;
             }
@@ -192,36 +182,35 @@ namespace Ewoms
 
             {
                 std::set< int > send, recv;
+		typedef typename GridManager::EquilGrid::LeafGridView EquilGridView;
+                const EquilGridView equilGridView = gridManager.equilGrid().leafGridView() ;
+
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
+                typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView> EquilElementMapper;
+                EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
+#else
+                typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView, Dune::MCMGElementLayout> EquilElementMapper;
+                EquilElementMapper equilElemMapper(equilGridView);
+#endif
+
+                // We need a mapping from local to global grid, here we
+                // use equilGrid which represents a view on the global grid
+                const size_t globalSize = gridManager.equilGrid().leafGridView().size( 0 );
+                // reserve memory
+                globalCartesianIndex_.resize(globalSize, -1);
+
+                // loop over all elements (global grid) and store Cartesian index
+                auto elemIt = gridManager.equilGrid().leafGridView().template begin<0>();
+                const auto& elemEndIt = gridManager.equilGrid().leafGridView().template end<0>();
+                for (; elemIt != elemEndIt; ++elemIt) {
+                    int elemIdx = equilElemMapper.index(*elemIt );
+                    int cartElemIdx = gridManager.equilCartesianIndexMapper().cartesianIndex(elemIdx);
+                    globalCartesianIndex_[elemIdx] = cartElemIdx;
+                }
+
                 // the I/O rank receives from all other ranks
                 if( isIORank() )
                 {
-                    typedef typename GridManager::EquilGrid::LeafGridView EquilGridView;
-                    const EquilGridView equilGridView = gridManager.equilGrid().leafGridView() ;
-
-#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
-                    typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView> EquilElementMapper;
-                    EquilElementMapper equilElemMapper(equilGridView, Dune::mcmgElementLayout());
-#else
-                    typedef Dune::MultipleCodimMultipleGeomTypeMapper<EquilGridView, Dune::MCMGElementLayout> EquilElementMapper;
-                    EquilElementMapper equilElemMapper(equilGridView);
-#endif
-
-
-                    // the I/O rank needs a picture of the global grid, here we
-                    // use equilGrid which represents a view on the global grid
-                    const size_t globalSize = gridManager.equilGrid().leafGridView().size( 0 );
-                    // reserve memory
-                    globalCartesianIndex_.resize(globalSize, -1);
-
-                    // loop over all elements (global grid) and store Cartesian index
-                    auto elemIt = gridManager.equilGrid().leafGridView().template begin<0>();
-                    const auto& elemEndIt = gridManager.equilGrid().leafGridView().template end<0>();
-                    for (; elemIt != elemEndIt; ++elemIt) {
-                        int elemIdx = equilElemMapper.index(*elemIt );
-                        int cartElemIdx = gridManager.equilCartesianIndexMapper().cartesianIndex(elemIdx);
-                        globalCartesianIndex_[elemIdx] = cartElemIdx;
-                    }
-
                     for(int i=0; i<comm.size(); ++i)
                     {
                         if( i != ioRank )
@@ -254,15 +243,16 @@ namespace Ewoms
                 ElementMapper elemMapper(localGridView);
 #endif
 
-                for( auto it = localGridView.template begin< 0, Dune::Interior_Partition >(),
+                // A mapping for the whole grid (including the ghosts) is needed for restarts
-                     end = localGridView.template end< 0, Dune::Interior_Partition >(); it != end; ++it )
+                for( auto it = localGridView.template begin< 0 >(),
+                     end = localGridView.template end< 0 >(); it != end; ++it )
                 {
                     const auto element = *it ;
                     int elemIdx = elemMapper.index( element );
                     distributedCartesianIndex[elemIdx] = gridManager.cartesianIndex( elemIdx );
 
                     // only store interior element for collection
-                    assert( element.partitionType() == Dune :: InteriorEntity );
+                    //assert( element.partitionType() == Dune :: InteriorEntity );
 
                     localIndexMap_.push_back( elemIdx );
                 }
@@ -270,12 +260,9 @@ namespace Ewoms
                 // insert send and recv linkage to communicator
                 toIORankComm_.insertRequest( send, recv );
 
-                if( isIORank() )
+                // need an index map for each rank
-                {
+                indexMaps_.clear();
-                    // need an index map for each rank
+                indexMaps_.resize( comm.size() );
-                    indexMaps_.clear();
-                    indexMaps_.resize( comm.size() );
-                }
 
                 // distribute global id's to io rank for later association of dof's
                 DistributeIndexMapping distIndexMapping( globalCartesianIndex_, distributedCartesianIndex, localIndexMap_, indexMaps_ );
@@ -283,33 +270,40 @@ namespace Ewoms
             }
         }
 
-        template <class BufferList>
+        class PackUnPack : public P2PCommunicatorType::DataHandleInterface
-        class PackUnPackOutputBuffers : public P2PCommunicatorType::DataHandleInterface
         {
-            BufferList& bufferList_;
+            const Opm::data::Solution& localCellData_;
+            Opm::data::Solution& globalCellData_;
 
             const IndexMapType& localIndexMap_;
             const IndexMapStorageType& indexMaps_;
 
         public:
-            PackUnPackOutputBuffers( BufferList& bufferList,
+            PackUnPack( const Opm::data::Solution& localCellData,
-                                     const IndexMapType& localIndexMap,
+                        Opm::data::Solution& globalCellData,
-                                     const IndexMapStorageType& indexMaps,
+                        const IndexMapType& localIndexMap,
-                                     const size_t globalSize,
+                        const IndexMapStorageType& indexMaps,
-                                     const bool isIORank )
+                        const size_t globalSize,
-            : bufferList_( bufferList ),
+                        const bool isIORank )
+            : localCellData_( localCellData ),
+              globalCellData_( globalCellData ),
               localIndexMap_( localIndexMap ),
               indexMaps_( indexMaps )
             {
                 if( isIORank )
                 {
+                    // add missing data to global cell data
+                    for (const auto& pair : localCellData_) {
+                        const std::string& key = pair.first;
+                        std::size_t container_size = globalSize;
+                        auto OPM_OPTIM_UNUSED ret = globalCellData_.insert(key, pair.second.dim,
+                                                                           std::vector<double>(container_size),
+                                                                           pair.second.target);
+                        assert(ret.second);
+                    }
+
                     MessageBufferType buffer;
                     pack( 0, buffer );
-                    // resize all buffers
-                    for (auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
-                    {
-                      it->second->resize( globalSize );
-                    }
 
                     // the last index map is the local one
                     doUnpack( indexMaps.back(), buffer );
@@ -324,22 +318,26 @@ namespace Ewoms
                     OPM_THROW(std::logic_error,"link in method pack is not 0 as expected");
                 }
 
-                size_t buffers = bufferList_.size();
+                // write all cell data registered in local state
-                buffer.write( buffers );
+                for (const auto& pair : localCellData_) {
-                for (auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
+                    const auto& data = pair.second.data;
-                {
+
-                    write( buffer, localIndexMap_, *(it->second) );
+                    // write all data from local data to buffer
+                    write( buffer, localIndexMap_, data);
                 }
+
             }
 
             void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
             {
-                size_t buffers = 0;
+                // we loop over the data  as
-                buffer.read( buffers );
+                // its order governs the order the data got received.
-                assert( buffers == bufferList_.size() );
+                for (auto& pair : localCellData_) {
-                for( auto it = bufferList_.begin(), end = bufferList_.end(); it != end; ++it )
+                    const std::string& key = pair.first;
-                {
+                    auto& data = globalCellData_.data(key);
-                    read( buffer, indexMap, *(it->second) );
+
+                    //write all data from local cell data to buffer
+                    read( buffer, indexMap, data);
                 }
             }
 
@@ -351,71 +349,63 @@ namespace Ewoms
 
         protected:
             template <class Vector>
-            void write( MessageBufferType& buffer, const IndexMapType& localIndexMap, const Vector& data ) const
+            void write( MessageBufferType& buffer,
+                        const IndexMapType& localIndexMap,
+                        const Vector& vector,
+                        const unsigned int offset = 0,
+                        const unsigned int stride = 1 ) const
             {
-                const size_t size = localIndexMap.size();
+                unsigned int size = localIndexMap.size();
-                assert( size <= data.size() );
                 buffer.write( size );
-                for( size_t i=0; i<size; ++i )
+                assert( vector.size() >= stride * size );
+                for( unsigned int i=0; i<size; ++i )
                 {
-                    buffer.write( data[ localIndexMap[ i ] ] );
+                    const unsigned int index = localIndexMap[ i ] * stride + offset;
+                    assert( index < vector.size() );
+                    buffer.write( vector[ index ] );
                 }
             }
 
             template <class Vector>
-            void read( MessageBufferType& buffer, const IndexMapType& indexMap, Vector& data ) const
+            void read( MessageBufferType& buffer,
+                       const IndexMapType& indexMap,
+                       Vector& vector,
+                       const unsigned int offset = 0, const unsigned int stride = 1 ) const
             {
-                size_t size = indexMap.size();
+                unsigned int size = 0;
-                assert( size <= data.size() );
                 buffer.read( size );
                 assert( size == indexMap.size() );
-                for( size_t i=0; i<size; ++i )
+                for( unsigned int i=0; i<size; ++i )
                 {
-                    buffer.read( data[ indexMap[ i ] ] );
+                    const unsigned int index = indexMap[ i ] * stride + offset;
+                    assert( index < vector.size() );
+                    buffer.read( vector[ index ] );
                 }
             }
 
-            void writeString( MessageBufferType& buffer, const std::string& s) const
-            {
-                const int size = s.size();
-                buffer.write( size );
-                for( int i=0; i<size; ++i )
-                {
-                    buffer.write( s[ i ] );
-                }
-            }
 
-            void readString( MessageBufferType& buffer, std::string& s) const
-            {
-                int size = -1;
-                buffer.read( size );
-                s.resize( size );
-                for( int i=0; i<size; ++i )
-                {
-                    buffer.read( s[ i ] );
-                }
-            }
         };
 
         // gather solution to rank 0 for EclipseWriter
-        template <class BufferList>
+        void collect( const Opm::data::Solution& localCellData )
-        void collect( BufferList& bufferList ) const
         {
+            globalCellData_ = {};
             // index maps only have to be build when reordering is needed
-            if( ! needsReordering && ! isParallel_ )
+            if( ! needsReordering && ! isParallel() )
             {
                 return ;
             }
 
             // this also packs and unpacks the local buffers one ioRank
-            PackUnPackOutputBuffers< BufferList >
+            PackUnPack
-                packUnpack( bufferList,
+                packUnpack( localCellData,
+                            globalCellData_,
                             localIndexMap_,
                             indexMaps_,
                             numCells(),
                             isIORank() );
 
-            if ( ! isParallel_ )
+            if ( ! isParallel() )
             {
                 // no need to collect anything.
                 return;
@@ -430,9 +420,14 @@ namespace Ewoms
 #endif
         }
 
+        const Opm::data::Solution& globalCellData() const
+        {
+            return globalCellData_;
+        }
+
         bool isIORank() const
         {
-            return isIORank_;
+            return toIORankComm_.rank() == ioRank;
         }
 
         bool isParallel() const
@@ -440,6 +435,23 @@ namespace Ewoms
             return toIORankComm_.size() > 1;
         }
 
+        int localIdxToGlobalIdx(const unsigned localIdx) {
+
+            if ( ! isParallel() )
+            {
+                return localIdx;
+            }
+            // the last indexMap is the local one
+            IndexMapType& indexMap = indexMaps_.back();
+            if( indexMap.empty() )
+                OPM_THROW(std::logic_error,"index map is not created on this rank");
+
+            if (localIdx > indexMap.size())
+                OPM_THROW(std::logic_error,"local index is outside map range");
+
+            return indexMap[localIdx];
+        }
+
         size_t numCells () const { return globalCartesianIndex_.size(); }
 
     protected:
@@ -447,10 +459,7 @@ namespace Ewoms
         IndexMapType                    globalCartesianIndex_;
         IndexMapType                    localIndexMap_;
         IndexMapStorageType             indexMaps_;
-        // true if we are on I/O rank
+        Opm::data::Solution             globalCellData_;
-        bool                            isIORank_;
-        /// \brief True if there is more than one MPI process
-        bool                            isParallel_;
     };
 
 } // end namespace Opm

ebos/ecloutputblackoilmodule.hh

@@ -28,13 +28,15 @@
 #define EWOMS_ECL_OUTPUT_BLACK_OIL_MODULE_HH
 
 #include "eclwriter.hh"
-#include "ecldeckunits.hh"
 
-#include <ewoms/io/baseoutputmodule.hh>
 #include <ewoms/common/propertysystem.hh>
 #include <ewoms/common/parametersystem.hh>
 
 #include <opm/common/Valgrind.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
+#include <opm/output/data/Cells.hpp>
+#include <opm/output/eclipse/EclipseIO.hpp>
+
 
 #include <dune/common/fvector.hh>
 
@@ -42,26 +44,9 @@
 
 namespace Ewoms {
 namespace Properties {
-// create new type tag for the VTK multi-phase output
+// create new type tag for the Ecl-output
 NEW_TYPE_TAG(EclOutputBlackOil);
 
-// create the property tags needed for the multi phase module
-NEW_PROP_TAG(EclOutputWriteSaturations);
-NEW_PROP_TAG(EclOutputWritePressures);
-NEW_PROP_TAG(EclOutputWriteGasDissolutionFactor);
-NEW_PROP_TAG(EclOutputWriteOilVaporizationFactor);
-NEW_PROP_TAG(EclOutputWriteGasFormationVolumeFactor);
-NEW_PROP_TAG(EclOutputWriteOilFormationVolumeFactor);
-NEW_PROP_TAG(EclOutputWriteOilSaturationPressure);
-
-// set default values for what quantities to output
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteSaturations, true);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWritePressures, true);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasDissolutionFactor, true);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilVaporizationFactor, false);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteGasFormationVolumeFactor, true);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilFormationVolumeFactor, true);
-SET_BOOL_PROP(EclOutputBlackOil, EclOutputWriteOilSaturationPressure, false);
 } // namespace Properties
 
 // forward declaration
@@ -75,20 +60,18 @@ class EcfvDiscretization;
  *        ECL binary format.
  */
 template <class TypeTag>
-class EclOutputBlackOilModule : public BaseOutputModule<TypeTag>
+class EclOutputBlackOilModule
 {
-    typedef BaseOutputModule<TypeTag> ParentType;
 
     typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
     typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
     typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
-
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
 
-    typedef Ewoms::EclWriter<TypeTag> EclWriter;
-
     enum { numPhases = FluidSystem::numPhases };
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
@@ -96,68 +79,164 @@ class EclOutputBlackOilModule : public BaseOutputModule<TypeTag>
     enum { gasCompIdx = FluidSystem::gasCompIdx };
     enum { oilCompIdx = FluidSystem::oilCompIdx };
 
-    typedef typename ParentType::ScalarBuffer ScalarBuffer;
+    typedef std::vector<Scalar> ScalarBuffer;
 
 public:
     EclOutputBlackOilModule(const Simulator& simulator)
-        : ParentType(simulator)
+        : simulator_(simulator)
     { }
 
-    /*!
-     * \brief Register all run-time parameters for the multi-phase VTK output
-     * module.
-     */
-    static void registerParameters()
-    {
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteSaturations,
-                             "Include the saturations of all fluid phases in the "
-                             "ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWritePressures,
-                             "Include the absolute pressures of all fluid phases in the "
-                             "ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor,
-                             "Include the gas dissolution factor of saturated oil in the "
-                             "ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilVaporizationFactor,
-                             "Include the oil vaporization factor of saturated gas in the "
-                             "ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor,
-                             "Include the gas formation volume factor in the "
-                             "ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilFormationVolumeFactor,
-                             "Include the oil formation volume factor of saturated oil "
-                             "in the ECL output files");
-        EWOMS_REGISTER_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure,
-                             "Include the saturation pressure of oil in the "
-                             "ECL output files");
-    }
-
     /*!
      * \brief Allocate memory for the scalar fields we would like to
-     *        write to the VTK file.
+     *        write to ECL output files
      */
-    void allocBuffers()
+    void allocBuffers(unsigned bufferSize, unsigned reportStepNum, const Opm::RestartConfig& restartConfig, const bool substep, const bool log)
     {
+
         if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
             return;
 
-        auto bufferType = ParentType::ElementBuffer;
+        std::map<std::string, int> rstKeywords = restartConfig.getRestartKeywords(reportStepNum);
-        if (saturationsOutput_()) {
+        for (auto& keyValue : rstKeywords) {
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
+            keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
-                this->resizeScalarBuffer_(saturation_[phaseIdx], bufferType);
         }
-        if (pressuresOutput_()) {
+
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                this->resizeScalarBuffer_(pressure_[phaseIdx], bufferType);
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            saturation_[phaseIdx].resize(bufferSize,0.0);
         }
-        if (gasDissolutionFactorOutput_())
+
-            this->resizeScalarBuffer_(gasDissolutionFactor_, bufferType);
+        oilPressure_.resize(bufferSize,0.0);
-        if (oilVaporizationFactorOutput_())
+        temperature_.resize(bufferSize,0.0);
-            this->resizeScalarBuffer_(oilVaporizationFactor_, bufferType);
+
-        if (gasFormationVolumeFactorOutput_())
+        // Output the same as legacy
-            this->resizeScalarBuffer_(gasFormationVolumeFactor_, bufferType);
+        // TODO: Only needed if DISGAS or VAPOIL
-        if (oilSaturationPressureOutput_())
+        if (true)
-            this->resizeScalarBuffer_(oilSaturationPressure_, bufferType);
+            rs_.resize(bufferSize,0.0);
+        if (true)
+            rv_.resize(bufferSize,0.0);
+
+        if (GET_PROP_VALUE(TypeTag, EnableSolvent)) {
+            sSol_.resize(bufferSize,0.0);
+        }
+        if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
+            cPolymer_.resize(bufferSize,0.0);
+        }
+
+        // Output the same as legacy
+        // TODO: Only needed if Vappars or hysteresis.
+        soMax_.resize(bufferSize,0.0);
+        pcSwMdcOw_.resize(bufferSize,0.0);
+        krnSwMdcOw_.resize(bufferSize,0.0);
+        pcSwMdcGo_.resize(bufferSize,0.0);
+        krnSwMdcGo_.resize(bufferSize,0.0);
+
+        // Only provide RESTART_AUXILIARY if it is asked for by the user
+        if (!restartConfig.getWriteRestartFile(reportStepNum) || substep)
+            return;
+
+        // Output the same as legacy
+        // TODO: Only needed if DISGAS or VAPOIL
+        if (FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx)) {
+            if (rstKeywords["RSSAT"] > 0) {
+                rstKeywords["RSSAT"] = 0;
+                gasDissolutionFactor_.resize(bufferSize,0.0);
+            }
+            if (rstKeywords["RVSAT"] > 0) {
+                rstKeywords["RVSAT"] = 0;
+                oilVaporizationFactor_.resize(bufferSize,0.0);
+            }
+        }
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["BW"] > 0)
+        {
+            rstKeywords["BW"] = 0;
+            invB_[waterPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(oilPhaseIdx) && rstKeywords["BO"] > 0)
+        {
+            rstKeywords["BO"] = 0;
+            invB_[oilPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(gasPhaseIdx) && rstKeywords["BG"] > 0)
+        {
+            rstKeywords["BG"] = 0;
+            invB_[gasPhaseIdx].resize(bufferSize,0.0);
+        }
+
+        if (rstKeywords["DEN"] > 0) {
+            rstKeywords["DEN"] = 0;
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (!FluidSystem::phaseIsActive(phaseIdx))
+                    continue;
+                density_[phaseIdx].resize(bufferSize,0.0);
+            }
+        }
+        const bool hasVWAT = (rstKeywords["VISC"] > 0) || (rstKeywords["VWAT"] > 0);
+        const bool hasVOIL = (rstKeywords["VISC"] > 0) || (rstKeywords["VOIL"] > 0);
+        const bool hasVGAS = (rstKeywords["VISC"] > 0) || (rstKeywords["VGAS"] > 0);
+        rstKeywords["VISC"] = 0;
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx) && hasVWAT)
+        {
+            rstKeywords["VWAT"] = 0;
+            viscosity_[waterPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(oilPhaseIdx) && hasVOIL > 0)
+        {
+            rstKeywords["VOIL"] = 0;
+            viscosity_[oilPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(gasPhaseIdx) && hasVGAS > 0)
+        {
+            rstKeywords["VGAS"] = 0;
+            viscosity_[gasPhaseIdx].resize(bufferSize,0.0);
+        }
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx) && rstKeywords["KRW"] > 0)
+        {
+            rstKeywords["KRW"] = 0;
+            relativePermeability_[waterPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(oilPhaseIdx) && rstKeywords["KRO"] > 0)
+        {
+            rstKeywords["KRO"] = 0;
+            relativePermeability_[oilPhaseIdx].resize(bufferSize,0.0);
+        }
+        if (FluidSystem::phaseIsActive(gasPhaseIdx) && rstKeywords["KRG"] > 0)
+        {
+            rstKeywords["KRG"] = 0;
+            relativePermeability_[gasPhaseIdx].resize(bufferSize,0.0);
+        }
+
+        if (rstKeywords["PBPD"] > 0)  {
+            rstKeywords["PBPD"] = 0;
+            bubblePointPressure_.resize(bufferSize,0.0);
+            dewPointPressure_.resize(bufferSize,0.0);
+        }
+
+        //Warn for any unhandled keyword
+        if (log) {
+            for (auto& keyValue : rstKeywords) {
+                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);
+                }
+            }
+        }
+
+        failedCellsPb_.clear();
+        failedCellsPd_.clear();
+
+        // Not supported in flow legacy
+        if (false)
+            saturatedOilFormationVolumeFactor_.resize(bufferSize,0.0);
+        if (false)
+            oilSaturationPressure_.resize(bufferSize,0.0);
     }
 
     /*!
@@ -166,158 +245,556 @@ public:
      */
     void processElement(const ElementContext& elemCtx)
     {
-        if (!EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput))
-            return;
-
         typedef Opm::MathToolbox<Evaluation> Toolbox;
 
         if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
             return;
 
         for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
-            const auto& fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();
+            const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
+            const auto& fs = intQuants.fluidState();
             typedef typename std::remove_const<typename std::remove_reference<decltype(fs)>::type>::type FluidState;
             unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
             unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
 
-            if (saturationsOutput_()) {
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (saturation_[phaseIdx].size() == 0)
-                    if (!FluidSystem::phaseIsActive(phaseIdx))
+                    continue;
-                        continue;
 
-                    saturation_[phaseIdx][globalDofIdx] = Toolbox::value(fs.saturation(phaseIdx));
+                saturation_[phaseIdx][globalDofIdx] = Toolbox::value(fs.saturation(phaseIdx));
-                    Opm::Valgrind::CheckDefined(saturation_[phaseIdx][globalDofIdx]);
+                Opm::Valgrind::CheckDefined(saturation_[phaseIdx][globalDofIdx]);
-                }
             }
-            if (pressuresOutput_()) {
-                for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                    if (!FluidSystem::phaseIsActive(phaseIdx))
-                        continue;
 
-                    pressure_[phaseIdx][globalDofIdx] = Toolbox::value(fs.pressure(phaseIdx));
+            if (oilPressure_.size() > 0) {
-                    Opm::Valgrind::CheckDefined(pressure_[phaseIdx][globalDofIdx]);
+                oilPressure_[globalDofIdx] = Toolbox::value(fs.pressure(oilPhaseIdx));
-                }
+                Opm::Valgrind::CheckDefined(oilPressure_[globalDofIdx]);
+
             }
-            if (gasDissolutionFactorOutput_()) {
+            if (gasDissolutionFactor_.size() > 0) {
-                Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx);
-                gasDissolutionFactor_[globalDofIdx] =
-                    FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx, SoMax);
-                Opm::Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]);
-            }
-            if (oilVaporizationFactorOutput_()) {
                 Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx);
                 gasDissolutionFactor_[globalDofIdx] =
                     FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx, SoMax);
                 Opm::Valgrind::CheckDefined(gasDissolutionFactor_[globalDofIdx]);
+
             }
-            if (gasFormationVolumeFactorOutput_()) {
+            if (oilVaporizationFactor_.size() > 0) {
+                Scalar SoMax = elemCtx.model().maxOilSaturation(globalDofIdx);
+                oilVaporizationFactor_[globalDofIdx] =
+                    FluidSystem::template saturatedDissolutionFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx, SoMax);
+                Opm::Valgrind::CheckDefined(oilVaporizationFactor_[globalDofIdx]);
+
+            }
+            if (gasFormationVolumeFactor_.size() > 0) {
                 gasFormationVolumeFactor_[globalDofIdx] =
                     1.0/FluidSystem::template inverseFormationVolumeFactor<FluidState, Scalar>(fs, gasPhaseIdx, pvtRegionIdx);
                 Opm::Valgrind::CheckDefined(gasFormationVolumeFactor_[globalDofIdx]);
+
             }
-            if (oilSaturationPressureOutput_()) {
+            if (saturatedOilFormationVolumeFactor_.size() > 0) {
+                saturatedOilFormationVolumeFactor_[globalDofIdx] =
+                    1.0/FluidSystem::template saturatedInverseFormationVolumeFactor<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
+                Opm::Valgrind::CheckDefined(saturatedOilFormationVolumeFactor_[globalDofIdx]);
+
+            }
+            if (oilSaturationPressure_.size() > 0) {
                 oilSaturationPressure_[globalDofIdx] =
                     FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
                 Opm::Valgrind::CheckDefined(oilSaturationPressure_[globalDofIdx]);
+
             }
+
+            if (rs_.size()) {
+                rs_[globalDofIdx] = Toolbox::value(fs.Rs());
+                Opm::Valgrind::CheckDefined(rs_[globalDofIdx]);
+            }
+
+            if (rv_.size()) {
+                rv_[globalDofIdx] = Toolbox::value(fs.Rv());
+                Opm::Valgrind::CheckDefined(rv_[globalDofIdx]);
+            }
+
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (invB_[phaseIdx].size() == 0)
+                    continue;
+
+                invB_[phaseIdx][globalDofIdx] = Toolbox::value(fs.invB(phaseIdx));
+                Opm::Valgrind::CheckDefined(invB_[phaseIdx][globalDofIdx]);
+            }
+
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (density_[phaseIdx].size() == 0)
+                    continue;
+
+                density_[phaseIdx][globalDofIdx] = Toolbox::value(fs.density(phaseIdx));
+                Opm::Valgrind::CheckDefined(density_[phaseIdx][globalDofIdx]);
+            }
+
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (viscosity_[phaseIdx].size() == 0)
+                    continue;
+
+                viscosity_[phaseIdx][globalDofIdx] = Toolbox::value(fs.viscosity(phaseIdx));
+                Opm::Valgrind::CheckDefined(viscosity_[phaseIdx][globalDofIdx]);
+            }
+
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                if (relativePermeability_[phaseIdx].size() == 0)
+                    continue;
+
+                relativePermeability_[phaseIdx][globalDofIdx] = Toolbox::value(intQuants.relativePermeability(phaseIdx));
+                Opm::Valgrind::CheckDefined(relativePermeability_[phaseIdx][globalDofIdx]);
+            }
+
+            if (sSol_.size() > 0) {
+                sSol_[globalDofIdx] = intQuants.solventSaturation().value();
+            }
+
+            if (cPolymer_.size() > 0) {
+                cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value();
+            }
+
+            if (bubblePointPressure_.size() > 0)
+            {
+                try {
+                    bubblePointPressure_[globalDofIdx] = Toolbox::value(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
+                }
+                catch (const Opm::NumericalProblem& e) {
+                    const auto globalIdx = elemCtx.simulator().gridManager().grid().globalCell()[globalDofIdx];
+                    failedCellsPb_.push_back(globalIdx);
+                }
+            }
+            if (dewPointPressure_.size() > 0)
+            {
+                try {
+                    dewPointPressure_[globalDofIdx] = Toolbox::value(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
+                }
+                catch (const Opm::NumericalProblem& e) {
+                    const auto globalIdx = elemCtx.simulator().gridManager().grid().globalCell()[globalDofIdx];
+                    failedCellsPd_.push_back(globalIdx);
+                }
+            }
+
+            if (soMax_.size() > 0)
+                soMax_[globalDofIdx] = elemCtx.simulator().model().maxOilSaturation(globalDofIdx);
+
+            const auto& matLawManager = elemCtx.simulator().problem().materialLawManager();
+            if (matLawManager->enableHysteresis()) {
+                if (pcSwMdcOw_.size() > 0 && krnSwMdcOw_.size() > 0) {
+                    matLawManager->oilWaterHysteresisParams(
+                                pcSwMdcOw_[globalDofIdx],
+                                krnSwMdcOw_[globalDofIdx],
+                                globalDofIdx);
+                }
+                if (pcSwMdcGo_.size() > 0 && krnSwMdcGo_.size() > 0) {
+                    matLawManager->gasOilHysteresisParams(
+                                pcSwMdcGo_[globalDofIdx],
+                                krnSwMdcGo_[globalDofIdx],
+                                globalDofIdx);
+                }
+            }
+
+            // hack to make the intial output of rs and rv Ecl compatible.
+            // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
+            // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
+            // rs and rv with the values computed in the initially.
+            // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
+            // This can be removed when ebos has 100% controll over output
+            if (elemCtx.simulator().episodeIndex() < 0 && FluidSystem::phaseIsActive(oilPhaseIdx) && FluidSystem::phaseIsActive(gasPhaseIdx) ) {
+
+                const auto& fs_initial = elemCtx.simulator().problem().initialFluidState(globalDofIdx);
+
+                // use initial rs and rv values
+                if (rv_.size() > 0)
+                    rv_[globalDofIdx] = fs_initial.Rv();
+
+                if (rs_.size() > 0)
+                    rs_[globalDofIdx] = fs_initial.Rs();
+
+                // re-compute the volume factors, viscosities and densities if asked for
+                if (density_[oilPhaseIdx].size() > 0)
+                    density_[oilPhaseIdx][globalDofIdx] = FluidSystem::density(fs_initial,
+                                                                                    oilPhaseIdx,
+                                                                                    intQuants.pvtRegionIndex());
+                if (density_[gasPhaseIdx].size() > 0)
+                    density_[gasPhaseIdx][globalDofIdx] = FluidSystem::density(fs_initial,
+                                                                                    gasPhaseIdx,
+                                                                                    intQuants.pvtRegionIndex());
+
+                if (invB_[oilPhaseIdx].size() > 0)
+                    invB_[oilPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fs_initial,
+                                                                                                      oilPhaseIdx,
+                                                                                                      intQuants.pvtRegionIndex());
+                if (invB_[gasPhaseIdx].size() > 0)
+                    invB_[gasPhaseIdx][globalDofIdx] = FluidSystem::inverseFormationVolumeFactor(fs_initial,
+                                                                                                      gasPhaseIdx,
+                                                                                                      intQuants.pvtRegionIndex());
+                if (viscosity_[oilPhaseIdx].size() > 0)
+                    viscosity_[oilPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fs_initial,
+                                                                                        oilPhaseIdx,
+                                                                                        intQuants.pvtRegionIndex());
+                if (viscosity_[gasPhaseIdx].size() > 0)
+                    viscosity_[gasPhaseIdx][globalDofIdx] = FluidSystem::viscosity(fs_initial,
+                                                                                        gasPhaseIdx,
+                                                                                        intQuants.pvtRegionIndex());
+            }
+
+        }
+
+    }
+
+
+    void outputErrorLog()
+    {
+        const size_t maxNumCellsFaillog = 20;
+
+        int pbSize = failedCellsPb_.size(), pd_size = failedCellsPd_.size();
+        std::vector<int> displPb, displPd, recvLenPb, recvLenPd;
+        const auto& comm = simulator_.gridView().comm();
+
+        if ( isIORank_() )
+        {
+            displPb.resize(comm.size()+1, 0);
+            displPd.resize(comm.size()+1, 0);
+            recvLenPb.resize(comm.size());
+            recvLenPd.resize(comm.size());
+        }
+
+        comm.gather(&pbSize, recvLenPb.data(), 1, 0);
+        comm.gather(&pd_size, recvLenPd.data(), 1, 0);
+        std::partial_sum(recvLenPb.begin(), recvLenPb.end(), displPb.begin()+1);
+        std::partial_sum(recvLenPd.begin(), recvLenPd.end(), displPd.begin()+1);
+        std::vector<int> globalFailedCellsPb, globalFailedCellsPd;
+
+        if ( isIORank_() )
+        {
+            globalFailedCellsPb.resize(displPb.back());
+            globalFailedCellsPd.resize(displPd.back());
+        }
+
+        comm.gatherv(failedCellsPb_.data(), static_cast<int>(failedCellsPb_.size()),
+                     globalFailedCellsPb.data(), recvLenPb.data(),
+                     displPb.data(), 0);
+        comm.gatherv(failedCellsPd_.data(), static_cast<int>(failedCellsPd_.size()),
+                     globalFailedCellsPd.data(),  recvLenPd.data(),
+                     displPd.data(), 0);
+        std::sort(globalFailedCellsPb.begin(), globalFailedCellsPb.end());
+        std::sort(globalFailedCellsPd.begin(), globalFailedCellsPd.end());
+
+        if (globalFailedCellsPb.size() > 0) {
+            std::stringstream errlog;
+            errlog << "Finding the bubble point pressure failed for " << globalFailedCellsPb.size() << " cells [";
+            errlog << globalFailedCellsPb[0];
+            const size_t max_elems = std::min(maxNumCellsFaillog, globalFailedCellsPb.size());
+            for (size_t i = 1; i < max_elems; ++i) {
+                errlog << ", " << globalFailedCellsPb[i];
+            }
+            if (globalFailedCellsPb.size() > maxNumCellsFaillog) {
+                errlog << ", ...";
+            }
+            errlog << "]";
+            Opm::OpmLog::warning("Bubble point numerical problem", errlog.str());
+        }
+        if (globalFailedCellsPd.size() > 0) {
+            std::stringstream errlog;
+            errlog << "Finding the dew point pressure failed for " << globalFailedCellsPd.size() << " cells [";
+            errlog << globalFailedCellsPd[0];
+            const size_t max_elems = std::min(maxNumCellsFaillog, globalFailedCellsPd.size());
+            for (size_t i = 1; i < max_elems; ++i) {
+                errlog << ", " << globalFailedCellsPd[i];
+            }
+            if (globalFailedCellsPd.size() > maxNumCellsFaillog) {
+                errlog << ", ...";
+            }
+            errlog << "]";
+            Opm::OpmLog::warning("Dew point numerical problem", errlog.str());
         }
     }
 
     /*!
-     * \brief Add all buffers to the VTK output writer.
+     * \brief Move all buffers to data::Solution.
      */
-    void commitBuffers(BaseOutputWriter& writer)
+    void assignToSolution(Opm::data::Solution& sol)
     {
         if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
             return;
 
-        if (!dynamic_cast<EclWriter*>(&writer))
+        if ( oilPressure_.size() > 0 ) {
-            return; // this module only consideres ecl writers...
+            sol.insert( "PRESSURE", Opm::UnitSystem::measure::pressure, std::move(oilPressure_), Opm::data::TargetType::RESTART_SOLUTION);
+        }
 
-        typedef EclDeckUnits<TypeTag> DeckUnits;
+        if ( temperature_.size() > 0 ) {
-        const DeckUnits& deckUnits = this->simulator_.problem().deckUnits();
+            sol.insert( "TEMP", Opm::UnitSystem::measure::temperature, std::move(temperature_), Opm::data::TargetType::RESTART_SOLUTION);
+        }
 
-        typename ParentType::BufferType bufferType = ParentType::ElementBuffer;
+        if( FluidSystem::phaseIsActive(waterPhaseIdx) && saturation_[waterPhaseIdx].size() > 0 ) {
-        if (pressuresOutput_()) {
+            sol.insert( "SWAT", Opm::UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION );
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
+        }
-                deckUnits.siToDeck(pressure_[phaseIdx], DeckUnits::pressure);
+        if( FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size() > 0) {
+            sol.insert( "SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION );
+        }
 
-            this->commitScalarBuffer_(writer, "PRESSURE", pressure_[oilPhaseIdx], bufferType);
+        if ( gasDissolutionFactor_.size() > 0 ) {
-            this->commitScalarBuffer_(writer, "PGAS", pressure_[gasPhaseIdx], bufferType);
+            sol.insert( "RSSAT", Opm::UnitSystem::measure::gas_oil_ratio, std::move(gasDissolutionFactor_), Opm::data::TargetType::RESTART_AUXILIARY );
-            this->commitScalarBuffer_(writer, "PWAT", pressure_[waterPhaseIdx], bufferType);
-        }
-        if (saturationsOutput_()) {
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx)
-                deckUnits.siToDeck(saturation_[phaseIdx], DeckUnits::saturation);
 
-            this->commitScalarBuffer_(writer, "SWAT", saturation_[waterPhaseIdx], bufferType);
-            this->commitScalarBuffer_(writer, "SGAS", saturation_[gasPhaseIdx], bufferType);
-            // the oil saturation is _NOT_ written to disk. Instead, it is calculated by
-            // the visualization tool. Wondering why is probably a waste of time...
         }
-        if (gasDissolutionFactorOutput_()) {
+        if ( oilVaporizationFactor_.size() > 0 ) {
-            deckUnits.siToDeck(gasDissolutionFactor_, DeckUnits::gasDissolutionFactor);
+            sol.insert( "RVSAT", Opm::UnitSystem::measure::oil_gas_ratio, std::move(oilVaporizationFactor_) , Opm::data::TargetType::RESTART_AUXILIARY );
-            this->commitScalarBuffer_(writer, "RS", gasDissolutionFactor_, bufferType);
         }
-        if (oilVaporizationFactorOutput_()) {
+        if ( rs_.size() > 0 ) {
-            deckUnits.siToDeck(oilVaporizationFactor_, DeckUnits::oilVaporizationFactor);
+            sol.insert( "RS", Opm::UnitSystem::measure::gas_oil_ratio, std::move(rs_), Opm::data::TargetType::RESTART_SOLUTION );
-            this->commitScalarBuffer_(writer, "RV", oilVaporizationFactor_, bufferType);
+
         }
-        if (gasFormationVolumeFactorOutput_()) {
+        if (rv_.size() > 0 ) {
-            // no unit conversion required
+            sol.insert( "RV", Opm::UnitSystem::measure::oil_gas_ratio, std::move(rv_) , Opm::data::TargetType::RESTART_SOLUTION );
-            this->commitScalarBuffer_(writer, "BG", gasFormationVolumeFactor_, bufferType);
         }
-        if (oilSaturationPressureOutput_()) {
+        if (invB_[waterPhaseIdx].size() > 0 ) {
-            deckUnits.siToDeck(oilSaturationPressure_, DeckUnits::pressure);
+            sol.insert( "1OVERBW", Opm::UnitSystem::measure::water_inverse_formation_volume_factor, std::move(invB_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
-            this->commitScalarBuffer_(writer, "PSAT", oilSaturationPressure_);
+        }
+        if (invB_[oilPhaseIdx].size() > 0 ) {
+            sol.insert( "1OVERBO", Opm::UnitSystem::measure::oil_inverse_formation_volume_factor, std::move(invB_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+        if (invB_[gasPhaseIdx].size() > 0 ) {
+            sol.insert( "1OVERBG", Opm::UnitSystem::measure::gas_inverse_formation_volume_factor, std::move(invB_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+
+        if (density_[waterPhaseIdx].size() > 0 ) {
+            sol.insert( "WAT_DEN", Opm::UnitSystem::measure::density, std::move(density_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+        if (density_[oilPhaseIdx].size() > 0 ) {
+            sol.insert( "OIL_DEN", Opm::UnitSystem::measure::density, std::move(density_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+        if (density_[gasPhaseIdx].size() > 0 ) {
+            sol.insert( "GAS_DEN", Opm::UnitSystem::measure::density, std::move(density_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+
+        if (viscosity_[waterPhaseIdx].size() > 0 ) {
+            sol.insert( "WAT_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
+        }
+        if (viscosity_[oilPhaseIdx].size() > 0 ) {
+            sol.insert( "OIL_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
+        }
+        if (viscosity_[gasPhaseIdx].size() > 0 ) {
+            sol.insert( "GAS_VISC", Opm::UnitSystem::measure::viscosity, std::move(viscosity_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY );
+        }
+
+        if (relativePermeability_[waterPhaseIdx].size() > 0) {
+            sol.insert( "WATKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[waterPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
+        }
+        if (relativePermeability_[oilPhaseIdx].size() > 0 ) {
+            sol.insert( "OILKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[oilPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
+        }
+        if (relativePermeability_[gasPhaseIdx].size() > 0 ) {
+            sol.insert( "GASKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
+        }
+
+        if (pcSwMdcOw_.size() > 0 )
+            sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, std::move(pcSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+        if (krnSwMdcOw_.size() > 0)
+            sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, std::move(krnSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+        if (pcSwMdcGo_.size() > 0)
+            sol.insert ("PCSWM_GO", Opm::UnitSystem::measure::identity, std::move(pcSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+        if (krnSwMdcGo_.size() > 0)
+            sol.insert ("KRNSW_GO", Opm::UnitSystem::measure::identity, std::move(krnSwMdcGo_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+        if (soMax_.size() > 0)
+            sol.insert ("SOMAX", Opm::UnitSystem::measure::identity, std::move(soMax_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (sSol_.size() > 0)
+            sol.insert ("SSOL", Opm::UnitSystem::measure::identity, std::move(sSol_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (cPolymer_.size() > 0)
+            sol.insert ("POLYMER", Opm::UnitSystem::measure::identity, std::move(cPolymer_), Opm::data::TargetType::RESTART_SOLUTION);
+
+        if (dewPointPressure_.size() > 0)
+            sol.insert ("PDEW", Opm::UnitSystem::measure::pressure, std::move(dewPointPressure_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+        if (bubblePointPressure_.size() > 0)
+            sol.insert ("PBUB", Opm::UnitSystem::measure::pressure, std::move(bubblePointPressure_), Opm::data::TargetType::RESTART_AUXILIARY);
+
+}
+
+    void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex) 
+    {
+
+        Scalar so = 1.0;
+        if( sol.has( "SWAT" ) ) {
+            saturation_[waterPhaseIdx][elemIdx] = sol.data("SWAT")[globalDofIndex];
+            so -= sol.data("SWAT")[globalDofIndex];
+        }
+
+        if( sol.has( "SGAS" ) ) {
+            saturation_[gasPhaseIdx][elemIdx] = sol.data("SGAS")[globalDofIndex];
+            so -= sol.data("SGAS")[globalDofIndex];
+        }
+
+        saturation_[oilPhaseIdx][elemIdx] = so;
+
+        if( sol.has( "PRESSURE" ) ) {
+            oilPressure_[elemIdx] = sol.data( "PRESSURE" )[globalDofIndex];
+        }
+
+        if( sol.has( "TEMP" ) ) {
+            temperature_[elemIdx] = sol.data( "TEMP" )[globalDofIndex];
+        }
+
+        if( sol.has( "RS" ) ) {
+            rs_[elemIdx] = sol.data("RS")[globalDofIndex];
+        }
+
+        if( sol.has( "RV" ) ) {
+            rv_[elemIdx] = sol.data("RV")[globalDofIndex];
+        }
+
+        if ( sol.has( "SSOL" ) ) {
+            sSol_[elemIdx] = sol.data("SSOL")[globalDofIndex];
+        }
+
+        if ( sol.has("POLYMER" ) ) {
+            cPolymer_[elemIdx] = sol.data("POLYMER")[globalDofIndex];
+        }
+
+        if ( sol.has("SOMAX" ) ) {
+            soMax_[elemIdx] = sol.data("SOMAX")[globalDofIndex];
+        }
+
+        if ( sol.has("PCSWM_OW" ) ) {
+            pcSwMdcOw_[elemIdx] = sol.data("PCSWM_OW")[globalDofIndex];
+        }
+
+        if ( sol.has("KRNSW_OW" ) ) {
+            krnSwMdcOw_[elemIdx] = sol.data("KRNSW_OW")[globalDofIndex];
+        }
+
+        if ( sol.has("PCSWM_GO" ) ) {
+            pcSwMdcGo_[elemIdx] = sol.data("PCSWM_GO")[globalDofIndex];
+        }
+
+        if ( sol.has("KRNSW_GO" ) ) {
+            krnSwMdcGo_[elemIdx] = sol.data("KRNSW_GO")[globalDofIndex];
         }
     }
 
+
+    template <class FluidState>
+    void assignToFluidState(FluidState& fs, unsigned elemIdx) const 
+    {
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (saturation_[phaseIdx].size() == 0)
+                continue;
+
+            fs.setSaturation(phaseIdx, saturation_[phaseIdx][elemIdx]);
+        }
+
+        if (oilPressure_.size() > 0) {
+
+            // this assumes that capillary pressures only depend on the phase saturations
+            // and possibly on temperature. (this is always the case for ECL problems.)
+            Dune::FieldVector< Scalar, numPhases > pc( 0 );
+            const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
+            MaterialLaw::capillaryPressures(pc, matParams, fs);
+            Opm::Valgrind::CheckDefined(oilPressure_[elemIdx]);
+            Opm::Valgrind::CheckDefined(pc);
+            assert(FluidSystem::phaseIsActive(oilPhaseIdx));
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                if (!FluidSystem::phaseIsActive(phaseIdx))
+                    continue;
+
+                fs.setPressure(phaseIdx, oilPressure_[elemIdx] + (pc[phaseIdx] - pc[oilPhaseIdx]));
+            }
+        }
+
+        if (temperature_.size() > 0) {
+            fs.setTemperature( temperature_[elemIdx]);
+        }
+
+        if (rs_.size() > 0) {
+           fs.setRs(rs_[elemIdx]);
+
+        }
+        if (rv_.size() > 0) {
+           fs.setRv(rv_[elemIdx]);
+        }
+    }
+
+    void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const {
+
+        if (soMax_.size() > 0)
+            simulator.model().setMaxOilSaturation(soMax_[elemIdx], elemIdx);
+
+        if (simulator.problem().materialLawManager()->enableHysteresis()) {
+            auto matLawManager = simulator.problem().materialLawManager();
+
+            if (pcSwMdcOw_.size() > 0 && krnSwMdcOw_.size() > 0) {
+                matLawManager->setOilWaterHysteresisParams(
+                            pcSwMdcOw_[elemIdx],
+                            krnSwMdcOw_[elemIdx],
+                            elemIdx);
+            }
+            if (pcSwMdcGo_.size() > 0 && krnSwMdcGo_.size() > 0) {
+                matLawManager->setGasOilHysteresisParams(
+                            pcSwMdcGo_[elemIdx],
+                            krnSwMdcGo_[elemIdx],
+                            elemIdx);
+            }
+        }
+
+    }
+
+    Scalar getSolventSaturation(unsigned elemIdx) const {
+        if(sSol_.size() > 0)
+            return sSol_[elemIdx];
+
+        return 0;
+    }
+
+    Scalar getPolymerConcentration(unsigned elemIdx) const {
+        if(cPolymer_.size() > 0)
+            return cPolymer_[elemIdx];
+
+        return 0;
+    }
+
+
 private:
-    static bool saturationsOutput_()
-    { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteSaturations); }
 
-    static bool pressuresOutput_()
+    bool isIORank_() const
-    { return EWOMS_GET_PARAM(TypeTag, bool, EclOutputWritePressures); }
-
-    static bool gasDissolutionFactorOutput_()
     {
-        return
+        const auto& comm = simulator_.gridView().comm();
-            FluidSystem::phaseIsActive(oilPhaseIdx) &&
+        return comm.rank() == 0;
-            FluidSystem::phaseIsActive(gasPhaseIdx) &&
-            EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasDissolutionFactor);
     }
 
-    static bool gasFormationVolumeFactorOutput_()
+    const Simulator& simulator_;
-    {
-        return
-            FluidSystem::phaseIsActive(oilPhaseIdx) &&
-            FluidSystem::phaseIsActive(gasPhaseIdx) &&
-            EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteGasFormationVolumeFactor);
-    }
-
-    static bool oilVaporizationFactorOutput_()
-    {
-        return
-            FluidSystem::phaseIsActive(oilPhaseIdx) &&
-            FluidSystem::phaseIsActive(gasPhaseIdx) &&
-            EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilVaporizationFactor);
-    }
-
-    static bool oilSaturationPressureOutput_()
-    {
-        return
-            FluidSystem::phaseIsActive(oilPhaseIdx) &&
-            FluidSystem::phaseIsActive(gasPhaseIdx) &&
-            EWOMS_GET_PARAM(TypeTag, bool, EclOutputWriteOilSaturationPressure);
-    }
 
     ScalarBuffer saturation_[numPhases];
-    ScalarBuffer pressure_[numPhases];
+    ScalarBuffer oilPressure_;
+    ScalarBuffer temperature_;
     ScalarBuffer gasDissolutionFactor_;
     ScalarBuffer oilVaporizationFactor_;
     ScalarBuffer gasFormationVolumeFactor_;
+    ScalarBuffer saturatedOilFormationVolumeFactor_;
     ScalarBuffer oilSaturationPressure_;
+    ScalarBuffer rs_;
+    ScalarBuffer rv_;
+    ScalarBuffer invB_[numPhases];
+    ScalarBuffer density_[numPhases];
+    ScalarBuffer viscosity_[numPhases];
+    ScalarBuffer relativePermeability_[numPhases];
+    ScalarBuffer sSol_;
+    ScalarBuffer cPolymer_;
+    ScalarBuffer soMax_;
+    ScalarBuffer pcSwMdcOw_;
+    ScalarBuffer krnSwMdcOw_;
+    ScalarBuffer pcSwMdcGo_;
+    ScalarBuffer krnSwMdcGo_;
+    ScalarBuffer bubblePointPressure_;
+    ScalarBuffer dewPointPressure_;
+    std::vector<int> failedCellsPb_;
+    std::vector<int> failedCellsPd_;
+
 };
 } // namespace Ewoms
 

ebos/eclproblem.hh

@@ -54,13 +54,11 @@
 #include "eclwellmanager.hh"
 #include "eclequilinitializer.hh"
 #include "eclwriter.hh"
-#include "eclsummarywriter.hh"
 #include "ecloutputblackoilmodule.hh"
 #include "ecltransmissibility.hh"
 #include "eclthresholdpressure.hh"
 #include "ecldummygradientcalculator.hh"
 #include "eclfluxmodule.hh"
-#include "ecldeckunits.hh"
 
 #include <ewoms/common/pffgridvector.hh>
 #include <ewoms/models/blackoil/blackoilmodel.hh>
@@ -80,6 +78,7 @@
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/Eqldims.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/common/Exceptions.hpp>
 
@@ -87,6 +86,8 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
+#include <opm/output/eclipse/EclipseIO.hpp>
+
 #include <boost/date_time.hpp>
 
 #include <vector>
@@ -205,9 +206,6 @@ SET_BOOL_PROP(EclBaseProblem, EnableVtkOutput, false);
 // ... but enable the ECL output by default
 SET_BOOL_PROP(EclBaseProblem, EnableEclOutput, true);
 
-// also enable the summary output.
-SET_BOOL_PROP(EclBaseProblem, EnableEclSummaryOutput, true);
-
 // the cache for intensive quantities can be used for ECL problems and also yields a
 // decent speedup...
 SET_BOOL_PROP(EclBaseProblem, EnableIntensiveQuantityCache, true);
@@ -292,11 +290,13 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
             /*enableTemperature=*/true> InitialFluidState;
 
     typedef Opm::MathToolbox<Evaluation> Toolbox;
-    typedef Ewoms::EclSummaryWriter<TypeTag> EclSummaryWriter;
     typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
 
     typedef EclWriter<TypeTag> EclWriterType;
 
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+
     struct RockParams {
         Scalar referencePressure;
         Scalar compressibility;
@@ -310,8 +310,6 @@ public:
     {
         ParentType::registerParameters();
 
-        Ewoms::EclOutputBlackOilModule<TypeTag>::registerParameters();
-
         EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWriteAllSolutions,
                              "Write all solutions to disk instead of only the ones for the "
                              "report steps");
@@ -330,20 +328,17 @@ public:
         , transmissibilities_(simulator.gridManager())
         , thresholdPressures_(simulator)
         , wellManager_(simulator)
-        , deckUnits_(simulator)
         , eclWriter_( EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput)
                         ? new EclWriterType(simulator) : nullptr )
-        , summaryWriter_(simulator)
         , pffDofData_(simulator.gridView(), this->elementMapper())
     {
-        // add the output module for the Ecl binary output
+        // Tell the extra modules to initialize its internal data structures
-        simulator.model().addOutputModule(new Ewoms::EclOutputBlackOilModule<TypeTag>(simulator));
-
-        // Tell the solvent module to initialize its internal data structures
         const auto& gridManager = simulator.gridManager();
         SolventModule::initFromDeck(gridManager.deck(), gridManager.eclState());
         PolymerModule::initFromDeck(gridManager.deck(), gridManager.eclState());
 
+        // Hack to compute the initial thpressure values for restarts
+        restartApplied = false;
     }
 
     /*!
@@ -502,14 +497,6 @@ public:
     {
         if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
             wellManager_.beginTimeStep();
-
-            // this is a little hack to write the initial condition, which we need to do
-            // before the first time step has finished.
-            static bool initialWritten = false;
-            if (this->simulator().episodeIndex() == 0 && !initialWritten) {
-                summaryWriter_.write(wellManager_, /*isInitial=*/true);
-                initialWritten = true;
-            }
         }
     }
 
@@ -548,9 +535,6 @@ public:
 
         if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
             wellManager_.endTimeStep();
-
-            // write the summary information after each time step
-            summaryWriter_.write(wellManager_);
         }
 
         // we no longer need the initial soluiton
@@ -613,28 +597,41 @@ public:
      */
     void writeOutput(bool verbose = true)
     {
-        // calculate the time _after_ the time was updated
         Scalar t = this->simulator().time() + this->simulator().timeStepSize();
 
-        // prepare the ECL and the VTK writers
+        Opm::data::Wells dw;
-        if ( eclWriter_ )
+        if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
-            eclWriter_->beginWrite(t);
+            using rt = Opm::data::Rates::opt;
+            for (unsigned wellIdx = 0; wellIdx < wellManager_.numWells(); ++wellIdx) {
+                const auto& well = wellManager_.well(wellIdx);
+                auto& wellOut = dw[ well->name() ];
 
+                wellOut.bhp = well->bottomHolePressure();
+                wellOut.thp = well->tubingHeadPressure();
+                wellOut.temperature = 0;
+                wellOut.rates.set( rt::wat, well->surfaceRate(waterPhaseIdx) );
+                wellOut.rates.set( rt::oil, well->surfaceRate(oilPhaseIdx) );
+                wellOut.rates.set( rt::gas, well->surfaceRate(gasPhaseIdx) );
+            }
+        }
+        Scalar totalSolverTime = 0.0;
+        Scalar nextstep = this->simulator().timeStepSize();
+        Opm::data::Solution fip;
+        writeOutput(dw, t, false, totalSolverTime, nextstep, fip, verbose);
+    }
+
+    void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip, bool verbose = true)
+    {
         // use the generic code to prepare the output fields and to
         // write the desired VTK files.
         ParentType::writeOutput(verbose);
 
+        // output using eclWriter if enabled
         if ( eclWriter_ ) {
-            this->model().appendOutputFields(*eclWriter_);
+            eclWriter_->writeOutput(dw, t, substep, totalSolverTime, nextstep, fip);
-            eclWriter_->endWrite();
         }
-    }
 
-    /*!
+    }
-     * \brief Returns the object which converts between SI and deck units.
-     */
-    const EclDeckUnits<TypeTag>& deckUnits() const
-    { return deckUnits_; }
 
     /*!
      * \copydoc FvBaseMultiPhaseProblem::intrinsicPermeability
@@ -960,6 +957,7 @@ public:
      */
     void initialSolutionApplied()
     {
+
         if (!GET_PROP_VALUE(TypeTag, DisableWells)) {
             // initialize the wells. Note that this needs to be done after initializing the
             // intrinsic permeabilities and the after applying the initial solution because
@@ -967,11 +965,24 @@ public:
             wellManager_.init(this->simulator().gridManager().eclState(), this->simulator().gridManager().schedule());
         }
 
+        // the initialSolutionApplied is called recursively by readEclRestartSolution_()
+        // in order to setup the inital threshold pressures correctly
+        if (restartApplied)
+            return;
+
         // let the object for threshold pressures initialize itself. this is done only at
         // this point, because determining the threshold pressures may require to access
         // the initial solution.
         thresholdPressures_.finishInit();
 
+        const auto& eclState = this->simulator().gridManager().eclState();
+        const auto& initconfig = eclState.getInitConfig();
+        if(initconfig.restartRequested()) {
+            restartApplied = true;
+            this->simulator().setEpisodeIndex(initconfig.getRestartStep());
+            readEclRestartSolution_();
+        }
+
         // release the memory of the EQUIL grid since it's no longer needed after this point
         this->simulator().gridManager().releaseEquilGrid();
     }
@@ -1013,6 +1024,13 @@ public:
         return initialFluidStates_[globalDofIdx];
     }
 
+    void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO) {
+        eclWriter_->setEclIO(std::move(eclIO));
+    }
+
+    const Opm::EclipseIO& eclIO() const
+    {return eclWriter_->eclIO();}
+
 private:
     Scalar cellCenterDepth( const Element& element ) const
     {
@@ -1205,14 +1223,15 @@ private:
     void readInitialCondition_()
     {
         const auto& gridManager = this->simulator().gridManager();
-        const auto& deck = gridManager.deck();
 
+        const auto& deck = gridManager.deck();
         if (!deck.hasKeyword("EQUIL"))
             readExplicitInitialCondition_();
         else
             readEquilInitialCondition_();
 
         readBlackoilExtentionsInitialConditions_();
+
     }
 
 
@@ -1235,6 +1254,35 @@ private:
         }
     }
 
+    void readEclRestartSolution_()
+    {
+        // since the EquilInitializer provides fluid states that are consistent with the
+        // black-oil model, we can use naive instead of mass conservative determination
+        // of the primary variables.
+        useMassConservativeInitialCondition_ = false;
+
+        eclWriter_->restartBegin();
+
+        size_t numElems = this->model().numGridDof();
+        initialFluidStates_.resize(numElems);
+        if (enableSolvent)
+            solventSaturation_.resize(numElems,0.0);
+
+        if (enablePolymer)
+            polymerConcentration_.resize(numElems,0.0);
+
+        for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
+            auto& elemFluidState = initialFluidStates_[elemIdx];
+            eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
+            eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
+            if (enableSolvent)
+                 solventSaturation_[elemIdx] = eclWriter_->eclOutputModule().getSolventSaturation(elemIdx);
+            if (enablePolymer)
+                 polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
+        }
+        this->model().applyInitialSolution();
+    }
+
     void readExplicitInitialCondition_()
     {
         const auto& gridManager = this->simulator().gridManager();
@@ -1572,12 +1620,12 @@ private:
 
     EclWellManager<TypeTag> wellManager_;
 
-    EclDeckUnits<TypeTag> deckUnits_;
-
     std::unique_ptr< EclWriterType > eclWriter_;
-    EclSummaryWriter summaryWriter_;
 
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
+
+    bool restartApplied;
+
 };
 } // namespace Ewoms
 

ebos/eclwriter.hh

@@ -30,11 +30,12 @@
 
 #include <opm/material/densead/Evaluation.hpp>
 
-#include "ertwrappers.hh"
 #include "collecttoiorank.hh"
+#include "ecloutputblackoilmodule.hh"
 
 #include <ewoms/disc/ecfv/ecfvdiscretization.hh>
 #include <ewoms/io/baseoutputwriter.hh>
+#include <opm/output/eclipse/EclipseIO.hpp>
 
 #include <opm/common/Valgrind.hpp>
 #include <opm/common/ErrorMacros.hpp>
@@ -58,302 +59,170 @@ NEW_PROP_TAG(EnableEclOutput);
 template <class TypeTag>
 class EclWriter;
 
-template <class TypeTag, class GridManagerType>
-class EclWriterHelper
-{
-    friend class EclWriter<TypeTag>;
-
-    static void writeHeaders_(EclWriter<TypeTag>& writer)
-    {
-        typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
-        if (!std::is_same<Discretization, Ewoms::EcfvDiscretization<TypeTag> >::value)
-            OPM_THROW(std::logic_error,
-                      "Ecl binary output only works for the element centered "
-                      "finite volume discretization.");
-
-#if ! HAVE_ERT
-        OPM_THROW(std::logic_error,
-                  "Ecl binary output requires the ERT libraries");
-#else
-        // set the index of the first time step written to 0...
-        writer.reportStepIdx_ = 0;
-
-        char* egridRawFileName = ecl_util_alloc_filename(/*outputDir=*/"./",
-                                                         writer.caseName().c_str(),
-                                                         ECL_EGRID_FILE,
-                                                         /*formatted=*/false, // -> write binary output
-                                                         writer.reportStepIdx_);
-        std::string egridFileName(egridRawFileName);
-        std::free(egridRawFileName);
-
-        ErtGrid ertGrid(writer.simulator_.gridManager().eclState().getInputGrid(),
-                        writer.simulator_.gridManager().grid(),
-                        writer.simulator_.gridManager().cartesianIndexMapper(),
-                        writer.simulator_.problem().deckUnits());
-        ertGrid.write(egridFileName, writer.reportStepIdx_);
-#endif
-    }
-};
 
 /*!
  * \ingroup EclBlackOilSimulator
  *
- * \brief Implements writing Ecl binary output files.
+ * \brief Collects necessary output values and pass it to opm-output.
  *
  * Caveats:
- * - For this class to do do anything meaningful, you need to have the
+ * - For this class to do do anything meaningful, you will have to
- *   ERT libraries with development headers installed and the ERT
+ *   have the OPM module opm-output.
- *   build system test must pass sucessfully.
  * - The only DUNE grid which is currently supported is Dune::CpGrid
- *   from the OPM module "opm-core". Using another grid won't
+ *   from the OPM module "opm-grid". Using another grid won't
  *   fail at compile time but you will provoke a fatal exception as
  *   soon as you try to write an ECL output file.
  * - This class requires to use the black oil model with the element
  *   centered finite volume discretization.
- * - MPI-parallel computations are not (yet?) supported.
  */
 template <class TypeTag>
-class EclWriter : public BaseOutputWriter
+class EclWriter
 {
-    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
     typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
     typedef typename GET_PROP_TYPE(TypeTag, GridManager) GridManager;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
 
     typedef CollectDataToIORank< GridManager > CollectDataToIORankType;
 
-    typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
+    typedef std::vector<Scalar> ScalarBuffer;
-    typedef BaseOutputWriter::VectorBuffer VectorBuffer;
-    typedef BaseOutputWriter::TensorBuffer TensorBuffer;
 
-    friend class EclWriterHelper<TypeTag, GridManager>;
 
 public:
     EclWriter(const Simulator& simulator)
         : simulator_(simulator)
-        , gridView_(simulator_.gridView())
+        , eclOutputModule_(simulator)
-#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
-        , elementMapper_(gridView_, Dune::mcmgElementLayout())
-        , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
-#else
-        , elementMapper_(gridView_)
-        , vertexMapper_(gridView_)
-#endif
         , collectToIORank_( simulator_.gridManager() )
     {
-        reportStepIdx_ = 0;
+        Grid globalGrid = simulator_.gridManager().grid();
+        globalGrid.switchToGlobalView();
+        eclIO_.reset(new Opm::EclipseIO(simulator_.gridManager().eclState(),
+                                        Opm::UgGridHelpers::createEclipseGrid( globalGrid , simulator_.gridManager().eclState().getInputGrid() ),
+                                        simulator_.gridManager().schedule(),
+                                        simulator_.gridManager().summaryConfig()));
     }
 
     ~EclWriter()
     { }
 
-    /*!
+    void setEclIO(std::unique_ptr<Opm::EclipseIO>&& eclIO) {
-     * \brief Returns the name of the simulation.
+        eclIO_ = std::move(eclIO);
-     *
+    }
-     * This is the prefix of the files written to disk.
+
-     */
+    const Opm::EclipseIO& eclIO() const
-    std::string caseName() const
+    {return *eclIO_;}
-    { return boost::to_upper_copy(simulator_.problem().name()); }
 
     /*!
-     * \brief Updates the internal data structures after mesh
+     * \brief collect and pass data and pass it to eclIO writer
-     *        refinement.
-     *
-     * If the grid changes between two calls of beginWrite(), this
-     * method _must_ be called before the second beginWrite()!
      */
-    void gridChanged()
+    void writeOutput(const Opm::data::Wells& dw, Scalar t, bool substep, Scalar totalSolverTime, Scalar nextstep, const Opm::data::Solution& fip)
     {
-        elementMapper_.update();
-        vertexMapper_.update();
-    }
 
-    /*!
+        #if !HAVE_OPM_OUTPUT
-     * \brief Called whenever a new time step must be written.
+                OPM_THROW(std::runtime_error,
-     */
+                          "Opm-output must be available to write ECL output!");
-    void beginWrite(double t OPM_UNUSED)
+        #else
-    {
-        if (enableEclOutput_() && reportStepIdx_ == 0 && collectToIORank_.isIORank() )
-            EclWriterHelper<TypeTag, GridManager>::writeHeaders_(*this);
-    }
 
-    /*
+        int episodeIdx = simulator_.episodeIndex() + 1;
-     * \brief Add a vertex-centered scalar field to the output.
+        const auto& gridView = simulator_.gridManager().gridView();
-     *
+        int numElements = gridView.size(/*codim=*/0);
-     * For the EclWriter, this method is a no-op which throws a
+        bool log = collectToIORank_.isIORank();
-     * std::logic_error exception
+        eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), substep, log);
-     */
-    void attachScalarVertexData(ScalarBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
-    {
-        OPM_THROW(std::logic_error,
-                  "The EclWriter can only write element based quantities!");
-    }
 
-    /*
+        ElementContext elemCtx(simulator_);
-     * \brief Add a vertex-centered vector field to the output.
+        ElementIterator elemIt = gridView.template begin</*codim=*/0>();
-     *
+        const ElementIterator& elemEndIt = gridView.template end</*codim=*/0>();
-     * For the EclWriter, this method is a no-op which throws a
+        for (; elemIt != elemEndIt; ++elemIt) {
-     * std::logic_error exception
+            const Element& elem = *elemIt;
-     */
+            elemCtx.updatePrimaryStencil(elem);
-    void attachVectorVertexData(VectorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
+            elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
-    {
+            eclOutputModule_.processElement(elemCtx);
-        OPM_THROW(std::logic_error,
-                  "The EclWriter can only write element based quantities!");
-    }
-
-    /*
-     * \brief Add a vertex-centered tensor field to the output.
-     */
-    void attachTensorVertexData(TensorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
-    {
-        OPM_THROW(std::logic_error,
-                  "The EclWriter can only write element based quantities!");
-    }
-
-    /*!
-     * \brief Add a scalar quantity to the output.
-     *
-     * The buffer must exist at least until the call to endWrite()
-     * finishes. Modifying the buffer between the call to this method
-     * and endWrite() results in _undefined behavior_.
-     */
-    void attachScalarElementData(ScalarBuffer& buf, std::string name)
-    {
-        attachedBuffers_.push_back(std::pair<std::string, ScalarBuffer*>(name, &buf));
-    }
-
-    /*
-     * \brief Add a element-centered vector field to the output.
-     *
-     * For the EclWriter, this method is a no-op which throws a
-     * std::logic_error exception
-     */
-    void attachVectorElementData(VectorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
-    {
-        OPM_THROW(std::logic_error,
-                  "Currently, the EclWriter can only write scalar quantities!");
-    }
-
-    /*
-     * \brief Add a element-centered tensor field to the output.
-     */
-    void attachTensorElementData(TensorBuffer& buf OPM_UNUSED, std::string name OPM_UNUSED)
-    {
-        OPM_THROW(std::logic_error,
-                  "Currently, the EclWriter can only write scalar quantities!");
-    }
-
-    /*!
-     * \brief Finalizes the current writer.
-     *
-     * This means that everything will be written to disk, except if
-     * the onlyDiscard argument is true. In this case, no output is
-     * written but the 'janitorial' jobs at the end of a time step are
-     * still done.
-     */
-    void endWrite(bool onlyDiscard = false)
-    {
-        if (onlyDiscard || !enableEclOutput_() || !simulator_.episodeWillBeOver()) {
-            // detach all buffers
-            attachedBuffers_.clear();
-            return;
         }
+        eclOutputModule_.outputErrorLog();
 
-#if !HAVE_ERT
+        // collect all data to I/O rank and assign to sol
-        OPM_THROW(std::runtime_error,
+        Opm::data::Solution localCellData = fip;
-                  "The ERT libraries must be available to write ECL output!");
+        eclOutputModule_.assignToSolution(localCellData);
-#else
+        if (collectToIORank_.isParallel())
-
+            collectToIORank_.collect(localCellData);
-        // collect all data to I/O rank and store in attachedBuffers_
-        // this also reorders the data such that it fits the underlying eclGrid
-        collectToIORank_.collect( attachedBuffers_ );
 
         // write output on I/O rank
         if (collectToIORank_.isIORank()) {
-            ErtRestartFile restartFile(simulator_, reportStepIdx_);
-            restartFile.writeHeader(simulator_, reportStepIdx_);
 
-            ErtSolution solution(restartFile);
+            std::map<std::string, std::vector<double>> extraRestartData;
-            auto bufIt = attachedBuffers_.begin();
+            std::map<std::string, double> miscSummaryData;
-            const auto& bufEndIt = attachedBuffers_.end();
-            for (; bufIt != bufEndIt; ++ bufIt) {
-                const std::string& name = bufIt->first;
-                const ScalarBuffer& buffer = *bufIt->second;
 
-                std::shared_ptr<const ErtKeyword<float>>
+            // Add suggested next timestep to extra data.
-                    bufKeyword(new ErtKeyword<float>(name, buffer));
+            extraRestartData["OPMEXTRA"] = std::vector<double>(1, nextstep);
-                solution.add(bufKeyword);
+
+            // Add TCPU if simulatorReport is not defaulted.
+            if (totalSolverTime != 0.0) {
+                miscSummaryData["TCPU"] = totalSolverTime;
             }
+
+            const Opm::data::Solution& cellData = collectToIORank_.isParallel() ? collectToIORank_.globalCellData() : localCellData;
+            eclIO_->writeTimeStep(episodeIdx,
+                                  substep,
+                                  t,
+                                  cellData,
+                                  dw,
+                                  miscSummaryData,
+                                  extraRestartData,
+                                  false);
         }
 
-        // detach all buffers
-        attachedBuffers_.clear();
-
-        // next time we take the next report step
-        ++ reportStepIdx_;
 #endif
     }
 
-    /*!
+    void restartBegin() {
-     * \brief Write the multi-writer's state to a restart file.
+        std::map<std::string, Opm::RestartKey> solution_keys {{"PRESSURE" , Opm::RestartKey(Opm::UnitSystem::measure::pressure)},
-     */
+                                                         {"SWAT" , Opm::RestartKey(Opm::UnitSystem::measure::identity)},
-    template <class Restarter>
+                                                         {"SGAS" , Opm::RestartKey(Opm::UnitSystem::measure::identity)},
-    void serialize(Restarter& res)
+                                                         {"TEMP" , Opm::RestartKey(Opm::UnitSystem::measure::temperature)},
-    {
+                                                         {"RS" , Opm::RestartKey(Opm::UnitSystem::measure::gas_oil_ratio)},
-        res.serializeSectionBegin("EclWriter");
+                                                         {"RV" , Opm::RestartKey(Opm::UnitSystem::measure::oil_gas_ratio)},
-        res.serializeStream() << reportStepIdx_ << "\n";
+                                                         {"SOMAX", {Opm::UnitSystem::measure::identity, false}},
-        res.serializeSectionEnd();
+                                                         {"PCSWM_OW", {Opm::UnitSystem::measure::identity, false}},
+                                                         {"KRNSW_OW", {Opm::UnitSystem::measure::identity, false}},
+                                                         {"PCSWM_GO", {Opm::UnitSystem::measure::identity, false}},
+                                                         {"KRNSW_GO", {Opm::UnitSystem::measure::identity, false}}};
+
+        std::map<std::string, bool> extra_keys {
+            {"OPMEXTRA" , false}
+        };
+
+        unsigned episodeIdx = simulator_.episodeIndex();
+        const auto& gridView = simulator_.gridManager().gridView();
+        unsigned numElements = gridView.size(/*codim=*/0);
+        eclOutputModule_.allocBuffers(numElements, episodeIdx, simulator_.gridManager().eclState().getRestartConfig(), true, false);
+
+        auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
+        for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
+            unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
+            eclOutputModule_.setRestart(restart_values.solution, elemIdx, globalIdx);
+        }
     }
 
-    /*!
+
-     * \brief Read the multi-writer's state from a restart file.
+    const EclOutputBlackOilModule<TypeTag>& eclOutputModule() const {
-     */
+        return eclOutputModule_;
-    template <class Restarter>
-    void deserialize(Restarter& res)
-    {
-        res.deserializeSectionBegin("EclWriter");
-        res.deserializeStream() >> reportStepIdx_;
-        res.deserializeSectionEnd();
     }
 
+
 private:
     static bool enableEclOutput_()
     { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
 
-    // make sure the field is well defined if running under valgrind
-    // and make sure that all values can be displayed by paraview
-    void sanitizeBuffer_(std::vector<float>& b)
-    {
-        static bool warningPrinted = false;
-        for (size_t i = 0; i < b.size(); ++i) {
-            Opm::Valgrind::CheckDefined(b[i]);
-
-            if (!warningPrinted && !std::isfinite(b[i])) {
-                std::cerr << "WARNING: data field written to disk contains non-finite entries!\n";
-                warningPrinted = true;
-            }
-
-            // set values which are too small to 0 to avoid possible
-            // problems
-            if (std::abs(b[i]) < std::numeric_limits<float>::min()) {
-                b[i] = 0.0;
-            }
-        }
-    }
-
     const Simulator& simulator_;
-    const GridView gridView_;
+    EclOutputBlackOilModule<TypeTag> eclOutputModule_;
-
-    ElementMapper elementMapper_;
-    VertexMapper vertexMapper_;
     CollectDataToIORankType collectToIORank_;
+    std::unique_ptr<Opm::EclipseIO> eclIO_;
 
-    double curTime_;
-    unsigned reportStepIdx_;
-
-    std::list<std::pair<std::string, ScalarBuffer*> > attachedBuffers_;
 };
 } // namespace Ewoms