Merge pull request #772 from andlaus/use_deck_units_for_summary_output

Use deck units for summary output

opm/core/io/OutputWriter.cpp

@@ -47,10 +47,12 @@ private:
 /// Psuedo-constructor, can appear in template
 template <typename Format> unique_ptr <OutputWriter>
 create (const ParameterGroup& params,
+        std::shared_ptr <const Deck> deck,
         std::shared_ptr <const EclipseState> eclipseState,
         const Opm::PhaseUsage &phaseUsage,
         std::shared_ptr <const UnstructuredGrid> grid) {
     return unique_ptr <OutputWriter> (new Format (params,
+                                                  deck,
                                                   eclipseState,
                                                   phaseUsage,
                                                   grid->number_of_cells,
@@ -65,6 +67,7 @@ create (const ParameterGroup& params,
 /// to the list below!
 typedef map <const char*, unique_ptr <OutputWriter> (*)(
         const ParameterGroup&,
+        std::shared_ptr <const Deck> deck,
         std::shared_ptr <const EclipseState> eclipseState,
         const Opm::PhaseUsage &phaseUsage,
         std::shared_ptr <const UnstructuredGrid>)> map_t;
@@ -76,6 +79,7 @@ map_t FORMATS = {
 
 unique_ptr <OutputWriter>
 OutputWriter::create (const ParameterGroup& params,
+                      std::shared_ptr <const Deck> deck,
                       std::shared_ptr <const EclipseState> eclipseState,
                       const Opm::PhaseUsage &phaseUsage,
                       std::shared_ptr <const UnstructuredGrid> grid) {
@@ -93,7 +97,7 @@ OutputWriter::create (const ParameterGroup& params,
         // invoke the constructor for the type if we found the keyword
         // and put the pointer to this writer onto the list
         if (params.getDefault <bool> (name, false)) {
-            list->push_front (it->second (params, eclipseState, phaseUsage, grid));
+            list->push_front (it->second (params, deck, eclipseState, phaseUsage, grid));
         }
     }
 

opm/core/io/OutputWriter.hpp

@@ -21,6 +21,7 @@
 #define OPM_OUTPUT_WRITER_HPP
 
 #include <memory>  // unique_ptr, shared_ptr
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
 
 struct UnstructuredGrid;
@@ -107,6 +108,7 @@ public:
      */
     static std::unique_ptr <OutputWriter>
     create (const parameter::ParameterGroup& params,
+            std::shared_ptr <const Deck> deck,
             std::shared_ptr <const EclipseState> eclipseState,
             const Opm::PhaseUsage &phaseUsage,
             std::shared_ptr <const UnstructuredGrid> grid);

opm/core/io/eclipse/EclipseWriter.cpp

@@ -1,5 +1,5 @@
 /*
-  Copyright (c) 2013-2014 Andreas Lauser
+  Copyright (c) 2013-2015 Andreas Lauser
   Copyright (c) 2013 SINTEF ICT, Applied Mathematics.
   Copyright (c) 2013 Uni Research AS
   Copyright (c) 2015 IRIS AS
@@ -95,10 +95,10 @@ void restrictAndReorderToActiveCells(std::vector<double> &data,
 }
 
 // convert the units of an array
-void convertFromSiTo(std::vector<double> &siValues, double toSiConversionFactor)
+void convertFromSiTo(std::vector<double> &siValues, double toSiConversionFactor, double toSiOffset = 0)
 {
     for (size_t curIdx = 0; curIdx < siValues.size(); ++curIdx) {
-        siValues[curIdx] = unit::convert::to(siValues[curIdx], toSiConversionFactor);
+        siValues[curIdx] = unit::convert::to(siValues[curIdx] - toSiOffset, toSiConversionFactor);
     }
 }
 
@@ -417,11 +417,17 @@ public:
         boost::filesystem::path casePath(outputDir);
         casePath /= boost::to_upper_copy(baseName);
 
+        // convert the start time to seconds since 1970-1-1@00:00:00
+        boost::posix_time::ptime startTime
+            = timer.startDateTime();
+        tm t = boost::posix_time::to_tm(startTime);
+        double secondsSinceEpochStart = std::mktime(&t);
+
         ertHandle_ = ecl_sum_alloc_writer(casePath.string().c_str(),
                                           false, /* formatted   */
                                           true,  /* unified     */
                                           ":",    /* join string */
-                                          timer.simulationTimeElapsed(),
+                                          secondsSinceEpochStart,
                                           nx,
                                           ny,
                                           nz);
@@ -443,7 +449,8 @@ public:
     // on the classes defined in the following.
 
     // add rate variables for each of the well in the input file
-    void addAllWells(Opm::EclipseStateConstPtr eclipseState,
+    void addAllWells(Opm::DeckConstPtr deck,
+                     Opm::EclipseStateConstPtr eclipseState,
                      const PhaseUsage& uses);
     void writeTimeStep(int writeStepIdx,
                        const SimulatorTimerInterface& timer,
@@ -612,6 +619,8 @@ public:
     { return ertHandle_; }
 
 protected:
+
+
     void updateTimeStepWellIndex_(const std::map<std::string, int>& nameToIdxMap)
     {
         const std::string& wellName = well_->name();
@@ -724,7 +733,8 @@ public:
              Opm::WellConstPtr well,
              PhaseUsage uses,
              BlackoilPhases::PhaseIndex phase,
-             WellType type)
+             WellType type,
+             bool useFieldUnits)
         : WellReport(summary,
                      eclipseState,
                      well,
@@ -732,8 +742,9 @@ public:
                      phase,
                      type,
                      'R',
-                     "SM3/DAY" /* surf. cub. m. per day */)
+                     handleUnit_(phase, useFieldUnits))
-    { }
+    {
+    }
 
     virtual double retrieveValue(const int /* writeStepIdx */,
                                  const SimulatorTimerInterface& timer,
@@ -755,8 +766,40 @@ public:
         // TODO: Why only positive rates?
         using namespace Opm::unit;
         return convert::to(std::max(0., rate(wellState)),
-                           cubic(meter)/day);
+                           targetRateToSiConversionFactor_);
     }
+
+private:
+    const std::string handleUnit_(BlackoilPhases::PhaseIndex phase, bool useField) {
+        using namespace Opm::unit;
+        if (phase == BlackoilPhases::Liquid || phase == BlackoilPhases::Aqua) {
+            if (useField) {
+                unitName_ = "STB/DAY";
+                targetRateToSiConversionFactor_ = stb/day; // m^3/s -> STB/day
+            }
+            else {
+                unitName_ = "SM3/DAY";
+                targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day
+            }
+        }
+        else if (phase == BlackoilPhases::Vapour) {
+            if (useField) {
+                unitName_ = "MSCF/DAY";
+                targetRateToSiConversionFactor_ = 1000*cubic(feet)/day; // m^3/s -> MSCF^3/day
+            }
+            else {
+                unitName_ = "SM3/DAY";
+                targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day
+            }
+        }
+        else
+            OPM_THROW(std::logic_error,
+                      "Unexpected phase " << phase);
+        return unitName_;
+    }
+
+    const char* unitName_;
+    double targetRateToSiConversionFactor_;
 };
 
 /// Monitors the total production in a well.
@@ -768,7 +811,8 @@ public:
               Opm::WellConstPtr well,
               PhaseUsage uses,
               BlackoilPhases::PhaseIndex phase,
-              WellType type)
+              WellType type,
+              bool useFieldUnits)
         : WellReport(summary,
                      eclipseState,
                      well,
@@ -776,7 +820,7 @@ public:
                      phase,
                      type,
                      'T',
-                     "SM3" /* surface cubic meter */ )
+                     handleUnit_(phase, useFieldUnits))
           // nothing produced when the reporting starts
         , total_(0.)
     { }
