Merge pull request #286 from dr-robertk/PR/add-support-for-restart

BlackoilOutputWriter: added support for backup and restore.

CMakeLists_files.cmake

@@ -94,6 +94,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/autodiff/AutoDiffBlock.hpp
 	opm/autodiff/AutoDiffHelpers.hpp
 	opm/autodiff/AutoDiff.hpp
+	opm/autodiff/BackupRestore.hpp
 	opm/autodiff/BlackoilPropsAd.hpp
 	opm/autodiff/BlackoilPropsAdFromDeck.hpp
 	opm/autodiff/BlackoilPropsAdInterface.hpp

opm/autodiff/BackupRestore.hpp

@@ -0,0 +1,342 @@
+/*
+  Copyright 2014 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_BACKUPRESTORE_HEADER_INCLUDED
+#define OPM_BACKUPRESTORE_HEADER_INCLUDED
+
+#include <iostream>
+
+#include <opm/core/simulator/SimulatorState.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+namespace Opm {
+
+    namespace {
+        template <class T>
+        void writeValue( std::ostream& stream, const T& value )
+        {
+            stream.write( (const char *) &value, sizeof( T ) );
+        }
+
+        template <class T>
+        void readValue( std::istream& stream, T& value )
+        {
+            stream.read( (char *) &value, sizeof( T ) );
+        }
+
+        // write any container that provides a size method and a data method
+        template <class Vector>
+        void writeContainer( std::ostream& stream, const Vector& vector)
+        {
+            typedef typename Vector :: value_type T;
+            unsigned int size = vector.size() * sizeof( T );
+            writeValue( stream, size );
+            if( size > 0 ) {
+                stream.write( (const char *) vector.data(), size );
+            }
+        }
+
+        // write map where the key is a std::string
+        // and the value is a std::vector or std::array
+        template <class Map>
+        void writeMap( std::ostream& stream, const Map& m)
+        {
+            const unsigned int mapEntries = m.size();
+            writeValue( stream, mapEntries );
+            if( mapEntries > 0 )
+            {
+                const auto& end = m.end();
+                for(auto it = m.begin(); it != end; ++it )
+                {
+                    // write key (assume that key is a container)
+                    writeContainer( stream, (*it).first );
+                    // write value (assume that value is a container)
+                    writeContainer( stream, (*it).second );
+                }
+            }
+        }
+
+        template <class Container>
+        void resizeContainer( Container& container, size_t size )
+        {
+            container.resize( size );
+        }
+
+        template <class T, unsigned long n>
+        void resizeContainer( std::array<T, n>& a, size_t size )
+        {
+            assert( int(size) == int(n) );
+        }
+
+        template <class Container>
+        void readContainerImpl( std::istream& stream, Container& container, const bool adjustSize )
+        {
+            typedef typename Container :: value_type T;
+            unsigned int dataSize = 0;
+            readValue( stream, dataSize );
+            if( adjustSize && dataSize > 0 ) {
+                resizeContainer( container, dataSize/sizeof(T) );
+            }
+
+            if( dataSize != container.size() * sizeof( T ) )
+            {
+                OPM_THROW(std::logic_error,
+                        "Size of stored data and simulation data does not match"
+                        << dataSize << " " << (container.size() * sizeof( T )) );
+            }
+            if( dataSize > 0 ) {
+                stream.read( (char *) container.data(), dataSize );
+            }
+        }
+
+        template <class Container>
+        void readAndResizeContainer( std::istream& stream, Container& container )
+        {
+            readContainerImpl( stream, container, true );
+        }
+
+        template <class Container>
+        void readContainer( std::istream& stream, Container& container )
+        {
+            readContainerImpl( stream, container, false );
+        }
+
+        template <class Map>
+        void readMap( std::istream& stream, Map& m)
+        {
+            m.clear();
+            unsigned int mapEntries = 0;
+            readValue( stream, mapEntries );
+            for( unsigned int entry = 0; entry<mapEntries; ++entry )
+            {
+                std::pair< typename Map :: key_type, typename Map :: mapped_type > mapEntry;
+                // read key
+                readAndResizeContainer( stream, mapEntry.first );
+                // read values
+                readContainer( stream, mapEntry.second );
+
+                // insert entry into map
+                m.insert( mapEntry );
+            }
+        }
+
+        enum { SimulatorStateId = 0,
+               WellStateId = 1,
+               BlackoilStateId = 2,
+               WellStateFullyImplicitBackoilId = 3 };
+
+        inline int objectId( const SimulatorState& state ) {
+            return SimulatorStateId;
+        }
+        inline int objectId( const WellState& state ) {
+            return WellStateId;
+        }
+        inline int objectId( const BlackoilState& state ) {
+            return BlackoilStateId;
+        }
+        inline int objectId( const WellStateFullyImplicitBlackoil& state ) {
+            return WellStateFullyImplicitBackoilId;
+        }
+
+        template <class State>
+        void checkObjectId( std::istream& in, const State& state )
+        {
+            // read id and compare with object
+            int id = -1;
+            readValue( in, id );
+            if( id != objectId( state ) ) {
+                OPM_THROW(std::logic_error,"backup-restore object type mismatch");
+            }
+        }
+    }
+
+    // SimulatorState
+    inline
+    std::ostream& operator << (std::ostream& out, const SimulatorState& state )
+    {
+        // write id of object to stream
+        writeValue( out, objectId( state ) );
+
+        const int numPhases = state.numPhases();
+        writeValue( out, numPhases );
+
+        // write variables
+        writeContainer( out, state.pressure() );
+        writeContainer( out, state.temperature() );
+        writeContainer( out, state.facepressure() );
+        writeContainer( out, state.faceflux() );
+        writeContainer( out, state.saturation() );
+
+        return out;
+    }
+
+    inline
+    std::istream& operator >> (std::istream& in, SimulatorState& state )
+    {
+        // check id of stored object
+        checkObjectId( in, state );
+
+        int numPhases = 0;
+        readValue( in, numPhases );
+
+        if( numPhases != state.numPhases() )
+            OPM_THROW(std::logic_error,"num phases wrong");
+
+        // read variables
+        readContainer( in, state.pressure() );
+        readContainer( in, state.temperature() );
+        readContainer( in, state.facepressure() );
+        readContainer( in, state.faceflux() );
+        readContainer( in, state.saturation() );
+
+        return in;
+    }
+
+    // BlackoilState
+    inline
+    std::ostream& operator << (std::ostream& out, const BlackoilState& state )
+    {
+        // write id of object to stream
+        writeValue( out, objectId( state ) );
+
+        // backup simulator state
+        const SimulatorState& simstate = static_cast< const SimulatorState& > (state);
+        out << simstate;
+
+        // backup additional blackoil state variables
+        writeContainer( out, state.surfacevol() );
+        writeContainer( out, state.gasoilratio() );
+        writeContainer( out, state.rv() );
+
+        return out;
+    }
+
+    inline
+    std::istream& operator >> (std::istream& in, BlackoilState& state )
+    {
+        // check id of stored object
+        checkObjectId( in, state );
+
+        // restore simulator state
+        SimulatorState& simstate = static_cast< SimulatorState& > (state);
+        in >> simstate;
+
+        // restore additional blackoil state variables
+        readContainer( in, state.surfacevol() );
+        readContainer( in, state.gasoilratio() );
+        readContainer( in, state.rv() );
+
+        return in;
+    }
+
+    // WellState
+    inline
+    std::ostream& operator << (std::ostream& out, const WellState& state )
+    {
+        // write id of object to stream
+        writeValue( out, objectId( state ) );
+
+        // backup well state
+        writeContainer( out, state.bhp() );
+        writeContainer( out, state.temperature() );
+        writeContainer( out, state.wellRates() );
+        writeContainer( out, state.perfRates() );
+        writeContainer( out, state.perfPress() );
+
+        return out;
+    }
+
+    inline
+    std::istream& operator >> (std::istream& in, WellState& state )
+    {
+        // check id of stored object
+        checkObjectId( in, state );
+
+        // restore well state
+        readAndResizeContainer( in, state.bhp() );
+        readAndResizeContainer( in, state.temperature() );
+        readAndResizeContainer( in, state.wellRates() );
+        readAndResizeContainer( in, state.perfRates() );
+        readAndResizeContainer( in, state.perfPress() );
+
+        return in;
+    }
+
+    // WellStateFullyImplicitBlackoil
+    inline
+    std::ostream& operator << (std::ostream& out, const WellStateFullyImplicitBlackoil& state )
+    {
+        // write id of object to stream
+        writeValue( out, objectId( state ) );
+
+        // backup well state
+        const WellState& wellState = static_cast< const WellState& > (state);
+        out << wellState;
+
+        const int numWells = state.numWells();
+        writeValue( out, numWells );
+        if( numWells > 0 )
+        {
+            const int numPhases = state.numPhases();
+            writeValue( out, numPhases );
+
+            // backup additional variables
+            writeContainer( out, state.perfPhaseRates() );
+            writeContainer( out, state.currentControls() );
+            writeMap( out, state.wellMap() );
+        }
+
+        return out;
+    }
+
+    inline
+    std::istream& operator >> (std::istream& in, WellStateFullyImplicitBlackoil& state )
+    {
+        // check id of stored object
+        checkObjectId( in, state );
+
+        // restore well state
+        WellState& wellState = static_cast< WellState& > (state);
+        in >> wellState;
+
+        int numWells = 0;
+        readValue( in, numWells );
+        if( numWells != state.numWells() )
+            OPM_THROW(std::logic_error,"wrong numWells");
+
+        if( numWells > 0 )
+        {
+            int numPhases = 0;
+            readValue( in, numPhases );
+            if( numPhases != state.numPhases() )
+                OPM_THROW(std::logic_error,"wrong numPhases");
+
+            // restore additional variables
+            readAndResizeContainer( in, state.perfPhaseRates()  );
+            readAndResizeContainer( in, state.currentControls() );
+            readMap( in, state.wellMap() );
+        }
+
+        return in;
+    }
+
+} // namespace Opm
+
+#endif // OPM_BACKUPRESTORE_HEADER_INCLUDED

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -25,8 +25,10 @@
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/Units.hpp>
 
 #include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/BackupRestore.hpp>
 
 #include <sstream>
 #include <iomanip>
@@ -248,5 +250,100 @@ namespace Opm
         if ( eclWriter_ ) {
             eclWriter_->writeTimeStep(timer, state, wellState);
         }
+        // write backup file
+        if( backupfile_ )
+        {
+            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) );
+
+                const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state);
+                backupfile_ << boState;
+
+                const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState);
+                backupfile_ << boWellState;
+                /*
+                const WellStateFullyImplicitBlackoil* boWellState =
+                    dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState);
+                if( boWellState ) {
+                    backupfile_ << (*boWellState);
+                }
+                else
+                    OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed");
+                */
+                backupfile_ << std::flush;
+            }
+        }
+    }
+
+    void
+    BlackoilOutputWriter::
+    restore(SimulatorTimerInterface& timer,
+            BlackoilState& state,
+            WellStateFullyImplicitBlackoil& wellState,
+            const std::string& filename,
+            const int desiredResportStep )
+    {
+        std::ifstream restorefile( filename );
+        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;
+
+                writeTimeStep( 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;
+                }
+
+                // next step
+                timer.advance();
+
+                // read next report step
+                restorefile.read( (char *) &reportStep, sizeof(int) );
+                if( timer.reportStepNum() != reportStep ) {
+                    break;
+                }
+            }
+        }
+        else
+        {
+            std::cerr << "Warning: Couldn't open restore file '" << filename << "'" << std::endl;
+        }
     }
 }

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -32,6 +32,8 @@
 
 #include <opm/autodiff/GridHelpers.hpp>
 
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
 #include <string>
 #include <sstream>
 #include <iomanip>
