OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #116 from andlaus/simplify_simulator

Browse Source 
Simplify simulator
This commit is contained in:
Atgeirr Flø Rasmussen 2015-05-29 15:47:06 +02:00
commit bd51ce8dbe
11 changed files with 464 additions and 1276 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
examples/flow_polymer.cpp
View File

@ -201,8 +201,8 @@ try
// Init state variables (saturation and pressure). // Init state variables (saturation and pressure).
if (param.has("init_saturation")) { if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state.blackoilState()); initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state.blackoilState()); initBlackoilSurfvol(*grid->c_grid(), *props, state);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) { if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int np = props->numPhases(); const int np = props->numPhases();
@ -215,10 +215,10 @@ try
} else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) { } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
state.init(*grid->c_grid(), props->numPhases()); state.init(*grid->c_grid(), props->numPhases());
const double grav = param.getDefault("gravity", unit::gravity); const double grav = param.getDefault("gravity", unit::gravity);
initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState()); initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state);
state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0); state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
} else { } else {
initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState()); initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
} }
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props. // The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
@ -268,21 +268,22 @@ try
std::vector<double> threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid()); std::vector<double> threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid());
SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid> simulator(param, SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid>
*grid->c_grid(), simulator(param,
geology, *grid->c_grid(),
*new_props, geology,
polymer_props_ad, *new_props,
rock_comp->isActive() ? rock_comp.get() : 0, polymer_props_ad,
*fis_solver, rock_comp->isActive() ? rock_comp.get() : 0,
grav, *fis_solver,
deck->hasKeyword("DISGAS"), grav,
deck->hasKeyword("VAPOIL"), deck->hasKeyword("DISGAS"),
polymer, deck->hasKeyword("VAPOIL"),
eclipseState, polymer,
outputWriter, eclipseState,
deck, outputWriter,
threshold_pressures); deck,
threshold_pressures);
if (!schedule->initOnly()){ if (!schedule->initOnly()){
std::cout << "\n\n================ Starting main simulation loop ===============\n" std::cout << "\n\n================ Starting main simulation loop ===============\n"

104
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -58,8 +58,13 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp> #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp> #include <opm/autodiff/GeoProps.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp> #include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp> #include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp> #include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp> #include <boost/algorithm/string.hpp>
@ -112,9 +117,48 @@ try
double gravity[3] = { 0.0 }; double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename"); std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr newParser(new Opm::Parser()); // Write parameters used for later reference.
Opm::DeckConstPtr deck = newParser->parseFile(deck_filename); bool output = param.getDefault("output", true);
std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck)); std::string output_dir;
if (output) {
// Create output directory if needed.
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
std::cerr << "Creating directories failed: " << fpath << std::endl;
return EXIT_FAILURE;
}
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::string logFile = output_dir + "/LOGFILE.txt";
Opm::ParserPtr parser(new Opm::Parser());
{
std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
Opm::OpmLog::addBackend( "STREAM" , streamLog );
Opm::OpmLog::addBackend( "COUNTER" , counterLog );
}
Opm::DeckConstPtr deck;
std::shared_ptr<EclipseState> eclipseState;
try {
deck = parser->parseFile(deck_filename);
Opm::checkDeck(deck);
eclipseState.reset(new Opm::EclipseState(deck));
}
catch (const std::invalid_argument& e) {
std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
std::cerr << "Exception caught: " << e.what() << std::endl;
return EXIT_FAILURE;
}
// Grid init // Grid init
std::vector<double> porv; std::vector<double> porv;
if (eclipseState->hasDoubleGridProperty("PORV")) { if (eclipseState->hasDoubleGridProperty("PORV")) {
@ -123,20 +167,23 @@ try
grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
auto &cGrid = *grid->c_grid(); auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
Opm::EclipseWriter outputWriter(param, Opm::BlackoilOutputWriter outputWriter(cGrid,
eclipseState, param,
pu, eclipseState,
cGrid.number_of_cells, pu );
cGrid.global_cell);
// Rock and fluid init // Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid())); new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
PolymerProperties polymer_props(deck, eclipseState); PolymerProperties polymer_props(deck, eclipseState);
PolymerPropsAd polymer_props_ad(polymer_props); PolymerPropsAd polymer_props_ad(polymer_props);
// Rock compressibility. // Rock compressibility.
rock_comp.reset(new RockCompressibility(deck, eclipseState)); rock_comp.reset(new RockCompressibility(deck, eclipseState));
// Gravity. // Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure). // Init state variables (saturation and pressure).
if (param.has("init_saturation")) { if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
@ -155,22 +202,6 @@ try
fis_solver.reset(new NewtonIterationBlackoilSimple(param)); fis_solver.reset(new NewtonIterationBlackoilSimple(param));
} }
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::string output_dir;
if (output) {
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
param.writeParam(output_dir + "/simulation.param");
}
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap()); Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer; SimulatorTimer simtimer;
simtimer.init(timeMap); simtimer.init(timeMap);
@ -193,17 +224,18 @@ try
SimulatorReport fullReport; SimulatorReport fullReport;
// Create and run simulator. // Create and run simulator.
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
SimulatorFullyImplicitCompressiblePolymer simulator(param, SimulatorFullyImplicitCompressiblePolymer<UnstructuredGrid>
*grid->c_grid(), simulator(param,
geology, *grid->c_grid(),
*new_props, geology,
polymer_props_ad, *new_props,
rock_comp->isActive() ? rock_comp.get() : 0, polymer_props_ad,
eclipseState, rock_comp->isActive() ? rock_comp.get() : 0,
outputWriter, eclipseState,
deck, outputWriter,
*fis_solver, deck,
grav); *fis_solver,
grav);
fullReport= simulator.run(simtimer, state); fullReport= simulator.run(simtimer, state);
std::cout << "\n\n================ End of simulation ===============\n\n"; std::cout << "\n\n================ End of simulation ===============\n\n";