@@ -812,10 +856,41 @@ public:
         // add this timesteps production to the total
         total_ += intg;
         // report the new production total
-        return total_;
+        return unit::convert::to(total_, targetRateToSiConversionFactor_);
     }
 
 private:
+    const std::string handleUnit_(BlackoilPhases::PhaseIndex phase, bool useField) {
+        using namespace Opm::unit;
+        if (phase == BlackoilPhases::Liquid || phase == BlackoilPhases::Aqua) {
+            if (useField) {
+                unitName_ = "STB/DAY";
+                targetRateToSiConversionFactor_ = stb/day; // m^3/s -> STB/day
+            }
+            else {
+                unitName_ = "SM3/DAY";
+                targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day
+            }
+        }
+        else if (phase == BlackoilPhases::Vapour) {
+            if (useField) {
+                unitName_ = "MSCF/DAY";
+                targetRateToSiConversionFactor_ = 1000*cubic(feet)/day; // m^3/s -> MSCF^3/day
+            }
+            else {
+                unitName_ = "SM3/DAY";
+                targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day
+            }
+        }
+        else
+            OPM_THROW(std::logic_error,
+                      "Unexpected phase " << phase);
+        return unitName_;
+    }
+
+    const char* unitName_;
+    double targetRateToSiConversionFactor_;
+
     /// Aggregated value of the course of the simulation
     double total_;
 };
