Use data::Solution to gather cell data including user requested data.

Previously only the cell data registered with SimulationDataContainer
war gathered during parallel output. User requested data was neglected
and a warning was issued. With this commit we intialize the local view
of data::Solution on all processes with the data registered in
SimulationDataContainer and add cell data requested by the user. This is
then gathered on the IO process, and used for the output layer. To
rudimentarily support matlab we also create a global view of
SimulationDataContainer for it.

opm/autodiff/ParallelDebugOutput.hpp

@@ -48,13 +48,19 @@ namespace Opm
         virtual ~ParallelDebugOutputInterface() {}
 
         //! \brief gather solution to rank 0 for EclipseWriter
+        //! \param localReservoirState The reservoir state
         //! \param localWellState      The well state
+        //! \param localCellData       The cell data used for eclipse output
+        //!                            (needs to include the cell data of
+        //!                            localReservoirState)
         //! \param wellStateStepNumber The step number of the well state.
         virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
+                                      const data::Solution& localCellData,
                                       const int wellStateStepNumber ) = 0;
 
         virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
+        virtual const data::Solution& globalCellData() const = 0 ;
         virtual const WellState& globalWellState() const = 0 ;
         virtual bool isIORank() const = 0;
         virtual bool isParallel() const = 0;
@@ -70,25 +76,33 @@ namespace Opm
 
         const SimulationDataContainer* globalState_;
         const WellState*      wellState_;
+        const data::Solution* globalCellData_;
 
     public:
         ParallelDebugOutput ( const GridImpl& grid,
                               const EclipseState& /* eclipseState */,
                               const int,
-                              const double* )
+                              const double*,
+                              const Opm::PhaseUsage& )
             : grid_( grid ) {}
 
         // gather solution to rank 0 for EclipseWriter
         virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
+                                      const data::Solution& localCellData,
                                       const int /* wellStateStepNumber */)
         {
             globalState_ = &localReservoirState;
             wellState_   = &localWellState;
+            globalCellData_ = &localCellData;
             return true ;
         }
 
         virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; }
+        virtual const data::Solution& globalCellData() const
+        {
+            return *globalCellData_;
+        }
         virtual const WellState& globalWellState() const { return *wellState_; }
         virtual bool isIORank () const { return true; }
         virtual bool isParallel () const { return false; }
@@ -231,12 +245,16 @@ namespace Opm
         ParallelDebugOutput( const Dune::CpGrid& otherGrid,
                              const EclipseState& eclipseState,
                              const int numPhases,
-                             const double* permeability )
+                             const double* permeability,
+                             const Opm::PhaseUsage& phaseUsage)
             : grid_(),
               eclipseState_( eclipseState ),
               permeability_( permeability ),
               toIORankComm_( otherGrid.comm() ),
-              isIORank_( otherGrid.comm().rank() == ioRank )
+              globalCellData_(new data::Solution),
+              isIORank_( otherGrid.comm().rank() == ioRank ),
+              phaseUsage_(phaseUsage)
+
         {
             const CollectiveCommunication& comm = otherGrid.comm();
             if( comm.size() > 1 )
@@ -321,6 +339,8 @@ namespace Opm
         {
             const SimulationDataContainer& localState_;
             SimulationDataContainer& globalState_;
+            const data::Solution& localCellData_;
+            data::Solution& globalCellData_;
             const WellState& localWellState_;
             WellState& globalWellState_;
             const IndexMapType& localIndexMap_;
@@ -328,14 +348,18 @@ namespace Opm
 
         public:
             PackUnPackSimulationDataContainer( const SimulationDataContainer& localState,
-                                      SimulationDataContainer& globalState,
+                                               SimulationDataContainer& globalState,
-                                      const WellState& localWellState,
+                                               const data::Solution& localCellData,
-                                      WellState& globalWellState,
+                                               data::Solution& globalCellData,
-                                      const IndexMapType& localIndexMap,
+                                               const WellState& localWellState,
-                                      const IndexMapStorageType& indexMaps,
+                                               WellState& globalWellState,
-                                      const bool isIORank )
+                                               const IndexMapType& localIndexMap,
+                                               const IndexMapStorageType& indexMaps,
+                                               const bool isIORank )
             : localState_( localState ),
               globalState_( globalState ),
+              localCellData_( localCellData ),
+              globalCellData_( globalCellData ),
               localWellState_( localWellState ),
               globalWellState_( globalWellState ),
               localIndexMap_( localIndexMap ),
@@ -351,6 +375,16 @@ namespace Opm
                         }
                     }
 