2
opm/polymer/CompressibleTpfaPolymer.cpp
View File

@ -90,7 +90,7 @@ namespace Opm
{ {
c_ = &state.concentration(); c_ = &state.concentration();
cmax_ = &state.maxconcentration(); cmax_ = &state.maxconcentration();
CompressibleTpfa::solve(dt, state.blackoilState(), well_state); CompressibleTpfa::solve(dt, state, well_state);
} }
/// Compute per-solve dynamic properties. /// Compute per-solve dynamic properties.

39
opm/polymer/PolymerBlackoilState.hpp
View File

@ -30,7 +30,7 @@ namespace Opm
/// Simulator state for a compressible two-phase simulator with polymer. /// Simulator state for a compressible two-phase simulator with polymer.
/// We use the Blackoil state parameters. /// We use the Blackoil state parameters.
class PolymerBlackoilState class PolymerBlackoilState : public BlackoilState
{ {
public: public:
void init(const UnstructuredGrid& g, int num_phases) void init(const UnstructuredGrid& g, int num_phases)
@ -40,53 +40,18 @@ namespace Opm
void init(int number_of_cells, int number_of_faces, int num_phases) void init(int number_of_cells, int number_of_faces, int num_phases)
{ {
state_blackoil_.init(number_of_cells, number_of_faces, num_phases); BlackoilState::init(number_of_cells, number_of_faces, num_phases);
concentration_.resize(number_of_cells, 0.0); concentration_.resize(number_of_cells, 0.0);
cmax_.resize(number_of_cells, 0.0); cmax_.resize(number_of_cells, 0.0);
} }
int numPhases() const
{
return state_blackoil_.numPhases();
}
enum ExtremalSat { MinSat = BlackoilState::MinSat, MaxSat = BlackoilState::MaxSat };
void setFirstSat(const std::vector<int>& cells,
const Opm::BlackoilPropertiesInterface& props,
ExtremalSat es)
{
// A better solution for embedding BlackoilState::ExtremalSat could perhaps
// be found, to avoid the cast.
state_blackoil_.setFirstSat(cells, props, static_cast<BlackoilState::ExtremalSat>(es));
}
std::vector<double>& pressure () { return state_blackoil_.pressure(); }
std::vector<double>& temperature () { return state_blackoil_.temperature(); }
std::vector<double>& surfacevol () { return state_blackoil_.surfacevol(); }
std::vector<double>& facepressure() { return state_blackoil_.facepressure(); }
std::vector<double>& faceflux () { return state_blackoil_.faceflux(); }
std::vector<double>& saturation () { return state_blackoil_.saturation(); }
std::vector<double>& gasoilratio () { return state_blackoil_.gasoilratio(); }
std::vector<double>& rv () { return state_blackoil_.rv(); }
std::vector<double>& concentration() { return concentration_; } std::vector<double>& concentration() { return concentration_; }
std::vector<double>& maxconcentration() { return cmax_; } std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& pressure () const { return state_blackoil_.pressure(); }
const std::vector<double>& temperature () const { return state_blackoil_.temperature(); }
const std::vector<double>& surfacevol () const { return state_blackoil_.surfacevol(); }
const std::vector<double>& facepressure() const { return state_blackoil_.facepressure(); }
const std::vector<double>& faceflux () const { return state_blackoil_.faceflux(); }
const std::vector<double>& saturation () const { return state_blackoil_.saturation(); }
const std::vector<double>& gasoilratio() const { return state_blackoil_.gasoilratio(); }
const std::vector<double>& rv () const { return state_blackoil_.rv(); }
const std::vector<double>& concentration() const { return concentration_; } const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& maxconcentration() const { return cmax_; } const std::vector<double>& maxconcentration() const { return cmax_; }
BlackoilState& blackoilState() { return state_blackoil_; }
const BlackoilState& blackoilState() const { return state_blackoil_; }
private: private:
BlackoilState state_blackoil_;
std::vector<double> concentration_; std::vector<double> concentration_;
std::vector<double> cmax_; std::vector<double> cmax_;
}; };