@@ -829,7 +904,8 @@ public:
             Opm::WellConstPtr well,
             PhaseUsage uses,
             BlackoilPhases::PhaseIndex phase,
-            WellType type)
+            WellType type,
+            bool useFieldUnits)
         : WellReport(summary,
                      eclipseState,
                      well,
@@ -837,7 +913,7 @@ public:
                      phase,
                      type,
                      'B',
-                     "Pascal")
+                     handleUnit_(phase, useFieldUnits))
     { }
 
     virtual double retrieveValue(const int /* writeStepIdx */,
@@ -856,8 +932,27 @@ public:
             return 0.0;
         }
 
-        return bhp(wellState);
+        return unit::convert::to(bhp(wellState), targetRateToSiConversionFactor_);
     }
+
+private:
+    const std::string handleUnit_(BlackoilPhases::PhaseIndex phase, bool useField) {
+        using namespace Opm::unit;
+
+        if (useField) {
+            unitName_ = "PSIA";
+            targetRateToSiConversionFactor_ = psia; // Pa -> PSI
+        }
+        else {
+            unitName_ = "BARSA";
+            targetRateToSiConversionFactor_ = barsa; // Pa -> bar
+        }
+
+        return unitName_;
+    }
+
+    const char* unitName_;
+    double targetRateToSiConversionFactor_;
 };
 
 // no inline implementation of this since it depends on the
@@ -891,10 +986,13 @@ void Summary::writeTimeStep(int writeStepIdx,
     ecl_sum_fwrite(ertHandle());
 }
 
-void Summary::addAllWells(Opm::EclipseStateConstPtr eclipseState,
+void Summary::addAllWells(Opm::DeckConstPtr deck,
+                          Opm::EclipseStateConstPtr eclipseState,
                           const PhaseUsage& uses)
 {
     eclipseState_ = eclipseState;
+    bool useFieldUnits = !deck->hasKeyword("METRIC");
+
     // TODO: Only create report variables that are requested with keywords
     // (e.g. "WOPR") in the input files, and only for those wells that are
     // mentioned in those keywords
@@ -919,7 +1017,8 @@ void Summary::addAllWells(Opm::EclipseStateConstPtr eclipseState,
                                          wells[wellIdx],
                                          uses,
                                          ertPhaseIdx,
-                                         wellType)));
+                                         wellType,
+                                         useFieldUnits)));
                 // W{O,G,W}{I,P}T
                 addWell(std::unique_ptr <WellReport>(
                             new WellTotal(*this,
@@ -927,7 +1026,8 @@ void Summary::addAllWells(Opm::EclipseStateConstPtr eclipseState,
                                           wells[wellIdx],
                                           uses,
                                           ertPhaseIdx,
-                                          wellType)));
+                                          wellType,
+                                          useFieldUnits)));
             }
         }
     }