@@ -112,7 +114,7 @@ namespace Opm
         }
     }
 
-    /** Wrapper class satisfying the OutputWriter interface */
+    /** Wrapper class satisfying the OutputWriter interface and writing VTK output. */
     template <class Grid>
     class BlackoilVTKWriter : public OutputWriter
     {
@@ -140,6 +142,7 @@ namespace Opm
         const std::string outputDir_;
     };
 
+    /** Wrapper class satisfying the OutputWriter interface and writing Matlab output. */
     template <class Grid>
     class BlackoilMatlabWriter : public OutputWriter
     {
@@ -169,12 +172,12 @@ namespace Opm
             }
             outputWellStateMatlab(wellState, timer.currentStepNum(), outputDir_);
         }
-
     protected:
         const Grid& grid_;
         const std::string outputDir_;
     };
 
+    /** \brief Wrapper class for VTK, Matlab, and ECL output. */
     class BlackoilOutputWriter : public OutputWriter
     {
     public:
@@ -199,12 +202,21 @@ namespace Opm
         /** \brief return true if output is enabled */
         const bool output () const { return output_; }
 
+        void restore(SimulatorTimerInterface& timer,
+                     BlackoilState& state,
+                     WellStateFullyImplicitBlackoil& wellState,
+                     const std::string& filename,
+                     const int desiredReportStep);
+
     protected:
         // Parameters for output.
         const bool output_;
         const std::string outputDir_;
         const int output_interval_;
 
+        int lastBackupReportStep_;
+
+        std::ofstream backupfile_;
         std::unique_ptr< OutputWriter  > vtkWriter_;
         std::unique_ptr< OutputWriter  > matlabWriter_;
         std::unique_ptr< EclipseWriter > eclWriter_;
@@ -226,6 +238,7 @@ namespace Opm
       : output_( param.getDefault("output", true) ),
         outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ),
         output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ),
+        lastBackupReportStep_( -1 ),
         vtkWriter_( output_ && param.getDefault("output_vtk",false) ?
                      new BlackoilVTKWriter< Grid >( grid, outputDir_ ) : 0 ),
         matlabWriter_( output_ && param.getDefault("output_matlab", false) ?
@@ -246,6 +259,12 @@ namespace Opm
             catch (...) {
                 OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
             }
+
+            std::string backupfilename = param.getDefault("backupfile", std::string("") );
+            if( ! backupfilename.empty() )
+            {
+                backupfile_.open( backupfilename );
+            }
         }
     }
 }

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -227,6 +227,14 @@ namespace Opm
         // init output writer
         output_writer_.writeInit( timer );
 
+        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 );
+        }
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -186,6 +186,7 @@ namespace Opm
         }
 
         const WellMapType& wellMap() const { return wellMap_; }
+        WellMapType& wellMap() { return wellMap_; }
 
     private:
         std::vector<double> perfphaserates_;