26
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -195,13 +195,14 @@ namespace {
void int
FullyImplicitCompressiblePolymerSolver:: FullyImplicitCompressiblePolymerSolver::
step(const double dt, step(const double dt,
PolymerBlackoilState& x , PolymerBlackoilState& x ,
WellStateFullyImplicitBlackoil& xw, WellStateFullyImplicitBlackoilPolymer& xw)
const std::vector<double>& polymer_inflow)
{ {
const std::vector<double>& polymer_inflow = xw.polymerInflow();
// Initial max concentration of this time step from PolymerBlackoilState. // Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_)); cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
@ -215,7 +216,7 @@ namespace {
const double r0 = residualNorm(); const double r0 = residualNorm();
const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm<Eigen::Infinity>(); const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm<Eigen::Infinity>();
int it = 0; int it = 0;
std::cout << "\nIteration Residual Polymer Res\n" std::cout << "\nIteration Residual Polymer Res\n"
<< std::setw(9) << it << std::setprecision(9) << std::setw(9) << it << std::setprecision(9)
<< std::setw(18) << r0 << std::setprecision(9) << std::setw(18) << r0 << std::setprecision(9)
@ -224,6 +225,9 @@ namespace {
while (resTooLarge && (it < maxit)) { while (resTooLarge && (it < maxit)) {
const V dx = solveJacobianSystem(); const V dx = solveJacobianSystem();
// update the number of linear iterations used.
linearIterations_ += linsolver_.iterations();
updateState(dx, x, xw); updateState(dx, x, xw);
assemble(dt, x, xw, polymer_inflow); assemble(dt, x, xw, polymer_inflow);
@ -233,6 +237,7 @@ namespace {
resTooLarge = (r > atol) && (r > rtol*r0); resTooLarge = (r > atol) && (r > rtol*r0);
it += 1; it += 1;
newtonIterations_ += 1;
std::cout << std::setw(9) << it << std::setprecision(9) std::cout << std::setw(9) << it << std::setprecision(9)
<< std::setw(18) << r << std::setprecision(9) << std::setw(18) << r << std::setprecision(9)
<< std::setw(18) << rr_polymer << std::endl; << std::setw(18) << rr_polymer << std::endl;
@ -240,16 +245,25 @@ namespace {
if (resTooLarge) { if (resTooLarge) {
std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
return -1;
// OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
} }
// Update max concentration. // Update max concentration.
computeCmax(x); computeCmax(x);
return it;
} }
int FullyImplicitCompressiblePolymerSolver::newtonIterations() const
{
return newtonIterations_;
}
int FullyImplicitCompressiblePolymerSolver::linearIterations() const
{
return linearIterations_;
}
FullyImplicitCompressiblePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant() FullyImplicitCompressiblePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant()
: accum(2, ADB::null()) : accum(2, ADB::null())

13
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -27,6 +27,7 @@
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp> #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/polymer/PolymerProperties.hpp> #include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp> #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
struct UnstructuredGrid; struct UnstructuredGrid;
@ -79,11 +80,13 @@ namespace Opm {
/// \param[in] state reservoir state /// \param[in] state reservoir state
/// \param[in] wstate well state /// \param[in] wstate well state
/// \param[in] polymer_inflow polymer influx /// \param[in] polymer_inflow polymer influx
void int
step(const double dt, step(const double dt,
PolymerBlackoilState& state , PolymerBlackoilState& state ,
WellStateFullyImplicitBlackoil& wstate, WellStateFullyImplicitBlackoilPolymer& wstate);
const std::vector<double>& polymer_inflow);
int newtonIterations() const;
int linearIterations() const;
private: private:
typedef AutoDiffBlock<double> ADB; typedef AutoDiffBlock<double> ADB;
@ -142,6 +145,10 @@ namespace Opm {
// each of which has size equal to the number of cells. // each of which has size equal to the number of cells.
// The well_eq has size equal to the number of wells. // The well_eq has size equal to the number of wells.
LinearisedBlackoilResidual residual_; LinearisedBlackoilResidual residual_;
unsigned int newtonIterations_;
unsigned int linearIterations_;
// Private methods. // Private methods.
SolutionState SolutionState
constantState(const PolymerBlackoilState& x, constantState(const PolymerBlackoilState& x,

151
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
View File

@ -21,65 +21,92 @@
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#include <memory> #include <opm/autodiff/SimulatorBase.hpp>
#include <vector> #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
struct UnstructuredGrid; #include <opm/core/utility/parameters/ParameterGroup.hpp>
struct Wells; #include <opm/core/utility/ErrorMacros.hpp>
struct FlowBoundaryConditions;
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/NewtonSolver.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
//#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
//#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm namespace Opm
{ {
namespace parameter { class ParameterGroup; } template <class GridT>
class BlackoilPropsAdInterface; class SimulatorFullyImplicitBlackoilPolymer;
class RockCompressibility;
class DerivedGeology;
class NewtonIterationBlackoilInterface;
class SimulatorTimer;
class PolymerBlackoilState;
class WellStateFullyImplicitBlackoil;
class EclipseState;
class BlackoilOutputWriter;
class PolymerPropsAd;
class PolymerInflowInterface;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation. template<class GridT>
template<class T> struct SimulatorTraits<SimulatorFullyImplicitBlackoilPolymer<GridT> >
class SimulatorFullyImplicitBlackoilPolymer
{ {
typedef WellStateFullyImplicitBlackoilPolymer WellState;
typedef PolymerBlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilPolymerModel<Grid> Model;
typedef NewtonSolver<Model> Solver;
};
/// Class collecting all necessary components for a blackoil simulation with polymer
/// injection.
template <class GridT>
class SimulatorFullyImplicitBlackoilPolymer
: public SimulatorBase<SimulatorFullyImplicitBlackoilPolymer<GridT> >
{
typedef SimulatorFullyImplicitBlackoilPolymer<GridT> ThisType;
typedef SimulatorBase<ThisType> BaseType;
typedef SimulatorTraits<ThisType> Traits;
typedef typename Traits::Solver Solver;
public: public:
/// \brief The type of the grid that we use.
typedef T Grid;
/// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
/// parameter (default) effect
/// -----------------------------------------------------------
/// output (true) write output to files?
/// output_dir ("output") output directoty
/// output_interval (1) output every nth step
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal)
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure
/// nl_maxiter (30) max nonlinear iterations in transport
/// nl_tolerance (1e-9) transport solver absolute residual tolerance
/// num_transport_substeps (1) number of transport steps per pressure step
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] geo derived geological properties
/// \param[in] props fluid and rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
/// \param[in] disgas true for dissolved gas option
/// \param[in] vapoil true for vaporized oil option
/// \param[in] eclipse_state
/// \param[in] output_writer
/// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow
SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
const Grid& grid, const GridT& grid,
const DerivedGeology& geo, const DerivedGeology& geo,
BlackoilPropsAdInterface& props, BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props, const PolymerPropsAd& polymer_props,
@ -94,23 +121,21 @@ namespace Opm
Opm::DeckConstPtr& deck, Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face); const std::vector<double>& threshold_pressures_by_face);
/// Run the simulation. std::unique_ptr<Solver> createSolver(const Wells* wells);
/// This will run succesive timesteps until timer.done() is true. It will
/// modify the reservoir and well states. void handleAdditionalWellInflow(SimulatorTimer& timer,
/// \param[in,out] timer governs the requested reporting timesteps WellsManager& wells_manager,
/// \param[in,out] state state of reservoir: pressure, fluxes typename BaseType::WellState& well_state,
/// \param[in,out] well_state state of wells: bhp, perforation rates const Wells* wells);
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
PolymerBlackoilState& state);
private: private:
class Impl; const PolymerPropsAd& polymer_props_;
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. bool has_polymer_;
std::shared_ptr<Impl> pimpl_; DeckConstPtr deck_;
}; };
} // namespace Opm } // namespace Opm
#include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp" #include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED #endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED

660
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
View File

@ -19,599 +19,99 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>. along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/ */
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/NewtonSolver.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
//#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
//#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm namespace Opm
{ {
template<class T> template <class GridT>
class SimulatorFullyImplicitBlackoilPolymer<T>::Impl SimulatorFullyImplicitBlackoilPolymer<GridT>::
SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
const GridT& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face)
: BaseType(param,
grid,
geo,
props,
rock_comp_props,
linsolver,
gravity,
has_disgas,
has_vapoil,
eclipse_state,
output_writer,
threshold_pressures_by_face)
, polymer_props_(polymer_props)
, has_polymer_(has_polymer)
, deck_(deck)
{ {
public:
Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
bool has_disgas,
bool has_vapoil,
bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face);
SimulatorReport run(SimulatorTimer& timer,
PolymerBlackoilState& state);
private:
// Data.
typedef RateConverter::
SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
std::vector<int> > RateConverterType;
const parameter::ParameterGroup param_;
// Observed objects.
const Grid& grid_;
BlackoilPropsAdInterface& props_;
const PolymerPropsAd& polymer_props_;
const RockCompressibility* rock_comp_props_;
const double* gravity_;
// Solvers
const DerivedGeology& geo_;
NewtonIterationBlackoilInterface& solver_;
// Misc. data
std::vector<int> allcells_;
const bool has_disgas_;
const bool has_vapoil_;
const bool has_polymer_;
bool terminal_output_;
// eclipse_state
std::shared_ptr<EclipseState> eclipse_state_;
// output_writer
BlackoilOutputWriter& output_writer_;
Opm::DeckConstPtr& deck_;
RateConverterType rateConverter_;
// Threshold pressures.
std::vector<double> threshold_pressures_by_face_;
void
computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellStateFullyImplicitBlackoilPolymer& xw);
};
template<class T>
SimulatorFullyImplicitBlackoilPolymer<T>::SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face)
{
pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil, has_polymer,
eclipse_state, output_writer, deck, threshold_pressures_by_face));
} }
template <class GridT>
auto SimulatorFullyImplicitBlackoilPolymer<GridT>::
createSolver(const Wells* wells)
-> std::unique_ptr<Solver>
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::run(SimulatorTimer& timer,
PolymerBlackoilState& state)
{ {
return pimpl_->run(timer, state); typedef typename Traits::Model Model;
}
// \TODO: Treat bcs.
template<class T>
SimulatorFullyImplicitBlackoilPolymer<T>::Impl::Impl(const parameter::ParameterGroup& param,
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face)
: param_(param),
grid_(grid),
props_(props),
polymer_props_(polymer_props),
rock_comp_props_(rock_comp_props),
gravity_(gravity),
geo_(geo),
solver_(linsolver),
has_disgas_(has_disgas),
has_vapoil_(has_vapoil),
has_polymer_(has_polymer),
terminal_output_(param.getDefault("output_terminal", true)),
eclipse_state_(eclipse_state),
output_writer_(output_writer),
deck_(deck),
rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
threshold_pressures_by_face_(threshold_pressures_by_face)
{
// Misc init.
const int num_cells = AutoDiffGrid::numCells(grid);
allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells_[cell] = cell;
}
#if HAVE_MPI
if ( terminal_output_ ) {
if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
{
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
// Only rank 0 does print to std::cout
terminal_output_= (info.communicator().rank()==0);
}
}
#endif
}
template<class T>
SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::Impl::run(SimulatorTimer& timer,
PolymerBlackoilState& state)
{
WellStateFullyImplicitBlackoilPolymer prev_well_state;
// Create timers and file for writing timing info.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
typedef T Grid;
typedef BlackoilPolymerModel<Grid> Model;
typedef typename Model::ModelParameters ModelParams; typedef typename Model::ModelParameters ModelParams;
ModelParams modelParams( param_ ); ModelParams modelParams( BaseType::param_ );
typedef NewtonSolver<Model> Solver; typedef NewtonSolver<Model> Solver;
auto model = std::unique_ptr<Model>(new Model(modelParams,
BaseType::grid_,
BaseType::props_,
BaseType::geo_,
BaseType::rock_comp_props_,
polymer_props_,
wells,
BaseType::solver_,
BaseType::has_disgas_,
BaseType::has_vapoil_,
has_polymer_,
BaseType::terminal_output_));
if (!BaseType::threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(BaseType::threshold_pressures_by_face_);
}
typedef typename Solver::SolverParameters SolverParams; typedef typename Solver::SolverParameters SolverParams;
SolverParams solverParams( param_ ); SolverParams solverParams( BaseType::param_ );
return std::unique_ptr<Solver>(new Solver(solverParams, std::move(model)));
//adaptive time stepping
// std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
// if( param_.getDefault("timestep.adaptive", bool(false) ) )
// {
// adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) );
// }
// init output writer
output_writer_.writeInit( timer );
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
if( ! restorefilename.empty() )
{
// -1 means that we'll take the last report step that was written
const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep );
}
unsigned int totalNewtonIterations = 0;
unsigned int totalLinearIterations = 0;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
step_timer.start();
if ( terminal_output_ )
{
timer.report(std::cout);
}
// Create wells and well state.
WellsManager wells_manager(eclipse_state_,
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(grid_),
Opm::UgGridHelpers::globalCell(grid_),
Opm::UgGridHelpers::cartDims(grid_),
Opm::UgGridHelpers::dimensions(grid_),
Opm::UgGridHelpers::cell2Faces(grid_),
Opm::UgGridHelpers::beginFaceCentroids(grid_),
props_.permeability());
const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoilPolymer well_state;
well_state.init(wells, state.blackoilState(), prev_well_state);
// compute polymer inflow
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, eclipse_state_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
well_state.polymerInflow() = polymer_inflow_c;
// write simulation state at the report stage
output_writer_.writeTimeStep( timer, state.blackoilState(), well_state );
// Max oil saturation (for VPPARS), hysteresis update.
props_.updateSatOilMax(state.saturation());
props_.updateSatHyst(state.saturation(), allcells_);
// Compute reservoir volumes for RESV controls.
computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state);
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
Model model(modelParams, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_);
if (!threshold_pressures_by_face_.empty()) {
model.setThresholdPressures(threshold_pressures_by_face_);
}
Solver solver(solverParams, model);
// If sub stepping is enabled allow the solver to sub cycle
// in case the report steps are to large for the solver to converge
//
// \Note: The report steps are met in any case
// \Note: The sub stepping will require a copy of the state variables
// if( adaptiveTimeStepping ) {
// adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ );
// } else {
// solve for complete report step
solver.step(timer.currentStepLength(), state, well_state);
// }
// take time that was used to solve system for this reportStep
solver_timer.stop();
// accumulate the number of Newton and Linear Iterations
totalNewtonIterations += solver.newtonIterations();
totalLinearIterations += solver.linearIterations();
// Report timing.
const double st = solver_timer.secsSinceStart();
if ( terminal_output_ )
{
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
}
stime += st;
if ( output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
// Increment timer, remember well state.
++timer;
prev_well_state = well_state;
}
// Write final simulation state.
output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state );
// Stop timer and create timing report
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
report.total_newton_iterations = totalNewtonIterations;
report.total_linear_iterations = totalLinearIterations;
return report;
} }
namespace SimFIBODetails { template <class GridT>
typedef std::unordered_map<std::string, WellConstPtr> WellMap; void SimulatorFullyImplicitBlackoilPolymer<GridT>::
handleAdditionalWellInflow(SimulatorTimer& timer,
inline WellMap WellsManager& wells_manager,
mapWells(const std::vector<WellConstPtr>& wells) typename BaseType::WellState& well_state,
{ const Wells* wells)
WellMap wmap;
for (std::vector<WellConstPtr>::const_iterator
w = wells.begin(), e = wells.end();
w != e; ++w)
{
wmap.insert(std::make_pair((*w)->name(), *w));
}
return wmap;
}
inline int
resv_control(const WellControls* ctrl)
{
int i, n = well_controls_get_num(ctrl);
bool match = false;
for (i = 0; (! match) && (i < n); ++i) {
match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
}
if (! match) { i = 0; }
return i - 1; // -1 if no match, undo final "++" otherwise
}
inline bool
is_resv(const Wells& wells,
const int w)
{
return (0 <= resv_control(wells.ctrls[w]));
}
inline bool
is_resv(const WellMap& wmap,
const std::string& name,
const std::size_t step)
{
bool match = false;
WellMap::const_iterator i = wmap.find(name);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
match = (wp->isProducer(step) &&
wp->getProductionProperties(step)
.hasProductionControl(WellProducer::RESV))
|| (wp->isInjector(step) &&
wp->getInjectionProperties(step)
.hasInjectionControl(WellInjector::RESV));
}
return match;
}
inline std::vector<int>
resvWells(const Wells* wells,
const std::size_t step,
const WellMap& wmap)
{
std::vector<int> resv_wells;
if( wells )
{
for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
if (is_resv(*wells, w) ||
((wells->name[w] != 0) &&
is_resv(wmap, wells->name[w], step)))
{
resv_wells.push_back(w);
}
}
}
return resv_wells;
}
inline void
historyRates(const PhaseUsage& pu,
const WellProductionProperties& p,
std::vector<double>& rates)
{
assert (! p.predictionMode);
assert (rates.size() ==
std::vector<double>::size_type(pu.num_phases));
if (pu.phase_used[ BlackoilPhases::Aqua ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Aqua ];
rates[i] = p.WaterRate;
}
if (pu.phase_used[ BlackoilPhases::Liquid ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Liquid ];
rates[i] = p.OilRate;
}
if (pu.phase_used[ BlackoilPhases::Vapour ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Vapour ];
rates[i] = p.GasRate;
}
}
} // namespace SimFIBODetails
template <class T>
void
SimulatorFullyImplicitBlackoilPolymer<T>::
Impl::computeRESV(const std::size_t step,
const Wells* wells,
const BlackoilState& x,
WellStateFullyImplicitBlackoilPolymer& xw)
{ {
typedef SimFIBODetails::WellMap WellMap; // compute polymer inflow
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step); if (deck_->hasKeyword("WPOLYMER")) {
const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
if (! resv_wells.empty()) {
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x);
std::vector<double> distr (np);
std::vector<double> hrates(np);
std::vector<double> prates(np);
for (std::vector<int>::const_iterator
rp = resv_wells.begin(), e = resv_wells.end();
rp != e; ++rp)
{
WellControls* ctrl = wells->ctrls[*rp];
const bool is_producer = wells->type[*rp] == PRODUCER;
// RESV control mode, all wells
{
const int rctrl = SimFIBODetails::resv_control(ctrl);
if (0 <= rctrl) {
const std::vector<double>::size_type off = (*rp) * np;
if (is_producer) {
// Convert to positive rates to avoid issues
// in coefficient calculations.
std::transform(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin(), std::negate<double>());
} else {
std::copy(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin());
}
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(prates, fipreg, distr);
well_controls_iset_distr(ctrl, rctrl, & distr[0]);
}
}
// RESV control, WCONHIST wells. A bit of duplicate
// work, regrettably.
if (is_producer && wells->name[*rp] != 0) {
WellMap::const_iterator i = wmap.find(wells->name[*rp]);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
const WellProductionProperties& p =
wp->getProductionProperties(step);
if (! p.predictionMode) {
// History matching (WCONHIST/RESV)
SimFIBODetails::historyRates(pu, p, hrates);
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(hrates, fipreg, distr);
// WCONHIST/RESV target is sum of all
// observed phase rates translated to
// reservoir conditions. Recall sign
// convention: Negative for producers.
const double target =
- std::inner_product(distr.begin(), distr.end(),
hrates.begin(), 0.0);
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
& distr[0], ctrl);
// For WCONHIST/RESV the BHP limit is set to 1 atm.
// TODO: Make it possible to modify the BHP limit using
// the WELTARG keyword
const int ok_bhp =
well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm),
NULL, ctrl);
if (ok_resv != 0 && ok_bhp != 0) {
xw.currentControls()[*rp] = 0;
well_controls_set_current(ctrl, 0);
}
}
}
}
} }
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
} }
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
well_state.polymerInflow() = polymer_inflow_c;
} }
} // namespace Opm } // namespace Opm

