From 4aa4108367601b4f609bf82aa6c0c0fbdedf1a7b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 11 Mar 2014 15:31:46 +0100 Subject: [PATCH] EclipseWriter: Don't mingle multiple operations using arguments This means that EclipseKeyword now never processes the data it gets. Instead the data must be explicitly preprocessed by the calling site using the new auxiliary functions "convertUnit()", "extractFromStriped()", etc. This approach needs a few additional temporary copies of the data, but given the facts that readability of this code is much better using this approach, and that EclipseWriter is neither a performance- nor a memory critical codepath, I don't care too much about those temporary arrays... --- opm/core/io/eclipse/EclipseWriter.cpp | 333 +++++++++++++++----------- opm/core/io/eclipse/EclipseWriter.hpp | 5 + 2 files changed, 201 insertions(+), 137 deletions(-) diff --git a/opm/core/io/eclipse/EclipseWriter.cpp b/opm/core/io/eclipse/EclipseWriter.cpp index c13aef46..fdc17543 100644 --- a/opm/core/io/eclipse/EclipseWriter.cpp +++ b/opm/core/io/eclipse/EclipseWriter.cpp @@ -61,6 +61,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 +128,131 @@ private: static void no_delete (T*) { } }; + +// 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 +260,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 +285,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 @@ -238,23 +342,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); } /** @@ -369,39 +457,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 +479,21 @@ 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); @@ -541,15 +600,21 @@ struct EclipseInit : public EclipseHandle { timer.currentPosixTime()); } - void writeKeyword (const std::string& keyword, + 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) @@ -926,26 +991,6 @@ 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); -} - -/// 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, @@ -981,6 +1026,22 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer, 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) { @@ -999,7 +1060,9 @@ void EclipseWriter::writeSolution_(const SimulatorTimer& timer, // 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)); + 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 @@ -1008,14 +1071,15 @@ void EclipseWriter::writeSolution_(const SimulatorTimer& timer, continue; } if (uses_.phase_used [phase]) { - sol.add (EclipseKeyword (reservoirState.saturation(), - SAT_NAMES [phase], - &EclipseKeyword ::no_conversion, - uses_.phase_pos [phase], - uses_.num_phases)); + auto tmp = reservoirState.saturation(); + extractFromStripedData_(tmp, + /*offset=*/uses_.phase_pos[phase], + /*stride=*/uses_.num_phases); + sol.add (EclipseKeyword(SAT_NAMES[phase], tmp)); } } + ++outputTimeStepIdx_; } @@ -1053,7 +1117,6 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer, // the last step. summary_->writeTimeStep(timer, wellState); } -} // namespace Opm #else namespace Opm { @@ -1087,12 +1150,8 @@ 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, diff --git a/opm/core/io/eclipse/EclipseWriter.hpp b/opm/core/io/eclipse/EclipseWriter.hpp index a3451ac3..5ec99543 100644 --- a/opm/core/io/eclipse/EclipseWriter.hpp +++ b/opm/core/io/eclipse/EclipseWriter.hpp @@ -25,6 +25,7 @@ #include #include +#include #include // std::unique_ptr struct UnstructuredGrid; @@ -96,6 +97,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);