+                    // add missing data to global cell data
+                    for (const auto& pair : localCellData_) {
+                        const std::string& key = pair.first;
+                        std::size_t container_size = globalState_.numCells() *
+                            pair.second.data.size() / localState_.numCells();
+                        globalCellData_.insert(key, pair.second.dim,
+                                                std::vector<double>(container_size),
+                                                pair.second.target);
+                    }
+
                     MessageBufferType buffer;
                     pack( 0, buffer );
                     // the last index map is the local one
@@ -367,10 +401,10 @@ namespace Opm
                 }
 
                 // write all cell data registered in local state
-		for (const auto& pair : localState_.cellData()) {
+                for (const auto& pair : localCellData_) {
 		    const std::string& key = pair.first;
-                    const auto& data = pair.second;
+                    const auto& data = pair.second.data;
-                    const size_t stride = localState_.numCellDataComponents( key );
+                    const size_t stride = data.size()/localState_.numCells();
 
                     for( size_t i=0; i<stride; ++i )
                     {
@@ -388,10 +422,11 @@ namespace Opm
                 // write all cell data registered in local state
                 // we loop over the data of the local state as
                 // its order governs the order the data got received.
-                for (auto& pair : localState_.cellData()) {
+                for (auto& pair : localCellData_) {
                     const std::string& key = pair.first;
-                    auto& data = globalState_.getCellData(key);
+
-                    const size_t stride = globalState_.numCellDataComponents( key );
+                    auto& data = globalCellData_.data(key);
+                    const size_t stride = data.size() / globalState_.numCells();
 
                     for( size_t i=0; i<stride; ++i )
                     {
@@ -525,6 +560,7 @@ namespace Opm
         // gather solution to rank 0 for EclipseWriter
         bool collectToIORank( const SimulationDataContainer& localReservoirState,
                               const WellState& localWellState,
+                              const data::Solution& localCellData,
                               const int wellStateStepNumber )
         {
             if( isIORank() )
@@ -559,19 +595,32 @@ namespace Opm
             }
 
             PackUnPackSimulationDataContainer packUnpack( localReservoirState, *globalReservoirState_,
-                                                 localWellState, globalWellState_,
+                                                          localCellData, *globalCellData_,
-                                                 localIndexMap_, indexMaps_, isIORank() );
+                                                          localWellState, globalWellState_,
+                                                          localIndexMap_, indexMaps_,
+                                                          isIORank() );
 
             //toIORankComm_.exchangeCached( packUnpack );
             toIORankComm_.exchange( packUnpack );
 #ifndef NDEBUG
-            // mkae sure every process is on the same page
+            // make sure every process is on the same page
             toIORankComm_.barrier();
 #endif
+            if( isIORank() )
+            {
+                // Update values in the globalReservoirState
+                solutionToSim(*globalCellData_, phaseUsage_, *globalReservoirState_);
+            }
             return isIORank();
         }
 
         const SimulationDataContainer& globalReservoirState() const { return *globalReservoirState_; }
+
+        const data::Solution& globalCellData() const
+        {
+            return *globalCellData_;
+        }
+
         const WellState& globalWellState() const { return globalWellState_; }
 
         bool isIORank() const
@@ -600,10 +649,13 @@ namespace Opm
         IndexMapType                              localIndexMap_;
         IndexMapStorageType                       indexMaps_;
         std::unique_ptr<SimulationDataContainer>  globalReservoirState_;
+        std::unique_ptr<data::Solution>           globalCellData_;
         // this needs to be revised
         WellStateFullyImplicitBlackoil            globalWellState_;
         // true if we are on I/O rank
         const bool                                isIORank_;
+        // Phase usage needed to convert solution to simulation data container
+        Opm::PhaseUsage phaseUsage_;
     };
 #endif // #if HAVE_OPM_GRID
 

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -240,7 +240,13 @@ namespace Opm
                   const WellState& localWellState,
                   bool substep)
     {
-        writeTimeStepWithCellProperties(timer, localState, localWellState, {} , substep);
+        data::Solution localCellData{};
+        if( output_ )
+        {
+            localCellData = simToSolution(localState, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
+        }
+        writeTimeStepWithCellProperties(timer, localState, localCellData ,
+                                        localWellState, substep);
     }
 
 
@@ -252,8 +258,8 @@ namespace Opm
     writeTimeStepWithCellProperties(
                   const SimulatorTimerInterface& timer,
                   const SimulationDataContainer& localState,
+                  const data::Solution& localCellData,
                   const WellState& localWellState,
-                  const data::Solution& cellData,
                   bool substep)
     {
         // VTK output (is parallel if grid is parallel)
@@ -272,9 +278,12 @@ namespace Opm
             int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
                 (timer.reportStepNum() - 1) : timer.reportStepNum();
             // collect all solutions to I/O rank
-            isIORank = parallelOutput_->collectToIORank( localState, localWellState, wellStateStepNumber );
+            isIORank = parallelOutput_->collectToIORank( localState, localWellState,
+                                                         localCellData,
+                                                         wellStateStepNumber );
         }
 
+        const data::Solution& cellData = ( parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalCellData() : localCellData;
         const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
         const WellState& wellState  = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
 
@@ -314,12 +323,11 @@ namespace Opm
             if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
                 std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
             } else {
-                data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...)
+                // ... insert "extra" data (KR, VISC, ...)
-                combined_sol.insert(simProps.begin(), simProps.end());           // ... insert "extra" data (KR, VISC, ...)
                 eclWriter_->writeTimeStep(timer.reportStepNum(),
                                           substep,
                                           timer.simulationTimeElapsed(),
-                                          combined_sol,
+                                          simProps,
                                           wellState.report(phaseUsage_));
             }
         }

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -237,8 +237,8 @@ namespace Opm
         void writeTimeStepWithCellProperties(
                            const SimulatorTimerInterface& timer,
                            const SimulationDataContainer& reservoirState,
+                           const data::Solution& cellData,
                            const Opm::WellState& wellState,
-                           const data::Solution& solution,
                            bool substep = false);
 
         /*!
@@ -260,7 +260,7 @@ namespace Opm
         void writeTimeStepSerial(const SimulatorTimerInterface& timer,
                                  const SimulationDataContainer& reservoirState,
                                  const Opm::WellState& wellState,
-                                 const data::Solution& cellData,
+                                 const data::Solution& simProps,
                                  bool substep);
 
         /** \brief return output directory */
@@ -327,7 +327,7 @@ namespace Opm
                          const Opm::PhaseUsage &phaseUsage,
                          const double* permeability )
       : output_( param.getDefault("output", true) ),