473
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp
View File

@ -1,473 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <numeric>
#include <fstream>
#include <iostream>
namespace Opm
{
namespace
{
static void outputStateVtk(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir);
static void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir);
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir);
} // anonymous namespace
class SimulatorFullyImplicitCompressiblePolymer::Impl
{
public:
Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity);
SimulatorReport run(SimulatorTimer& timer,
PolymerBlackoilState& state);
private:
// Data.
// Parameters for output.
bool output_;
bool output_vtk_;
std::string output_dir_;
int output_interval_;
// Parameters for well control
bool check_well_controls_;
int max_well_control_iterations_;
// Observed objects.
const UnstructuredGrid& grid_;
const BlackoilPropsAdInterface& props_;
const PolymerPropsAd& polymer_props_;
const RockCompressibility* rock_comp_props_;
std::shared_ptr<EclipseState> eclipse_state_;
EclipseWriter& output_writer_;
Opm::DeckConstPtr& deck_;
NewtonIterationBlackoilInterface& linsolver_;
const double* gravity_;
// Solvers
DerivedGeology geo_;
// Misc. data
std::vector<int> allcells_;
};
SimulatorFullyImplicitCompressiblePolymer::
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
{
pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, eclipse_state, output_writer, deck, linsolver, gravity));
}
SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer,
PolymerBlackoilState& state)
{
return pimpl_->run(timer, state);
}
// \TODO: Treat bcs.
SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid,
const DerivedGeology& geo,
const BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
: grid_(grid),
props_(props),
polymer_props_(polymer_props),
rock_comp_props_(rock_comp_props),
eclipse_state_(eclipse_state),
output_writer_(output_writer),
deck_(deck),
linsolver_(linsolver),
gravity_(gravity),
geo_(geo)
{
// For output.
output_ = param.getDefault("output", true);
if (output_) {
output_vtk_ = param.getDefault("output_vtk", true);
output_dir_ = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir_);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
output_interval_ = param.getDefault("output_interval", 1);
}
// Well control related init.
check_well_controls_ = param.getDefault("check_well_controls", false);
max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
// Misc init.
const int num_cells = grid.number_of_cells;
allcells_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells_[cell] = cell;
}
}
SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer,
PolymerBlackoilState& state)
{
WellStateFullyImplicitBlackoil prev_well_state;
// Initialisation.
std::vector<double> porevol;
if (rock_comp_props_ && rock_comp_props_->isActive()) {
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
} else {
computePorevolume(grid_, props_.porosity(), porevol);
}
std::vector<double> initial_porevol = porevol;
std::vector<double> polymer_inflow_c(grid_.number_of_cells);
// Main simulation loop.
Opm::time::StopWatch solver_timer;
double stime = 0.0;
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
std::string tstep_filename = output_dir_ + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
//Main simulation loop.
while (!timer.done()) {
#if 0
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
if (wells_) {
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
}
std::fstream tstep_os;
if (output_) {
std::string filename = output_dir_ + "/step_timing.param";
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
}
#endif
// Report timestep and (optionally) write state to disk.
step_timer.start();
timer.report(std::cout);
WellsManager wells_manager(eclipse_state_,
timer.currentStepNum(),
Opm::UgGridHelpers::numCells(grid_),
Opm::UgGridHelpers::globalCell(grid_),
Opm::UgGridHelpers::cartDims(grid_),
Opm::UgGridHelpers::dimensions(grid_),
Opm::UgGridHelpers::cell2Faces(grid_),
Opm::UgGridHelpers::beginFaceCentroids(grid_),
props_.permeability());
const Wells* wells = wells_manager.c_wells();
WellStateFullyImplicitBlackoil well_state;
well_state.init(wells, state.blackoilState(), prev_well_state);
//Compute polymer inflow.
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, eclipse_state_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
}
if (output_) {
if (timer.currentStepNum() == 0) {
output_writer_.writeInit(timer);
}
output_writer_.writeTimeStep(timer, state.blackoilState(), well_state);
}
// Run solver.
solver_timer.start();
FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_);
solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
// Stop timer and report.
solver_timer.stop();
const double st = solver_timer.secsSinceStart();
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
stime += st;
// Update pore volumes if rock is compressible.
if (rock_comp_props_ && rock_comp_props_->isActive()) {
initial_porevol = porevol;
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
}
/*
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double polyinj = 0;
double polyprod = 0;
Opm::computeInjectedProduced(props_, polymer_props_,
state,
transport_src, polymer_inflow_c, timer.currentStepLength(),
injected, produced,
polyinj, polyprod);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
std::cout.precision(5);
const int width = 18;
std::cout << "\nMass balance report.\n";
std::cout << " Injected reservoir volumes: "
<< std::setw(width) << injected[0]
<< std::setw(width) << injected[1] << std::endl;
std::cout << " Produced reservoir volumes: "
<< std::setw(width) << produced[0]
<< std::setw(width) << produced[1] << std::endl;
std::cout << " Total inj reservoir volumes: "
<< std::setw(width) << tot_injected[0]
<< std::setw(width) << tot_injected[1] << std::endl;
std::cout << " Total prod reservoir volumes: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1] << std::endl;
*/
if (output_) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
step_report.reportParam(tstep_os);
}
++timer;
prev_well_state = well_state;
}
// Write final simulation state.
if (output_) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state);
}
total_timer.stop();
SimulatorReport report;
report.pressure_time = stime;
report.transport_time = 0.0;
report.total_time = total_timer.secsSinceStart();
return report;
}
namespace
{
static void outputStateVtk(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/vtk_files";
boost::filesystem::path fpath(vtkfilename.str());
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
}
static void outputStateMatlab(const UnstructuredGrid& grid,
const Opm::PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["cmax"] = &state.maxconcentration();
dm["concentration"] = &state.concentration();
dm["surfvolume"] = &state.surfacevol();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
#if 0
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
std::ofstream os(fname.c_str());
if (!os) {
OPM_THROW(std::runtime_error, "Failed to open " << fname);
}
watercut.write(os);
}
#endif
}
} // namespace Opm