@@ -949,7 +1049,8 @@ void Summary::addAllWells(Opm::EclipseStateConstPtr eclipseState,
                                 wells[wellIdx],
                                 uses,
                                 ertPhaseIdx,
-                                WELL_TYPES[0])));
+                                WELL_TYPES[0],
+                                useFieldUnits)));
     }
 }
 } // end namespace EclipseWriterDetails
@@ -1064,7 +1165,7 @@ void EclipseWriter::writeInit(const SimulatorTimerInterface &timer)
                                                      eclGrid->getNX(),
                                                      eclGrid->getNY(),
                                                      eclGrid->getNZ()));
-    summary_->addAllWells(eclipseState_, phaseUsage_);
+    summary_->addAllWells(deck_, eclipseState_, phaseUsage_);
 }
 
 // implementation of the writeTimeStep method
@@ -1149,6 +1250,12 @@ void EclipseWriter::writeTimeStep(const SimulatorTimerInterface& timer,
 
         sol.add(EclipseWriterDetails::Keyword<float>("PRESSURE", pressure));
 
+        // write the cell temperature
+        std::vector<double> temperature = reservoirState.temperature();
+        EclipseWriterDetails::convertFromSiTo(temperature, deckToSiTemperatureFactor_, deckToSiTemperatureOffset_);
+        EclipseWriterDetails::restrictAndReorderToActiveCells(temperature, gridToEclipseIdx_.size(), gridToEclipseIdx_.data());
+        sol.add(EclipseWriterDetails::Keyword<float>("TEMP", temperature));
+
         std::vector<double> saturation_water;
         std::vector<double> saturation_gas;
 
@@ -1225,11 +1332,13 @@ void EclipseWriter::writeTimeStep(const SimulatorTimerInterface& timer,
 
 
 EclipseWriter::EclipseWriter(const parameter::ParameterGroup& params,
+                             Opm::DeckConstPtr deck,
                              Opm::EclipseStateConstPtr eclipseState,
                              const Opm::PhaseUsage &phaseUsage,
                              int numCells,
                              const int* compressedToCartesianCellIdx)
-    : eclipseState_(eclipseState)
+    : deck_(deck)
+    , eclipseState_(eclipseState)
     , numCells_(numCells)
     , compressedToCartesianCellIdx_(compressedToCartesianCellIdx)
     , gridToEclipseIdx_(numCells, int(-1) )
@@ -1266,6 +1375,12 @@ EclipseWriter::EclipseWriter(const parameter::ParameterGroup& params,
     deckToSiPressure_ =
         eclipseState->getDeckUnitSystem()->parse("Pressure")->getSIScaling();
 
+    // factor and offset from the temperature values given in the deck to Kelvin
+    deckToSiTemperatureFactor_ =
+        eclipseState->getDeckUnitSystem()->parse("Temperature")->getSIScaling();
+    deckToSiTemperatureOffset_ =
+        eclipseState->getDeckUnitSystem()->parse("Temperature")->getSIOffset();
+
     init(params);
 }
 

opm/core/io/eclipse/EclipseWriter.hpp

@@ -65,6 +65,7 @@ public:
      *        binary files using ERT.
      */
     EclipseWriter(const parameter::ParameterGroup& params,
+                  Opm::DeckConstPtr deck,
                   Opm::EclipseStateConstPtr eclipseState,
                   const Opm::PhaseUsage &phaseUsage,
                   int numCells,
@@ -106,12 +107,15 @@ public:
     static ert_ecl_unit_enum convertUnitTypeErtEclUnitEnum(UnitSystem::UnitType unit);
 
 private:
+    Opm::DeckConstPtr deck_;
     Opm::EclipseStateConstPtr eclipseState_;
     int numCells_;
     std::array<int, 3> cartesianSize_;
     const int* compressedToCartesianCellIdx_;
     std::vector< int > gridToEclipseIdx_;
     double deckToSiPressure_;
+    double deckToSiTemperatureFactor_;
+    double deckToSiTemperatureOffset_;
     bool enableOutput_;
     int outputInterval_;
     int writeStepIdx_;

opm/core/simulator/SimulatorOutput.cpp

@@ -30,6 +30,7 @@ using namespace Opm;
 
 SimulatorOutputBase::SimulatorOutputBase (
         const parameter::ParameterGroup& params,
+        std::shared_ptr <const Deck> deck,
         std::shared_ptr <const EclipseState> eclipseState,
         const Opm::PhaseUsage &phaseUsage,
         std::shared_ptr <const UnstructuredGrid> grid,
@@ -45,7 +46,7 @@ SimulatorOutputBase::SimulatorOutputBase (
 
     // process parameters into a writer. we don't setup a new chain in
     // every timestep!
-    , writer_ (std::move (OutputWriter::create (params, eclipseState, phaseUsage, grid)))
+    , writer_ (std::move (OutputWriter::create (params, deck, eclipseState, phaseUsage, grid)))
 
     // always start from the first timestep
     , next_ (0) {

opm/core/simulator/SimulatorOutput.hpp

@@ -56,6 +56,7 @@ protected:
      * need to pick them up from the object members.
      */
     SimulatorOutputBase (const parameter::ParameterGroup& p,
+                         std::shared_ptr <const Deck> deck,
                          std::shared_ptr <const EclipseState> eclipseState,
                          const Opm::PhaseUsage &phaseUsage,
                          std::shared_ptr <const UnstructuredGrid> grid,
@@ -145,6 +146,7 @@ private:
 template <typename Simulator>
 struct SimulatorOutput : public SimulatorOutputBase {
 	SimulatorOutput (const parameter::ParameterGroup& params,
+                     std::shared_ptr <const Deck> deck,
                      std::shared_ptr <const EclipseState> eclipseState,
                      const Opm::PhaseUsage &phaseUsage,
                      std::shared_ptr <const UnstructuredGrid> grid,
@@ -153,7 +155,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
                      std::shared_ptr <const WellState> wellState,
                      std::shared_ptr <Simulator> sim)
         // send all other parameters to base class
-        : SimulatorOutputBase (params, eclipseState, phaseUsage,
+        : SimulatorOutputBase (params, deck, eclipseState, phaseUsage,
                                grid, timer, state, wellState)
 
         // store reference to simulator in derived class
@@ -169,6 +171,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
      * the arguments passed exceeds the lifetime of this object.
      */
     SimulatorOutput (const parameter::ParameterGroup& params,
+                     const Deck& deck,
                      const EclipseState& eclipseState,
                      const Opm::PhaseUsage &phaseUsage,
                      const UnstructuredGrid& grid,
@@ -178,6 +181,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
                      Simulator& sim)
         // send all other parameters to base class
         : SimulatorOutputBase (params,
+                               share_obj (deck),
                                share_obj (eclipseState),
                                phaseUsage,
                                share_obj (grid),

tests/test_EclipseWriteRFTHandler.cpp

@@ -136,6 +136,7 @@ std::shared_ptr<Opm::EclipseWriter> createEclipseWriter(std::shared_ptr<const Op
     const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid();
 
     std::shared_ptr<Opm::EclipseWriter> eclipseWriter = std::make_shared<Opm::EclipseWriter>(params,
+                                                                                             deck,
                                                                                              eclipseState,
                                                                                              phaseUsage,
                                                                                              ourFinerUnstructuredGrid.number_of_cells,

tests/test_EclipseWriter.cpp

@@ -85,6 +85,7 @@ void createEclipseWriter(const char *deckString)
 
     Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck);
     eclWriter.reset(new Opm::EclipseWriter(params,
+                                           deck,
                                            eclipseState,
                                            phaseUsage,
                                            ourFinerUnstructuredGrid.number_of_cells,

tests/test_writenumwells.cpp

@@ -151,6 +151,7 @@ Opm::EclipseWriterPtr createEclipseWriter(Opm::DeckConstPtr deck,
   const Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck);
 
   Opm::EclipseWriterPtr eclWriter(new Opm::EclipseWriter(params,
+                                                         deck,
                                                          eclipseState,
                                                          phaseUsage,
                                                          eclipseState->getEclipseGrid()->getCartesianSize(),