diff --git a/examples/sim_poly2p_comp_reorder.cpp b/examples/sim_poly2p_comp_reorder.cpp index 2a24dfbba..e91d9d13c 100644 --- a/examples/sim_poly2p_comp_reorder.cpp +++ b/examples/sim_poly2p_comp_reorder.cpp @@ -86,10 +86,10 @@ try // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); - boost::scoped_ptr deck; boost::scoped_ptr grid; boost::scoped_ptr props; boost::scoped_ptr rock_comp; + Opm::DeckConstPtr deck; EclipseStateConstPtr eclipseState; PolymerBlackoilState state; Opm::PolymerProperties poly_props; @@ -99,28 +99,28 @@ try if (use_deck) { std::string deck_filename = param.get("deck_filename"); ParserPtr parser(new Opm::Parser()); - eclipseState.reset(new Opm::EclipseState(parser->parseFile(deck_filename))); + deck = parser->parseFile(deck_filename); + eclipseState.reset(new Opm::EclipseState(deck)); - deck.reset(new EclipseGridParser(deck_filename)); // Grid init - grid.reset(new GridManager(*deck)); + grid.reset(new GridManager(deck)); // Rock and fluid init - props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid())); + props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid())); // check_well_controls = param.getDefault("check_well_controls", false); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. - rock_comp.reset(new RockCompressibility(*deck)); + rock_comp.reset(new RockCompressibility(deck)); // Gravity. - gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; // Init state variables (saturation and pressure). if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { - initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } initBlackoilSurfvol(*grid->c_grid(), *props, state); // Init polymer properties. - poly_props.readFromDeck(*deck); + poly_props.readFromDeck(deck); } else { // Grid init. const int nx = param.getDefault("nx", 100); @@ -246,8 +246,7 @@ try std::cout << "\n\n================ Starting main simulation loop ===============\n" - << " (number of epochs: " - << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + << std::flush; SimulatorReport rep; if (!use_deck) { @@ -277,49 +276,34 @@ try // With a deck, we may have more epochs etc. WellState well_state; int step = 0; + Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck)); SimulatorTimer simtimer; - // Use timer for last epoch to obtain total time. - deck->setCurrentEpoch(deck->numberOfEpochs() - 1); - simtimer.init(*deck); - const double total_time = simtimer.totalTime(); - // Check for WPOLYMER presence in last epoch to decide + simtimer.init(timeMap); + // Check for WPOLYMER presence in last report step to decide // polymer injection control type. - const bool use_wpolymer = deck->hasField("WPOLYMER"); + const bool use_wpolymer = deck->hasKeyword("WPOLYMER"); if (use_wpolymer) { if (param.has("poly_start_days")) { OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. " "You seem to be trying to control it via parameter poly_start_days (etc.) as well."); } } - for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { - // Set epoch index. - deck->setCurrentEpoch(epoch); + for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) { + simtimer.setCurrentStepNum(reportStepIdx); - // Update the timer. - if (deck->hasField("TSTEP")) { - simtimer.init(*deck); - } else { - if (epoch != 0) { - OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); - } - simtimer.init(param); - } - simtimer.setCurrentStepNum(step); - simtimer.setTotalTime(total_time); - - // Report on start of epoch. - std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" - << "\n (number of steps: " + // Report on start of report step. + std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------" + << "\n (number of remaining steps: " << simtimer.numSteps() - step << ")\n\n" << std::flush; // Create new wells, polymer inflow controls. - WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability()); + WellsManager wells(eclipseState , reportStepIdx , *grid->c_grid(), props->permeability()); boost::scoped_ptr polymer_inflow; if (use_wpolymer) { if (wells.c_wells() == 0) { OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); } - polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells())); + polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells())); } else { polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day, param.getDefault("poly_end_days", 800.0)*Opm::unit::day, @@ -327,9 +311,9 @@ try } // @@@ HACK: we should really make a new well state and - // properly transfer old well state to it every epoch, + // properly transfer old well state to it every report step, // since number of wells may change etc. - if (epoch == 0) { + if (reportStepIdx == 0) { well_state.init(wells.c_wells(), state); } @@ -345,7 +329,7 @@ try bcs.c_bcs(), linsolver, grav); - if (epoch == 0) { + if (reportStepIdx == 0) { warnIfUnusedParams(param); } SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); diff --git a/examples/sim_poly2p_incomp_reorder.cpp b/examples/sim_poly2p_incomp_reorder.cpp index bdd1797e6..98400a578 100644 --- a/examples/sim_poly2p_incomp_reorder.cpp +++ b/examples/sim_poly2p_incomp_reorder.cpp @@ -86,7 +86,7 @@ try // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); - boost::scoped_ptr deck; + Opm::DeckConstPtr deck; boost::scoped_ptr grid; boost::scoped_ptr props; boost::scoped_ptr rock_comp; @@ -99,27 +99,27 @@ try if (use_deck) { std::string deck_filename = param.get("deck_filename"); ParserPtr parser(new Opm::Parser()); - eclipseState.reset(new Opm::EclipseState(parser->parseFile(deck_filename))); + deck = parser->parseFile(deck_filename); + eclipseState.reset(new Opm::EclipseState(deck)); - deck.reset(new EclipseGridParser(deck_filename)); // Grid init - grid.reset(new GridManager(*deck)); + grid.reset(new GridManager(deck)); // Rock and fluid init - props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid())); + props.reset(new IncompPropertiesFromDeck(deck, *grid->c_grid())); // check_well_controls = param.getDefault("check_well_controls", false); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. - rock_comp.reset(new RockCompressibility(*deck)); + rock_comp.reset(new RockCompressibility(deck)); // Gravity. - gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; // Init state variables (saturation and pressure). if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); } else { - initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); } // Init polymer properties. - poly_props.readFromDeck(*deck); + poly_props.readFromDeck(deck); } else { // Grid init. const int nx = param.getDefault("nx", 100); @@ -250,8 +250,7 @@ try std::cout << "\n\n================ Starting main simulation loop ===============\n" - << " (number of epochs: " - << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + << std::flush; SimulatorReport rep; if (!use_deck) { @@ -282,49 +281,34 @@ try WellState well_state; int step = 0; + Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck)); SimulatorTimer simtimer; - // Use timer for last epoch to obtain total time. - deck->setCurrentEpoch(deck->numberOfEpochs() - 1); - simtimer.init(*deck); - const double total_time = simtimer.totalTime(); + simtimer.init(timeMap); // Check for WPOLYMER presence in last epoch to decide // polymer injection control type. - const bool use_wpolymer = deck->hasField("WPOLYMER"); + const bool use_wpolymer = deck->hasKeyword("WPOLYMER"); if (use_wpolymer) { if (param.has("poly_start_days")) { OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. " "You seem to be trying to control it via parameter poly_start_days (etc.) as well."); } } - for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { - // Set epoch index. - deck->setCurrentEpoch(epoch); + for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) { + simtimer.setCurrentStepNum(reportStepIdx); - // Update the timer. - if (deck->hasField("TSTEP")) { - simtimer.init(*deck); - } else { - if (epoch != 0) { - OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); - } - simtimer.init(param); - } - simtimer.setCurrentStepNum(step); - simtimer.setTotalTime(total_time); - - // Report on start of epoch. - std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" - << "\n (number of steps: " + // Report on start of report step. + std::cout << "\n\n-------------- Starting report step " << reportStepIdx << " --------------" + << "\n (number of remaining steps: " << simtimer.numSteps() - step << ")\n\n" << std::flush; // Create new wells, polymer inflow controls. - WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability()); + WellsManager wells(eclipseState , reportStepIdx , *grid->c_grid(), props->permeability()); boost::scoped_ptr polymer_inflow; if (use_wpolymer) { if (wells.c_wells() == 0) { OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); } - polymer_inflow.reset(new PolymerInflowFromDeck(*deck, *wells.c_wells(), props->numCells())); + polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells())); } else { polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day, param.getDefault("poly_end_days", 800.0)*Opm::unit::day, @@ -332,9 +316,9 @@ try } // @@@ HACK: we should really make a new well state and - // properly transfer old well state to it every epoch, + // properly transfer old well state to it every report step, // since number of wells may change etc. - if (epoch == 0) { + if (reportStepIdx == 0) { well_state.init(wells.c_wells(), state); } @@ -350,7 +334,7 @@ try bcs.c_bcs(), linsolver, grav); - if (epoch == 0) { + if (reportStepIdx == 0) { warnIfUnusedParams(param); } SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); diff --git a/examples/test_singlecellsolves.cpp b/examples/test_singlecellsolves.cpp index 4ff5fb1d4..90f402a76 100644 --- a/examples/test_singlecellsolves.cpp +++ b/examples/test_singlecellsolves.cpp @@ -69,8 +69,6 @@ try param.disableOutput(); // std::cout << "--------------- Reading parameters ---------------" << std::endl; - // If we have a "deck_filename", grid and props will be read from that. - boost::scoped_ptr deck; boost::scoped_ptr grid; boost::scoped_ptr props; PolymerState state; diff --git a/opm/polymer/PolymerInflow.cpp b/opm/polymer/PolymerInflow.cpp index 7c4e4b507..7f8a440dc 100644 --- a/opm/polymer/PolymerInflow.cpp +++ b/opm/polymer/PolymerInflow.cpp @@ -18,8 +18,8 @@ */ #include -#include #include +#include #include namespace Opm @@ -60,19 +60,19 @@ namespace Opm /// Constructor. /// @param[in] deck Input deck expected to contain WPOLYMER. - PolymerInflowFromDeck::PolymerInflowFromDeck(const EclipseGridParser& deck, + PolymerInflowFromDeck::PolymerInflowFromDeck(Opm::DeckConstPtr deck, const Wells& wells, const int num_cells) : sparse_inflow_(num_cells) { - if (!deck.hasField("WPOLYMER")) { + if (!deck->hasKeyword("WPOLYMER")) { OPM_MESSAGE("PolymerInflowFromDeck initialized without WPOLYMER in current epoch."); return; } // Extract concentrations and put into cell->concentration map. - const std::vector& wpl = deck.getWPOLYMER().wpolymer_; - const int num_wpl = wpl.size(); + Opm::DeckKeywordConstPtr wpolymerKeyword = deck->getKeyword("WPOLYMER"); + const int num_wpl = wpolymerKeyword->size(); std::map perfcell_conc; for (int i = 0; i < num_wpl; ++i) { // Only use well name and polymer concentration. @@ -80,16 +80,19 @@ namespace Opm // names. int wix = 0; for (; wix < wells.number_of_wells; ++wix) { - if (wpl[i].well_ == wells.name[wix]) { + if (wpolymerKeyword->getRecord(i)->getItem("WELL")->getString(0) == wells.name[wix]) { break; } } if (wix == wells.number_of_wells) { - OPM_THROW(std::runtime_error, "Could not find a match for well " << wpl[i].well_ << " from WPOLYMER."); + OPM_THROW(std::runtime_error, "Could not find a match for well " + << wpolymerKeyword->getRecord(i)->getItem("WELL")->getString(0) + << " from WPOLYMER."); } for (int j = wells.well_connpos[wix]; j < wells.well_connpos[wix+1]; ++j) { const int perf_cell = wells.well_cells[j]; - perfcell_conc[perf_cell] = wpl[i].polymer_concentration_; + perfcell_conc[perf_cell] = + wpolymerKeyword->getRecord(i)->getItem("POLYMER_CONCENTRATION")->getSIDouble(0); } } diff --git a/opm/polymer/PolymerInflow.hpp b/opm/polymer/PolymerInflow.hpp index f40348ec3..7f512c8ee 100644 --- a/opm/polymer/PolymerInflow.hpp +++ b/opm/polymer/PolymerInflow.hpp @@ -21,17 +21,13 @@ #define OPM_POLYMERINFLOW_HEADER_INCLUDED #include +#include #include struct Wells; namespace Opm { - - class EclipseGridParser; - - - /// @brief Interface for classes encapsulating polymer inflow information. class PolymerInflowInterface { @@ -90,7 +86,7 @@ namespace Opm /// \param[in] deck Input deck expected to contain WPOLYMER. /// \param[in] wells Wells structure. /// \param[in] num_cells Number of cells in grid. - PolymerInflowFromDeck(const EclipseGridParser& deck, + PolymerInflowFromDeck(Opm::DeckConstPtr deck, const Wells& wells, const int num_cells); diff --git a/opm/polymer/PolymerProperties.cpp b/opm/polymer/PolymerProperties.cpp index c2e1c4438..f11902836 100644 --- a/opm/polymer/PolymerProperties.cpp +++ b/opm/polymer/PolymerProperties.cpp @@ -21,6 +21,8 @@ #include #include #include +#include +#include namespace Opm { diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 0482fd46b..4ee4ba9fd 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -20,10 +20,15 @@ #ifndef OPM_POLYMERPROPERTIES_HEADER_INCLUDED #define OPM_POLYMERPROPERTIES_HEADER_INCLUDED +#include +#include +#include +#include +#include +#include #include #include -#include namespace Opm @@ -81,9 +86,9 @@ namespace Opm { } - PolymerProperties(const EclipseGridParser& gridparser) + PolymerProperties(Opm::DeckConstPtr deck) { - readFromDeck(gridparser); + readFromDeck(deck); } void set(double c_max, @@ -116,38 +121,48 @@ namespace Opm shear_vrf_vals_ = shear_vrf_vals; } - void readFromDeck(const EclipseGridParser& gridparser) + void readFromDeck(Opm::DeckConstPtr deck) { - // We assume NTMISC=1 - const std::vector& plymax = gridparser.getPLYMAX().plymax_; - c_max_ = plymax[0]; - const std::vector& tlmixpar = gridparser.getTLMIXPAR().tlmixpar_; - mix_param_ = tlmixpar[0]; + Opm::PlymaxTable plymaxTable(deck->getKeyword("PLYMAX"), /*tableIdx=*/0); + Opm::TlmixparTable tlmixparTable(deck->getKeyword("TLMIXPAR"), /*tableIdx=*/0); + + // We also assume that each table has exactly one row... + assert(plymaxTable.numRows() == 1); + assert(tlmixparTable.numRows() == 1); + + c_max_ = plymaxTable.getPolymerConcentrationColumn()[0]; + mix_param_ = tlmixparTable.getViscosityParameterColumn()[0]; // We assume NTSFUN=1 - const std::vector& plyrock = gridparser.getPLYROCK().plyrock_; - assert(plyrock.size() == 5); - dead_pore_vol_ = plyrock[0]; - res_factor_ = plyrock[1]; - rock_density_ = plyrock[2]; - ads_index_ = static_cast(plyrock[3]); - c_max_ads_ = plyrock[4]; + Opm::PlyrockTable plyrockTable(deck->getKeyword("PLYROCK"), /*tableIdx=*/0); + + // We also assume that each table has exactly one row... + assert(plyrockTable.numRows() == 1); + + dead_pore_vol_ = plyrockTable.getDeadPoreVolumeColumn()[0]; + res_factor_ = plyrockTable.getResidualResistanceFactorColumn()[0]; + rock_density_ = plyrockTable.getRockDensityFactorColumn()[0]; + ads_index_ = static_cast(plyrockTable.getAdsorbtionIndexColumn()[0]); + c_max_ads_ = plyrockTable.getMaxAdsorbtionColumn()[0]; // We assume NTPVT=1 - const PLYVISC& plyvisc = gridparser.getPLYVISC(); - c_vals_visc_ = plyvisc.concentration_; - visc_mult_vals_ = plyvisc.factor_; + Opm::PlyviscTable plyviscTable(deck->getKeyword("PLYVISC"), /*tableIdx=*/0); + + // We also assume that each table has exactly one row... + assert(plyviscTable.numRows() == 1); + + c_vals_visc_[0] = plyviscTable.getPolymerConcentrationColumn()[0]; + visc_mult_vals_[0] = plyviscTable.getViscosityMultiplierColumn()[0]; // We assume NTSFUN=1 - const PLYADS& plyads = gridparser.getPLYADS(); - c_vals_ads_ = plyads.local_concentration_; - ads_vals_ = plyads.adsorbed_concentration_; - - // We assume NTPVT=1 - const PLYSHEAR& plyshear = gridparser.getPLYSHEAR(); - water_vel_vals_ = plyshear.water_velocity_; - shear_vrf_vals_ = plyshear.vrf_; + Opm::PlyadsTable plyadsTable(deck->getKeyword("PLYADS"), /*tableIdx=*/0); + + // We also assume that each table has exactly one row... + assert(plyadsTable.numRows() == 1); + + c_vals_ads_[0] = plyadsTable.getPolymerConcentrationColumn()[0]; + ads_vals_[0] = plyadsTable.getAdsorbedPolymerColumn()[0]; } double cMax() const;