132
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp
View File

@ -21,86 +21,102 @@
#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#include <memory> #include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <vector> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
struct UnstructuredGrid; #include <opm/autodiff/GeoProps.hpp>
struct Wells; #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <numeric>
#include <fstream>
#include <iostream>
namespace Opm namespace Opm
{ {
namespace parameter { class ParameterGroup; } template <class GridT>
class BlackoilPropsAdInterface; class SimulatorFullyImplicitCompressiblePolymer;
class RockCompressibility;
class DerivedGeology; template <class GridT>
class WellStateFullyImplicitBlackoil; struct SimulatorTraits<SimulatorFullyImplicitCompressiblePolymer<GridT> >
class WellsManager; {
class EclipseWriter; typedef PolymerBlackoilState ReservoirState;
class EclipseState; typedef WellStateFullyImplicitBlackoilPolymer WellState;
class NewtonIterationBlackoilInterface; typedef BlackoilOutputWriter OutputWriter;
class SimulatorTimer; typedef GridT Grid;
class PolymerBlackoilState; typedef FullyImplicitCompressiblePolymerSolver Solver;
class PolymerPropsAd; };
class PolymerInflowInterface;
struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation. /// Class collecting all necessary components for a two-phase simulation.
template <class GridT>
class SimulatorFullyImplicitCompressiblePolymer class SimulatorFullyImplicitCompressiblePolymer
: public SimulatorBase<SimulatorFullyImplicitCompressiblePolymer<GridT> >
{ {
typedef SimulatorFullyImplicitCompressiblePolymer ThisType;
typedef SimulatorBase<ThisType> BaseType;
typedef typename BaseType::Solver Solver;
public: public:
/// Initialise from parameters and objects to observe. /// Initialise from parameters and objects to observe.
/// \param[in] param parameters, this class accepts the following:
/// parameter (default) effect
/// -----------------------------------------------------------
/// output (true) write output to files?
/// output_dir ("output") output directoty
/// output_interval (1) output every nth step
/// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
/// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal)
/// nl_pressure_maxiter (10) max nonlinear iterations in pressure
/// nl_maxiter (30) max nonlinear iterations in transport
/// nl_tolerance (1e-9) transport solver absolute residual tolerance
/// num_transport_substeps (1) number of transport steps per pressure step
/// use_segregation_split (false) solve for gravity segregation (if false,
/// segregation is ignored).
///
/// \param[in] grid grid data structure
/// \param[in] props fluid and rock properties
/// \param[in] polymer_props polymer properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] eclipse_state
/// \param[in] eclipse_writer
/// \param[in] deck
/// \param[in] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param, SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
const UnstructuredGrid& grid, const GridT& grid,
const DerivedGeology& geo, const DerivedGeology& geo,
const BlackoilPropsAdInterface& props, BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props, const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props, const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state, std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& eclipse_writer, BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck, Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver, NewtonIterationBlackoilInterface& linsolver,
const double* gravity); const double* gravity);
/// Run the simulation. Solver* createSolver(const Wells* wells);
/// This will run succesive timesteps until timer.done() is true. It will
/// modify the reservoir and well states.
/// \param[in,out] timer governs the requested reporting timesteps
/// \param[in,out] state state of reservoir: pressure, fluxes
/// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer,
PolymerBlackoilState& state);
private: void handleAdditionalWellInflow(SimulatorTimer& timer,
class Impl; WellsManager& wells_manager,
// Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. typename BaseType::WellState& well_state,
std::shared_ptr<Impl> pimpl_; const Wells* wells);
private:
Opm::DeckConstPtr deck_;
const PolymerPropsAd& polymer_props_;
}; };
} // namespace Opm } // namespace Opm
#include "SimulatorFullyImplicitCompressiblePolymer_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED

