diff --git a/opm/core/io/eclipse/EclipseWriter.cpp b/opm/core/io/eclipse/EclipseWriter.cpp index c13aef46..65b7e400 100644 --- a/opm/core/io/eclipse/EclipseWriter.cpp +++ b/opm/core/io/eclipse/EclipseWriter.cpp @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -34,6 +35,10 @@ #include #include // WellType +#include +#include +#include + #include // to_upper_copy #include #include // path @@ -61,6 +66,30 @@ using namespace Opm::parameter; // namespace start here since we don't want the ERT headers in it namespace { +namespace { +/// Helper function when we don't really want any transformation +/// (The C++ committee removed std::identity because it was "troublesome" (!?!) +static double noConversion (const double& u) { return u; } + + +/// Helper method that can be used in keyword transformation (must carry +/// the barsa argument) +static double toBar (const double& pressure) { + return Opm::unit::convert::to (pressure, Opm::unit::barsa); +} + +/// Helper method that can be used in keyword transformation (must carry +/// the milliDarcy argument) +static double toMilliDarcy (const double& permeability) { + return Opm::unit::convert::to (permeability, Opm::prefix::milli * Opm::unit::darcy); +} + +/// Names of the saturation property for each phase. The order of these +/// names are critical; they must be the same as the BlackoilPhases enum +static const char* SAT_NAMES[] = { "SWAT", "SOIL", "SGAS" }; + +} // anonymous namespace + /// Smart pointer/handle class for ERT opaque types, such as ecl_kw_type*. /// /// \tparam T Type of handle being wrapper @@ -104,53 +133,162 @@ private: static void no_delete (T*) { } }; +// retrieve all data fields in SI units of a deck keyword +std::vector getAllSiDoubles_(Opm::DeckKeywordConstPtr keywordPtr) +{ + std::vector retBuff; + for (unsigned i = 0; i < keywordPtr->size(); ++i) { + Opm::DeckRecordConstPtr recordPtr(keywordPtr->getRecord(i)); + for (unsigned j = 0; j < recordPtr->size(); ++j) { + Opm::DeckItemConstPtr itemPtr(recordPtr->getItem(j)); + for (unsigned k = 0; k < itemPtr->size(); ++k) { + retBuff.push_back(itemPtr->getSIDouble(k)); + } + } + } + return retBuff; +} + +// retrieve all integer data fields of a deck keyword +std::vector getAllIntegers_(Opm::DeckKeywordConstPtr keywordPtr) +{ + std::vector retBuff; + for (unsigned i = 0; i < keywordPtr->size(); ++i) { + Opm::DeckRecordConstPtr recordPtr(keywordPtr->getRecord(i)); + for (unsigned j = 0; j < recordPtr->size(); ++j) { + Opm::DeckItemConstPtr itemPtr(recordPtr->getItem(j)); + for (unsigned k = 0; k < itemPtr->size(); ++k) { + retBuff.push_back(itemPtr->getInt(k)); + } + } + } + return retBuff; +} + +// throw away the data for all non-active cells in an array +void restrictToActiveCells_(std::vector &data, const std::vector &actnumData) +{ + assert(actnumData.size() == data.size()); + + size_t curActiveIdx = 0; + for (size_t curIdx = 0; curIdx < data.size(); ++curIdx) { + if (!actnumData[curIdx]) + continue; // ignore non-active cells + + assert(curActiveIdx <= curIdx); + data[curActiveIdx] = data[curIdx]; + ++ curActiveIdx; + } + + data.resize(curActiveIdx); +} + +// throw away the data for all non-active cells in an array. (this is +// the variant of the function which takes an UnstructuredGrid object.) +void restrictToActiveCells_(std::vector &data, const UnstructuredGrid &grid) +{ + if (!grid.global_cell) + // if there is no active -> global mapping, all cells + // are considered active + return; + + // activate those cells that are actually there + for (int i = 0; i < grid.number_of_cells; ++i) { + // make sure that global cell indices are always at least as + // large as the active one and that the global cell indices + // are in increasing order. the latter might become + // problematic if cells are extensively re-ordered, but that + // does not seem to be the case so far + assert(grid.global_cell[i] >= i); + assert(i == 0 || grid.global_cell[i - 1] < grid.global_cell[i]); + + data[i] = data[grid.global_cell[i]]; + } + data.resize(grid.number_of_cells); +} + +// convert the units of an array +template +void convertUnit_(std::vector &data, TransferFunction &transferFn) +{ + for (size_t curIdx = 0; curIdx < data.size(); ++curIdx) { + data[curIdx] = transferFn(data[curIdx]); + } +} + +// extract a sub-array of a larger one which represents multiple +// striped ones +void extractFromStripedData_(std::vector &data, + int offset, + int stride) +{ + size_t tmpIdx = 0; + for (size_t curIdx = offset; curIdx < data.size(); curIdx += stride) { + assert(tmpIdx <= curIdx); + data[tmpIdx] = data[curIdx]; + ++tmpIdx; + } + // shirk the result + data.resize(tmpIdx); +} + +// enclosure of the current grid in a Cartesian space +int getCartesianSize_(const UnstructuredGrid& grid) { + const int nx = grid.cartdims[0]; + const int ny = grid.cartdims[1]; + const int nz = grid.cartdims[2]; + return nx * ny * nz; +} + +void getActiveCells_(const UnstructuredGrid& grid, + std::vector & actnum) +{ + // we must fill the Cartesian grid with flags + const int size = getCartesianSize_(grid); + + // if we don't have a global_cells field, then assume that all + // grid cells is active + if (!grid.global_cell) { + if (grid.number_of_cells != size) { + OPM_THROW (std::runtime_error, + "No ACTNUM map but grid size != Cartesian size"); + } + actnum.assign (size, 1); + } + else { + // start out with entire map being inactive + actnum.assign (size, 0); + + // activate those cells that are actually there + for (int i = 0; i < grid.number_of_cells; ++i) { + actnum[grid.global_cell[i]] = 1; + } + } +} + /** * Eclipse "keyword" (i.e. named data) for a vector. (This class is * different from EclKW in the constructors it provide). */ template struct EclipseKeyword : public EclipseHandle { - EclipseKeyword (const std::string& name, /// identification - const std::vector& data, /// array holding values - const int offset = 0, /// distance to first - const int stride = 1) /// distance between each - - // allocate handle and put in smart pointer base class - : EclipseHandle ( - ecl_kw_alloc (name.c_str(), data.size (), type ()), - ecl_kw_free) { - copyData (data, &no_conversion, offset, stride); - } - - /// Special initialization from double-precision array which - /// automatically invokes a version of the copy function which - /// downcasts. This is really only applicable to the T = float - /// template instance. - /// The data and name parameters are switched in this version - /// so that it doesn't conflict with the one above in the case - /// of T = double. - EclipseKeyword (const std::vector& data, - const std::string& name, - double (* const transf)(const double&), - const int offset = 0, - const int stride = 1); - - /// Convenience constructor that gets the set of data - /// from the samely named item in the parser + /// Special initialization from double-precision array. EclipseKeyword (const std::string& name, - const EclipseGridParser& parser) - // allocate handle and put in smart pointer base class - // notice dataSize is called both here *and* in copyData, - // but GCC 4.4 doesn't support delegating constructors, so - // we cannot avoid this without otherwise using a member - : EclipseHandle ( - ecl_kw_alloc (name.c_str(), - dataSize (parser.getValue (name)), - type ()), - ecl_kw_free) { - const std::vector & data = parser.getValue (name); - copyData (data, &no_conversion, 0, 1); - } + const std::vector& data) + : EclipseHandle(ecl_kw_alloc(name.c_str(), + data.size(), + type()), + ecl_kw_free) + { copyData (data, &noConversion, /*offset=*/0, /*stride=*/1); } + + /// Initialization from integer array. + EclipseKeyword (const std::string& name, + const std::vector& data) + : EclipseHandle(ecl_kw_alloc(name.c_str(), + data.size(), + type()), + ecl_kw_free) + { copyData (data, &noConversion, /*offset=*/0, /*stride=*/1); } /// Constructor for optional fields EclipseKeyword (const std::string& name) @@ -158,22 +296,24 @@ struct EclipseKeyword : public EclipseHandle { static_cast (name); } + // constructor to keep compatibility with the old eclipse parser + EclipseKeyword (const std::string& name, + const EclipseGridParser& parser); + // GCC 4.4 doesn't generate these constructors for us; provide the // default implementation explicitly here instead EclipseKeyword (EclipseKeyword&& rhs) - : EclipseHandle (std::move (rhs)) { } - EclipseKeyword& operator= (EclipseKeyword&& rhs) { + : EclipseHandle (std::move (rhs)) + { } + + EclipseKeyword& operator= (EclipseKeyword&& rhs) + { EclipseHandle ::operator= (std::move(rhs)); return *this; } EclipseKeyword (const EclipseKeyword&) = delete; EclipseKeyword& operator= (const EclipseKeyword&) = delete; - /// Helper function when we don't really want any transformation - /// (The C++ committee removed std::identity because it was "troublesome" (!?!) - template - static U no_conversion (const U& u) { return u; } - private: /// Map the C++ data type (given by T) to an Eclipse type enum static ecl_type_enum type (); @@ -181,7 +321,7 @@ private: /// Helper function that is the meat of the constructor template void copyData (const std::vector & data, - U (* const transf)(const U&), + double (* const transf)(const double&), const int offset, const int stride) { // number of elements to take @@ -216,6 +356,32 @@ private: const int recs = (num - 1 - offset) / stride + 1; return recs; } + + int dataSize (Opm::DeckKeywordConstPtr keyword, + const int offset = 0, + const int stride = 1) + { + int numFlatItems = 0; + Opm::DeckRecordConstPtr record = keyword->getRecord(0); + for (unsigned itemIdx = 0; itemIdx < record->size(); ++itemIdx) { + numFlatItems += record->getItem(itemIdx)->size(); + } + + // range cannot start outside of data set + assert(offset >= 0 && offset < numFlatItems); + + // don't jump out of the set when trying to + assert(stride > 0 && stride < numFlatItems - offset); + + // number of (strided) entries it will provide. the last item + // in the array is num - 1. the last strided item we can pick + // (from recs number of records) is (recs - 1) * stride + offset, + // which must be <= num - 1. we are interested in the maximum + // case where it holds to equals. rearranging the above gives us: + const int recs = (numFlatItems - 1 - offset) / stride + 1; + return recs; + } + }; // specializations for known keyword types @@ -238,23 +404,7 @@ EclipseKeyword ::EclipseKeyword ( type ()), ecl_kw_free) { const std::vector & data = parser.getValue (name); - copyData (data, &no_conversion, 0, 1); -} - -/// Provide only the float version, since that is the one for which -/// we need this conversion (we don't want it for int, for instance) -template <> -EclipseKeyword ::EclipseKeyword ( - const std::vector& data, - const std::string& name, - double (* const transf)(const double&), - const int offset, - const int stride) - // allocate handle and put in smart pointer base class - : EclipseHandle ( - ecl_kw_alloc (name.c_str(), dataSize (data, offset, stride), type ()), - ecl_kw_free) { - copyData (data, transf, offset, stride); + copyData (data, &noConversion, 0, 1); } /** @@ -305,6 +455,30 @@ std::vector parserDim (const EclipseGridParser& parser) { return dim; } +/// Get dimensions of the grid from the parse of the input file +std::vector parserDim (Opm::DeckConstPtr newParserDeck) { + std::vector dim(/* n = */ 3); + // dimensions explicitly given + if (newParserDeck->hasKeyword("SPECGRID")) { + SpecgridWrapper specgrid(newParserDeck->getKeyword("SPECGRID")); + dim = specgrid.numBlocksVector(); + } + // dimensions implicitly given by number of deltas + else if (newParserDeck->hasKeyword("DXV")) { + assert(newParserDeck->hasKeyword("DYV")); + assert(newParserDeck->hasKeyword("DZV")); + dim[0] = newParserDeck->getKeyword("DXV")->getRawDoubleData().size(); + dim[1] = newParserDeck->getKeyword("DYV")->getRawDoubleData().size(); + dim[2] = newParserDeck->getKeyword("DZV")->getRawDoubleData().size(); + } + else { + OPM_THROW(std::runtime_error, + "Only decks featureing either the SPECGRID or the D[XYZ]V keywords " + "are currently supported"); + } + return dim; +} + /// Convert OPM phase usage to ERT bitmask static int phaseMask (const PhaseUsage uses) { return (uses.phase_used [BlackoilPhases::Liquid] ? ECL_OIL_PHASE : 0) @@ -344,6 +518,24 @@ struct EclipseRestart : public EclipseHandle { num_active_cells, phaseMask (uses)); } + + void writeHeader (const SimulatorTimer& timer, + int outputStepIdx, + const PhaseUsage uses, + Opm::DeckConstPtr newParserDeck, + const int num_active_cells) { + const std::vector dim = parserDim (newParserDeck); + ecl_rst_file_fwrite_header (*this, + outputStepIdx, + timer.currentPosixTime(), + Opm::unit::convert::to (timer.simulationTimeElapsed (), + Opm::unit::day), + dim[0], + dim[1], + dim[2], + num_active_cells, + phaseMask (uses)); + } }; /** @@ -369,39 +561,6 @@ private: } }; -// enclosure of the current grid in a Cartesian space -int cart_size (const UnstructuredGrid& grid) { - const int nx = grid.cartdims[0]; - const int ny = grid.cartdims[1]; - const int nz = grid.cartdims[2]; - return nx * ny * nz; -} - -void active_cells (const UnstructuredGrid& grid, - std::vector & actnum) { - // we must fill the Cartesian grid with flags - const int size = cart_size (grid); - - // if we don't have a global_cells field, then assume that all - // grid cells is active - if (!grid.global_cell) { - if (grid.number_of_cells != size) { - OPM_THROW (std::runtime_error, - "No ACTNUM map but grid size != Cartesian size"); - } - actnum.assign (size, 1); - } - else { - // start out with entire map being inactive - actnum.assign (size, 0); - - // activate those cells that are actually there - for (int i = 0; i < grid.number_of_cells; ++i) { - actnum[grid.global_cell[i]] = 1; - } - } -} // active_cells - /** * Representation of an Eclipse grid. */ @@ -424,17 +583,65 @@ struct EclipseGrid : public EclipseHandle { else if (parser.hasField("ZCORN")) { struct grdecl g = parser.get_grdecl (); - EclipseKeyword coord_kw (COORD_KW, parser); - EclipseKeyword zcorn_kw (ZCORN_KW, parser); + auto coordData = parser.getFloatingPointValue(COORD_KW); + auto zcornData = parser.getFloatingPointValue(ZCORN_KW); + + EclipseKeyword coord_kw (COORD_KW, coordData); + EclipseKeyword zcorn_kw (ZCORN_KW, zcornData); // get the actually active cells, after processing std::vector actnum; - active_cells (grid, actnum); + getActiveCells_(grid, actnum); EclipseKeyword actnum_kw (ACTNUM_KW, actnum); EclipseKeyword mapaxes_kw (MAPAXES_KW); if (g.mapaxes) { - mapaxes_kw = std::move (EclipseKeyword (MAPAXES_KW, parser)); + auto mapaxesData = parser.getFloatingPointValue(MAPAXES_KW); + mapaxes_kw = std::move (EclipseKeyword (MAPAXES_KW, mapaxesData)); + } + + return EclipseGrid (g.dims, zcorn_kw, coord_kw, actnum_kw, mapaxes_kw); + } + else { + OPM_THROW(std::runtime_error, + "Can't create an ERT grid (no supported keywords found in deck)"); + } + } + + /// Create a grid based on the keywords available in input file + static EclipseGrid make (Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + if (newParserDeck->hasKeyword("DXV")) { + // make sure that the DYV and DZV keywords are present if the + // DXV keyword is used in the deck... + assert(newParserDeck->hasKeyword("DYV")); + assert(newParserDeck->hasKeyword("DZV")); + + const auto& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData(); + const auto& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData(); + const auto& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData(); + + return EclipseGrid (dxv, dyv, dzv); + } + else if (newParserDeck->hasKeyword("ZCORN")) { + struct grdecl g; + GridManager::createGrdecl(newParserDeck, g); + + auto coordData = getAllSiDoubles_(newParserDeck->getKeyword(COORD_KW)); + auto zcornData = getAllSiDoubles_(newParserDeck->getKeyword(ZCORN_KW)); + EclipseKeyword coord_kw (COORD_KW, coordData); + EclipseKeyword zcorn_kw (ZCORN_KW, zcornData); + + // get the actually active cells, after processing + std::vector actnum; + getActiveCells_(grid, actnum); + EclipseKeyword actnum_kw (ACTNUM_KW, actnum); + + EclipseKeyword mapaxes_kw (MAPAXES_KW); + if (g.mapaxes) { + auto mapaxesData = getAllSiDoubles_(newParserDeck->getKeyword(MAPAXES_KW)); + mapaxes_kw = std::move (EclipseKeyword (MAPAXES_KW, mapaxesData)); } return EclipseGrid (g.dims, zcorn_kw, coord_kw, actnum_kw, mapaxes_kw); @@ -541,15 +748,39 @@ struct EclipseInit : public EclipseHandle { timer.currentPosixTime()); } - void writeKeyword (const std::string& keyword, + void writeHeader (const UnstructuredGrid& grid, + const SimulatorTimer& timer, + Opm::DeckConstPtr newParserDeck, + const PhaseUsage uses) + { + auto dataField = getAllSiDoubles_(newParserDeck->getKeyword(PORO_KW)); + restrictToActiveCells_(dataField, grid); + + EclipseGrid eclGrid = EclipseGrid::make (newParserDeck, grid); + + EclipseKeyword poro (PORO_KW, dataField); + ecl_init_file_fwrite_header (*this, + eclGrid, + poro, + phaseMask (uses), + timer.currentPosixTime ()); + } + + void writeKeyword (const std::string& keywordName, const EclipseGridParser& parser, double (* const transf)(const double&)) { - EclipseKeyword kw (parser.getValue (keyword), - keyword, - transf); + auto data = parser.getValue (keywordName); + convertUnit_(data, transf); + EclipseKeyword kw (keywordName, data); ecl_kw_fwrite (kw, *this); } + void writeKeyword (const std::string& keywordName, const std::vector &data) + { + EclipseKeyword kw (keywordName, data); + ecl_kw_fwrite(kw, *this); + } + // GCC 4.4 doesn't generate these constructors for us; provide the // default implementation explicitly here instead EclipseInit (EclipseInit&& rhs) @@ -591,6 +822,14 @@ struct EclipseSummary : public EclipseHandle { alloc_writer (outputDir, baseName, timer, parser), ecl_sum_free) { } + EclipseSummary (const std::string& outputDir, + const std::string& baseName, + const SimulatorTimer& timer, + Opm::DeckConstPtr newParserDeck) + : EclipseHandle ( + alloc_writer (outputDir, baseName, timer, newParserDeck), + ecl_sum_free) { } + typedef std::unique_ptr var_t; typedef std::vector vars_t; @@ -610,6 +849,10 @@ struct EclipseSummary : public EclipseHandle { void addWells (const EclipseGridParser& parser, const PhaseUsage& uses); + // add rate variables for each of the well in the input file + void addWells (Opm::DeckConstPtr newParserDeck, + const PhaseUsage& uses); + // no inline implementation of this since it depends on the // EclipseWellReport type being completed first void writeTimeStep (const SimulatorTimer& timer, @@ -656,6 +899,27 @@ private: dim[1], dim[2]); } + + /// Helper routine that lets us use local variables to hold + /// intermediate results while filling out the allocations function's + /// argument list. + static ecl_sum_type* alloc_writer (const std::string& outputDir, + const std::string& baseName, + const SimulatorTimer& timer, + Opm::DeckConstPtr newParserDeck) { + boost::filesystem::path casePath (outputDir); + casePath /= boost::to_upper_copy (baseName); + + const std::vector dim = parserDim (newParserDeck); + return ecl_sum_alloc_writer (casePath.string ().c_str (), + false, /* formatted */ + true, /* unified */ + ":", /* join string */ + timer.simulationTimeElapsed (), + dim[0], + dim[1], + dim[2]); + } }; @@ -687,6 +951,29 @@ protected: // producers can be seen as negative injectors , sign_ (type == INJECTOR ? +1. : -1.) { } + EclipseWellReport (const EclipseSummary& summary, /* section to add to */ + Opm::DeckConstPtr newParserDeck, /* well names */ + int whichWell, /* index of well line */ + PhaseUsage uses, /* phases present */ + BlackoilPhases::PhaseIndex phase, /* oil, water or gas */ + WellType type, /* prod. or inj. */ + char aggregation, /* rate or total */ + std::string unit) + : EclipseHandle ( + ecl_sum_add_var (summary, + varName (phase, + type, + aggregation).c_str (), + wellName (newParserDeck, whichWell).c_str (), + /* num = */ 0, + unit.c_str(), + /* defaultValue = */ 0.)) + // save these for when we update the value in a timestep + , index_ (whichWell * uses.num_phases + uses.phase_pos [phase]) + + // producers can be seen as negative injectors + , sign_ (type == INJECTOR ? +1. : -1.) { } + public: /// Allows us to pass this type to ecl_sum_tstep_iset operator int () { @@ -711,6 +998,14 @@ private: return parser.getWELSPECS().welspecs[whichWell].name_; } + /// Get the name associated with this well + std::string wellName (Opm::DeckConstPtr newParserDeck, + int whichWell) + { + Opm::WelspecsWrapper welspecs(newParserDeck->getKeyword("WELSPECS")); + return welspecs.wellName(whichWell); + } + /// Compose the name of the summary variable, e.g. "WOPR" for /// well oil production rate. std::string varName (BlackoilPhases::PhaseIndex phase, @@ -772,6 +1067,22 @@ struct EclipseWellRate : public EclipseWellReport { type, 'R', "SM3/DAY" /* surf. cub. m. per day */ ) { } + + EclipseWellRate (const EclipseSummary& summary, + Opm::DeckConstPtr newParserDeck, + int whichWell, + PhaseUsage uses, + BlackoilPhases::PhaseIndex phase, + WellType type) + : EclipseWellReport (summary, + newParserDeck, + whichWell, + uses, + phase, + type, + 'R', + "SM3/DAY" /* surf. cub. m. per day */ ) { } + virtual double update (const SimulatorTimer& /*timer*/, const WellState& wellState) { // TODO: Why only positive rates? @@ -799,6 +1110,24 @@ struct EclipseWellTotal : public EclipseWellReport { // nothing produced when the reporting starts , total_ (0.) { } + EclipseWellTotal (const EclipseSummary& summary, + Opm::DeckConstPtr newParserDeck, + int whichWell, + PhaseUsage uses, + BlackoilPhases::PhaseIndex phase, + WellType type) + : EclipseWellReport (summary, + newParserDeck, + whichWell, + uses, + phase, + type, + 'T', + "SM3" /* surface cubic meter */ ) + + // nothing produced when the reporting starts + , total_ (0.) { } + virtual double update (const SimulatorTimer& timer, const WellState& wellState) { if (timer.currentStepNum() == 0) { @@ -926,26 +1255,49 @@ EclipseSummary::addWells (const EclipseGridParser& parser, } } -namespace { - -/// Helper method that can be used in keyword transformation (must curry -/// the barsa argument) -static double toBar (const double& pressure) { - return Opm::unit::convert::to (pressure, Opm::unit::barsa); +inline void +EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck, + const PhaseUsage& uses) { + // 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 + Opm::DeckKeywordConstPtr welspecsKeyword = newParserDeck->getKeyword("WELSPECS"); + const int numWells = welspecsKeyword->size(); + for (int phaseCounter = 0; + phaseCounter != BlackoilPhases::MaxNumPhases; + ++phaseCounter) { + const BlackoilPhases::PhaseIndex phase = + static_cast (phaseCounter); + // don't bother with reporting for phases that aren't there + if (!uses.phase_used [phaseCounter]) { + continue; + } + for (size_t typeIndex = 0; + typeIndex < sizeof (WELL_TYPES) / sizeof (WELL_TYPES[0]); + ++typeIndex) { + const WellType type = WELL_TYPES[typeIndex]; + for (int whichWell = 0; whichWell != numWells; ++whichWell) { + // W{O,G,W}{I,P}R + add (std::unique_ptr ( + new EclipseWellRate (*this, + newParserDeck, + whichWell, + uses, + phase, + type))); + // W{O,G,W}{I,P}T + add (std::unique_ptr ( + new EclipseWellTotal (*this, + newParserDeck, + whichWell, + uses, + phase, + type))); + } + } + } } -/// Helper method that can be used in keyword transformation (must curry -/// the milliDarcy argument) -static double toMilliDarcy (const double& permeability) { - return Opm::unit::convert::to (permeability, Opm::prefix::milli * Opm::unit::darcy); -} - -/// Names of the saturation property for each phase. The order of these -/// names are critical; they must be the same as the BlackoilPhases enum -static const char* SAT_NAMES[] = { "SWAT", "SOIL", "SGAS" }; - -} // anonymous namespace - namespace Opm { void EclipseWriter::writeInit(const SimulatorTimer &timer, @@ -957,62 +1309,155 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer, if (!enableOutput_) { return; } + if (newParserDeck_) { + /* Grid files */ + EclipseGrid eclGrid = EclipseGrid::make (newParserDeck_, *grid_); + eclGrid.write (outputDir_, baseName_, /*stepIdx=*/0); - /* Grid files */ - EclipseGrid ecl_grid = EclipseGrid::make (*parser_, *grid_); - ecl_grid.write (outputDir_, baseName_, /*stepIdx=*/0); + EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0); + fortio.writeHeader (*grid_, + timer, + newParserDeck_, + uses_); - EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0); - fortio.writeHeader (ecl_grid, - timer, - *parser_, - uses_); + if (newParserDeck_->hasKeyword("PERM")) { + auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERM")); + convertUnit_(data, toMilliDarcy); + fortio.writeKeyword ("PERM", data); + } - fortio.writeKeyword ("PERMX", *parser_, &toMilliDarcy); - fortio.writeKeyword ("PERMY", *parser_, &toMilliDarcy); - fortio.writeKeyword ("PERMZ", *parser_, &toMilliDarcy); + if (newParserDeck_->hasKeyword("PERMX")) { + auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMX")); + convertUnit_(data, toMilliDarcy); + fortio.writeKeyword ("PERMX", data); + } + if (newParserDeck_->hasKeyword("PERMY")) { + auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMY")); + convertUnit_(data, toMilliDarcy); + fortio.writeKeyword ("PERMY", data); + } + if (newParserDeck_->hasKeyword("PERMZ")) { + auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMZ")); + convertUnit_(data, toMilliDarcy); + fortio.writeKeyword ("PERMZ", data); + } - /* Initial solution (pressure and saturation) */ - writeSolution_(timer, reservoirState); + /* Initial solution (pressure and saturation) */ + writeSolution_(timer, reservoirState); - /* Create summary object (could not do it at construction time, - since it requires knowledge of the start time). */ - summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, *parser_)); - summary_->addWells (*parser_, uses_); + /* Create summary object (could not do it at construction time, + since it requires knowledge of the start time). */ + summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, newParserDeck_)); + summary_->addWells (newParserDeck_, uses_); + } + else { + /* Grid files */ + EclipseGrid ecl_grid = EclipseGrid::make (*parser_, *grid_); + ecl_grid.write (outputDir_, baseName_, /*stepIdx=*/0); + + EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0); + fortio.writeHeader (ecl_grid, + timer, + *parser_, + uses_); + + fortio.writeKeyword ("PERMX", *parser_, &toMilliDarcy); + fortio.writeKeyword ("PERMY", *parser_, &toMilliDarcy); + fortio.writeKeyword ("PERMZ", *parser_, &toMilliDarcy); + + /* Initial solution (pressure and saturation) */ + writeSolution_(timer, reservoirState); + + /* Create summary object (could not do it at construction time, + since it requires knowledge of the start time). */ + summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, *parser_)); + summary_->addWells (*parser_, uses_); + } +} + +void EclipseWriter::activeToGlobalCellData_(std::vector &globalCellsBuf, + const std::vector &activeCellsBuf, + const std::vector &inactiveCellsBuf) const +{ + globalCellsBuf = inactiveCellsBuf; + + // overwrite the values of active cells + for (int activeCellIdx = 0; + activeCellIdx < grid_->number_of_cells; + ++activeCellIdx) + { + int globalCellIdx = grid_->global_cell[activeCellIdx]; + globalCellsBuf[globalCellIdx] = activeCellsBuf[activeCellIdx]; + } } void EclipseWriter::writeSolution_(const SimulatorTimer& timer, const SimulatorState& reservoirState) { - // start writing to files - EclipseRestart rst (outputDir_, - baseName_, - timer, - outputTimeStepIdx_); - rst.writeHeader (timer, - outputTimeStepIdx_, - uses_, - *parser_, - reservoirState.pressure ().size ()); - EclipseSolution sol (rst); + if (newParserDeck_) { + // start writing to files + EclipseRestart rst(outputDir_, baseName_, timer, outputTimeStepIdx_); + rst.writeHeader (timer, outputTimeStepIdx_, uses_, newParserDeck_, reservoirState.pressure().size ()); + EclipseSolution sol (rst); - // write pressure and saturation fields (same as DataMap holds) - // convert the pressures from Pascals to bar because Eclipse - // seems to write bars - sol.add(EclipseKeyword(reservoirState.pressure(), "PRESSURE", &toBar)); + // write out the pressure of the reference phase (whatever + // phase that is...). this is not the most performant solution + // thinkable, but this is also not in the most performance + // critical code path! + std::vector tmp = reservoirState.pressure(); + restrictToActiveCells_(tmp, *grid_); + convertUnit_(tmp, toBar); - for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) { - // Eclipse never writes the oil saturation, so all post-processors - // must calculate this from the other saturations anyway - if (phase == BlackoilPhases::PhaseIndex::Liquid) { - continue; + sol.add(EclipseKeyword("PRESSURE", tmp)); + + for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) { + // Eclipse never writes the oil saturation, so all post-processors + // must calculate this from the other saturations anyway + if (phase == BlackoilPhases::PhaseIndex::Liquid) { + continue; + } + if (uses_.phase_used [phase]) { + tmp = reservoirState.saturation(); + extractFromStripedData_(tmp, + /*offset=*/uses_.phase_pos[phase], + /*stride=*/uses_.num_phases); + sol.add(EclipseKeyword(SAT_NAMES[phase], tmp)); + } } - if (uses_.phase_used [phase]) { - sol.add (EclipseKeyword (reservoirState.saturation(), - SAT_NAMES [phase], - &EclipseKeyword ::no_conversion, - uses_.phase_pos [phase], - uses_.num_phases)); + } + else { + // start writing to files + EclipseRestart rst (outputDir_, + baseName_, + timer, + outputTimeStepIdx_); + rst.writeHeader (timer, + outputTimeStepIdx_, + uses_, + *parser_, + reservoirState.pressure ().size ()); + EclipseSolution sol (rst); + + // write pressure and saturation fields (same as DataMap holds) + // convert the pressures from Pascals to bar because Eclipse + // seems to write bars + auto data = reservoirState.pressure(); + convertUnit_(data, toBar); + sol.add(EclipseKeyword("PRESSURE", data)); + + for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) { + // Eclipse never writes the oil saturation, so all post-processors + // must calculate this from the other saturations anyway + if (phase == BlackoilPhases::PhaseIndex::Liquid) { + continue; + } + if (uses_.phase_used [phase]) { + auto tmp = reservoirState.saturation(); + extractFromStripedData_(tmp, + /*offset=*/uses_.phase_pos[phase], + /*stride=*/uses_.num_phases); + sol.add (EclipseKeyword(SAT_NAMES[phase], tmp)); + } } } @@ -1053,7 +1498,6 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer, // the last step. summary_->writeTimeStep(timer, wellState); } -} // namespace Opm #else namespace Opm { @@ -1087,17 +1531,14 @@ void EclipseWriter::writeTimeStep( "The ERT libraries are required to write ECLIPSE output files."); } -} // namespace Opm - #endif // HAVE_ERT -namespace Opm { - EclipseWriter::EclipseWriter ( const ParameterGroup& params, std::shared_ptr parser, std::shared_ptr grid) : parser_ (parser) + , newParserDeck_(0) , grid_ (grid) , uses_ (phaseUsageFromDeck (*parser)) { @@ -1147,6 +1588,57 @@ EclipseWriter::EclipseWriter ( } } +EclipseWriter::EclipseWriter ( + const ParameterGroup& params, + Opm::DeckConstPtr newParserDeck, + std::shared_ptr grid) + : parser_ (0) + , newParserDeck_(newParserDeck) + , grid_ (grid) + , uses_ (phaseUsageFromDeck (newParserDeck)) +{ + // get the base name from the name of the deck + using boost::filesystem::path; + path deck (params.get ("deck_filename")); + if (boost::to_upper_copy (path (deck.extension ()).string ()) == ".DATA") { + baseName_ = path (deck.stem ()).string (); + } + else { + baseName_ = path (deck.filename ()).string (); + } + + // make uppercase of everything (or otherwise we'll get uppercase + // of some of the files (.SMSPEC, .UNSMRY) and not others + baseName_ = boost::to_upper_copy (baseName_); + + // retrieve the value of the "output" parameter + enableOutput_ = params.getDefault("output", /*defaultValue=*/true); + + // retrieve the interval at which something should get written to + // disk (once every N timesteps) + outputInterval_ = params.getDefault("output_interval", /*defaultValue=*/1); + + // store in current directory if not explicitly set + outputDir_ = params.getDefault("output_dir", "."); + + // set the index of the first time step written to 0... + outputTimeStepIdx_ = 0; + + if (enableOutput_) { + // make sure that the output directory exists, if not try to create it + if (!boost::filesystem::exists(outputDir_)) { + std::cout << "Trying to create directory \"" << outputDir_ << "\" for the simulation output\n"; + boost::filesystem::create_directories(outputDir_); + } + + if (!boost::filesystem::is_directory(outputDir_)) { + OPM_THROW(std::runtime_error, + "The path specified as output directory '" << outputDir_ + << "' is not a directory"); + } + } +} + // default destructor is OK, just need to be defined EclipseWriter::~EclipseWriter() { } diff --git a/opm/core/io/eclipse/EclipseWriter.hpp b/opm/core/io/eclipse/EclipseWriter.hpp index a3451ac3..8b98d3be 100644 --- a/opm/core/io/eclipse/EclipseWriter.hpp +++ b/opm/core/io/eclipse/EclipseWriter.hpp @@ -24,7 +24,10 @@ #include #include +#include + #include +#include #include // std::unique_ptr struct UnstructuredGrid; @@ -60,6 +63,14 @@ public: std::shared_ptr parser, std::shared_ptr grid); + /*! + * \brief Sets the common attributes required to write eclipse + * binary files using ERT. + */ + EclipseWriter(const parameter::ParameterGroup& params, + Opm::DeckConstPtr newParserDeck, + std::shared_ptr grid); + /** * We need a destructor in the compilation unit to avoid the * EclipseSummary being a complete type here. @@ -87,6 +98,7 @@ public: private: std::shared_ptr parser_; + Opm::DeckConstPtr newParserDeck_; std::shared_ptr grid_; bool enableOutput_; int outputInterval_; @@ -96,6 +108,10 @@ private: PhaseUsage uses_; // active phases in the input deck std::shared_ptr summary_; + void activeToGlobalCellData_(std::vector &globalCellsBuf, + const std::vector &activeCellsBuf, + const std::vector &inactiveCellsBuf) const; + /// Write solution field variables (pressure and saturation) void writeSolution_(const SimulatorTimer& timer, const SimulatorState& reservoirState);