-        parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability ) : 0 ),
+        parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, permeability, phaseUsage ) : 0 ),
         outputDir_( output_ ? param.getDefault("output_dir", std::string("output")) : "." ),
         output_interval_( output_ ? param.getDefault("output_interval", 1): 0 ),
         lastBackupReportStep_( -1 ),
@@ -787,57 +787,21 @@ namespace Opm
                   const Model& physicalModel,
                   bool substep)
     {
-        data::Solution cellData{};
+        data::Solution localCellData{};
         const RestartConfig& restartConfig = eclipseState_.getRestartConfig();
         const SummaryConfig& summaryConfig = eclipseState_.getSummaryConfig();
         const int reportStepNum = timer.reportStepNum();
         bool logMessages = output_ && parallelOutput_->isIORank();
         
-        if( output_ && !parallelOutput_->isParallel() )
+        if( output_ )
         {
-
+            localCellData = simToSolution(localState, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
-            detail::getRestartData( cellData, phaseUsage_, physicalModel,
+            detail::getRestartData( localCellData, phaseUsage_, physicalModel,
                                     restartConfig, reportStepNum, logMessages );
-            detail::getSummaryData( cellData, phaseUsage_, physicalModel, summaryConfig );
+            detail::getSummaryData( localCellData, phaseUsage_, physicalModel, summaryConfig );
         }
-        else
-        {
-            if ( logMessages )
-            {
-                std::map<std::string, int> rstKeywords = restartConfig.getRestartKeywords(reportStepNum);
-                std::vector<const char*> keywords = 
-                    { "WIP", "OIPL", "OIPG", "OIP", "GIPG", "GIPL", "GIP",
-                      "RPV", "FRPH", "RPRH"};
-                
-                std::ostringstream str;
-                str << "Output of restart/summary config not supported in parallel. Requested keywords were ";
-                std::size_t no_kw = 0;
-                    
-                auto func = [&] (const char* kw)
-                    {
-                        if ( detail::hasFRBKeyword(summaryConfig, kw) )
-                        {
-                            str << kw << " ";
-                            ++ no_kw;
-                        }
-                    };
-
-                std::for_each(keywords.begin(), keywords.end(), func);
-
-                for (auto& keyValue : rstKeywords)
-                {
-                        str << keyValue.first << " ";
-                        ++ no_kw;
-                }
-
-                if ( no_kw )
-                {
-                    Opm::OpmLog::warning("Unhandled ouput request", str.str());
-                }
-            }
-        }        
         
-        writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep);
+        writeTimeStepWithCellProperties(timer, localState, localCellData, localWellState, substep);
     }
 }
 #endif