101
opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp Normal file
View File

@ -0,0 +1,101 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_IMPL_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_IMPL_HEADER_INCLUDED
namespace Opm
{
/// Class collecting all necessary components for a two-phase simulation.
template <class GridT>
SimulatorFullyImplicitCompressiblePolymer<GridT>::
SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
const GridT& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
std::shared_ptr<EclipseState> eclipse_state,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity)
: BaseType(param,
grid,
geo,
props,
rock_comp_props,
linsolver,
gravity,
/*disgas=*/false,
/*vapoil=*/false,
eclipse_state,
output_writer,
/*threshold_pressures_by_face=*/std::vector<double>())
, deck_(deck)
, polymer_props_(polymer_props)
{
}
template <class GridT>
auto SimulatorFullyImplicitCompressiblePolymer<GridT>::
createSolver(const Wells* wells)
-> Solver*
{
return new Solver(BaseType::grid_,
BaseType::props_,
BaseType::geo_,
BaseType::rock_comp_props_,
polymer_props_,
*wells,
BaseType::solver_);
}
template <class GridT>
void SimulatorFullyImplicitCompressiblePolymer<GridT>::
handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,
typename BaseType::WellState& well_state,
const Wells* wells)
{
// compute polymer inflow
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
well_state.polymerInflow() = polymer_inflow_c;
}
} // namespace Opm
#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED