OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-16 07:51:56 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #96 from qilicun/updates

Browse Source 
Updates FIBOPOLYMER simulator based on flow simulator.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-03-24 20:29:46 +01:00
commit 62e8146424
7 changed files with 985 additions and 972 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

118
examples/sim_poly_fibo_ad.cpp → examples/flow_polymer.cpp
View File

@ -1,5 +1,5 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
@ -34,7 +34,6 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
@ -53,8 +52,13 @@
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.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/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
@ -87,7 +91,7 @@ try
{
using namespace Opm;
std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n";
std::cout << "\n================ Test program for fully implicit three-phase black-oil-polymer flow ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
@ -99,7 +103,7 @@ try
}
std::shared_ptr<GridManager> grid;
std::shared_ptr<BlackoilPropertiesInterface> props;
std::shared_ptr<BlackoilPropsAdInterface> new_props;
std::shared_ptr<BlackoilPropsAdFromDeck> new_props;
std::shared_ptr<RockCompressibility> rock_comp;
PolymerBlackoilState state;
// bool check_well_controls = false;
@ -107,23 +111,57 @@ try
double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr newParser(new Opm::Parser() );
Opm::DeckConstPtr deck = newParser->parseFile(deck_filename);
std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck));
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
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 (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
// 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
std::vector<double> porv;
if (eclipseState->hasDoubleGridProperty("PORV")) {
porv = eclipseState->getDoubleGridProperty("PORV")->getData();
}
std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
Opm::EclipseWriter outputWriter(param,
eclipseState,
pu,
cGrid.number_of_cells,
cGrid.global_cell);
Opm::BlackoilOutputWriter outputWriter(cGrid,
param,
eclipseState,
pu );
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
@ -142,8 +180,8 @@ try
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
initBlackoilSurfvol(*grid->c_grid(), *props, state);
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state.blackoilState());
initBlackoilSurfvol(*grid->c_grid(), *props, state.blackoilState());
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int np = props->numPhases();
@ -159,10 +197,19 @@ try
initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState());
state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
} else {
state.init(*grid->c_grid(), props->numPhases());
initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState());
}
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
if (deck->hasKeyword("SWATINIT")) {
const int nc = grid->c_grid()->number_of_cells;
std::vector<int> cells(nc);
for (int c = 0; c < nc; ++c) { cells[c] = c; }
std::vector<double> pc = state.saturation();
props->capPress(nc, state.saturation().data(), cells.data(), pc.data(),NULL);
new_props->setSwatInitScaling(state.saturation(),pc);
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
@ -174,24 +221,6 @@ try
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
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 (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
SimulatorTimer simtimer;
@ -210,13 +239,13 @@ try
if (use_wpolymer) {
OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario.");
}
}
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
SimulatorReport fullReport;
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
}
bool use_local_perm = param.getDefault("use_local_perm", true);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, use_local_perm, grav);
std::vector<double> threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid());
SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid> simulator(param,
*grid->c_grid(),
geology,
@ -233,8 +262,11 @@ try
deck,
threshold_pressures);
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
SimulatorReport fullReport = simulator.run(simtimer, state);
fullReport = simulator.run(simtimer, state);
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.report(std::cout);

179
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp
View File

@ -1,7 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
@ -22,6 +21,8 @@
#ifndef OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
#define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
#include <cassert>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
@ -30,6 +31,8 @@
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <array>
struct UnstructuredGrid;
struct Wells;
@ -56,29 +59,55 @@ namespace Opm {
class FullyImplicitBlackoilPolymerSolver
{
public:
// the Newton relaxation type
enum RelaxType { DAMPEN, SOR };
// class holding the solver parameters
struct SolverParameter
{
double dp_max_rel_;
double ds_max_;
double dr_max_rel_;
enum RelaxType relax_type_;
double relax_max_;
double relax_increment_;
double relax_rel_tol_;
double max_residual_allowed_;
double tolerance_mb_;
double tolerance_cnv_;
double tolerance_wells_;
int max_iter_;
SolverParameter( const parameter::ParameterGroup& param );
SolverParameter();
void reset();
};
/// \brief The type of the grid that we use.
typedef T Grid;
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param,
FullyImplicitBlackoilPolymerSolver(const SolverParameter& param,
const Grid& grid ,
const BlackoilPropsAdInterface& fluid,
const DerivedGeology& geo ,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const Wells& wells,
const Wells* wells,
const NewtonIterationBlackoilInterface& linsolver,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer);
const bool has_polymer,
const bool terminal_output);
/// \brief Set threshold pressures that prevent or reduce flow.
/// This prevents flow across faces if the potential
@ -99,12 +128,16 @@ namespace Opm {
/// \param[in] dt time step size
/// \param[in] state reservoir state
/// \param[in] wstate well state
void
/// \return number of linear iterations used
int
step(const double dt ,
PolymerBlackoilState& state ,
WellStateFullyImplicitBlackoil& wstate,
const std::vector<double>& polymer_inflow);
unsigned int newtonIterations () const { return newtonIterations_; }
unsigned int linearIterations () const { return linearIterations_; }
private:
// Types and enums
typedef AutoDiffBlock<double> ADB;
@ -133,21 +166,23 @@ namespace Opm {
ADB rv;
ADB concentration;
ADB qs;
ADB bhp;
ADB bhp;
// Below are quantities stored in the state for optimization purposes.
std::vector<ADB> canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active.
};
struct WellOps {
WellOps(const Wells& wells);
WellOps(const Wells* wells);
M w2p; // well -> perf (scatter)
M p2w; // perf -> well (gather)
};
enum { Water = BlackoilPropsAdInterface::Water,
Oil = BlackoilPropsAdInterface::Oil ,
Gas = BlackoilPropsAdInterface::Gas };
enum { Water = BlackoilPropsAdInterface::Water,
Oil = BlackoilPropsAdInterface::Oil ,
Gas = BlackoilPropsAdInterface::Gas ,
MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases
};
// the Newton relaxation type
enum RelaxType { DAMPEN, SOR };
enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
// Member data
@ -155,8 +190,8 @@ namespace Opm {
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
const PolymerPropsAd& polymer_props_ad_;
const Wells& wells_;
const PolymerPropsAd& polymer_props_ad_;
const Wells* wells_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
const std::vector<bool> active_;
@ -170,14 +205,8 @@ namespace Opm {
const bool has_vapoil_;
const bool has_polymer_;
const int poly_pos_;
double dp_max_rel_;
double ds_max_;
double drs_max_rel_;
enum RelaxType relax_type_;
double relax_max_;
double relax_increment_;
double relax_rel_tol_;
int max_iter_;
SolverParameter param_;
bool use_threshold_pressure_;
V threshold_pressures_by_interior_face_;
@ -187,16 +216,30 @@ namespace Opm {
LinearisedBlackoilResidual residual_;
/// \brief Whether we print something to std::cout
bool terminal_output_;
unsigned int newtonIterations_;
unsigned int linearIterations_;
std::vector<int> primalVariable_;
// Private methods.
// return true if wells are available
bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; }
// return wells object
const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
SolutionState
constantState(const PolymerBlackoilState& x,
const WellStateFullyImplicitBlackoil& xw);
const WellStateFullyImplicitBlackoil& xw) const;
void
makeConstantState(SolutionState& state) const;
SolutionState
variableState(const PolymerBlackoilState& x,
const WellStateFullyImplicitBlackoil& xw);
const WellStateFullyImplicitBlackoil& xw) const;
void
computeAccum(const SolutionState& state,
@ -223,6 +266,7 @@ namespace Opm {
void
assemble(const V& dtpv,
const PolymerBlackoilState& x,
const bool initial_assembly,
WellStateFullyImplicitBlackoil& xw,
const std::vector<double>& polymer_inflow);
@ -235,6 +279,18 @@ namespace Opm {
std::vector<ADB>
computePressures(const SolutionState& state) const;
std::vector<ADB>
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg) const;
V
computeGasPressure(const V& po,
const V& sw,
const V& so,
const V& sg) const;
std::vector<ADB>
computeRelPerm(const SolutionState& state) const;
@ -244,35 +300,26 @@ namespace Opm {
const std::vector<int>& well_cells) const;
void
computeMassFlux(const int actph ,
const V& transi,
const ADB& kr ,
const ADB& p ,
const SolutionState& state );
void
computeMassFlux(const V& trans,
const std::vector<ADB>& kr,
computeMassFlux(const V& transi,
const std::vector<ADB>& kr ,
const std::vector<ADB>& phasePressure,
const SolutionState& state);
const SolutionState& state );
void
computeCmax(PolymerBlackoilState& state);
ADB
computeMc(const SolutionState& state) const;
ADB
rockPorosity(const ADB& p) const;
ADB
rockPermeability(const ADB& p) const;
void applyThresholdPressures(ADB& dp);
double
residualNorm() const;
std::vector<double> residuals() const;
/// \brief Compute the residual norms of the mass balance for each phase,
/// the well flux, and the well equation.
/// \return a vector that contains for each phase the norm of the mass balance
/// and afterwards the norm of the residual of the well flux and the well equation.
std::vector<double> computeResidualNorms() const;
ADB
fluidViscosity(const int phase,
@ -321,7 +368,6 @@ namespace Opm {
const ADB& so,
const std::vector<int>& cells) const;
ADB
poroMult(const ADB& p) const;
@ -350,7 +396,35 @@ namespace Opm {
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
bool getConvergence(const double dt, const int it);
bool getConvergence(const double dt, const int iteration);
/// \brief Compute the reduction within the convergence check.
/// \param[in] B A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. B.col(i) contains the values
/// for phase i.
/// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. tempV.col(i) contains the
/// values
/// for phase i.
/// \param[in] R A matrix with MaxNumPhases columns and the same number rows
/// as the number of cells of the grid. B.col(i) contains the values
/// for phase i.
/// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum
/// of R for the phase i.
/// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the
/// maximum of tempV for the phase i.
/// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
/// of B for the phase i.
/// \param[in] nc The number of cells of the local grid.
/// \return The total pore volume over all cells.
double
convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases+1>& R,
std::array<double,MaxNumPhases+1>& R_sum,
std::array<double,MaxNumPhases+1>& maxCoeff,
std::array<double,MaxNumPhases+1>& B_avg,
int nc) const;
void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol,
@ -358,14 +432,15 @@ namespace Opm {
void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
double dpMaxRel() const { return dp_max_rel_; }
double dsMax() const { return ds_max_; }
double drsMaxRel() const { return drs_max_rel_; }
enum RelaxType relaxType() const { return relax_type_; }
double relaxMax() const { return relax_max_; };
double relaxIncrement() const { return relax_increment_; };
double relaxRelTol() const { return relax_rel_tol_; };
double maxIter() const { return max_iter_; }
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }
double drMaxRel() const { return param_.dr_max_rel_; }
enum RelaxType relaxType() const { return param_.relax_type_; }
double relaxMax() const { return param_.relax_max_; };
double relaxIncrement() const { return param_.relax_increment_; };
double relaxRelTol() const { return param_.relax_rel_tol_; };
double maxIter() const { return param_.max_iter_; }
double maxResidualAllowed() const { return param_.max_residual_allowed_; }
};
} // namespace Opm

1049
opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp
View File

File diff suppressed because it is too large Load Diff

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

@ -1,5 +1,5 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
@ -26,6 +26,7 @@
struct UnstructuredGrid;
struct Wells;
struct FlowBoundaryConditions;
namespace Opm
{
@ -37,9 +38,8 @@ namespace Opm
class SimulatorTimer;
class PolymerBlackoilState;
class WellStateFullyImplicitBlackoil;
class WellsManager;
class EclipseState;
class EclipseWriter;
class BlackoilOutputWriter;
class PolymerPropsAd;
class PolymerInflowInterface;
struct SimulatorReport;
@ -70,7 +70,6 @@ namespace Opm
/// \param[in] grid grid data structure
/// \param[in] geo derived geological properties
/// \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] linsolver linear solver
/// \param[in] gravity if non-null, gravity vector
@ -91,7 +90,7 @@ namespace Opm
const bool vapoil,
const bool polymer,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face);

211
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp
View File

@ -1,211 +0,0 @@
#include "config.h"
#include "SimulatorFullyImplicitBlackoilPolymerOutput.hpp"
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <boost/filesystem.hpp>
#ifdef HAVE_DUNE_CORNERPOINT
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
#endif
namespace Opm
{
void outputStateVtk(const UnstructuredGrid& grid,
const 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["concentration"] = &state.concentration();
dm["maxconcentration"] = &state.maxconcentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
AutoDiffGrid::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(grid),
state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
}
void outputStateMatlab(const UnstructuredGrid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["surfvolume"] = &state.surfacevol();
dm["rs"] = &state.gasoilratio();
dm["rv"] = &state.rv();
dm["concentration"] = &state.concentration();
dm["maxconcentration"] = &state.maxconcentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
UgGridHelpers::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(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"));
}
}
void outputWellStateMatlab(const Opm::WellState& well_state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["bhp"] = &well_state.bhp();
dm["wellrates"] = &well_state.wellRates();
// 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
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);
}
void outputWellReport(const Opm::WellReport& wellreport,
const std::string& output_dir)
{
// Write well report.
std::string fname = output_dir + "/wellreport.txt";
std::ofstream os(fname.c_str());
if (!os) {
OPM_THROW(std::runtime_error, "Failed to open " << fname);
}
wellreport.write(os);
}
#endif
#ifdef HAVE_DUNE_CORNERPOINT
void outputStateVtk(const Dune::CpGrid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
std::ostringstream vtkpath;
vtkpath << output_dir << "/vtk_files";
vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step;
boost::filesystem::path fpath(vtkpath.str());
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafGridView(), Dune::VTK::nonconforming);
#else
Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
#endif
writer.addCellData(state.saturation(), "saturation", state.numPhases());
writer.addCellData(state.pressure(), "pressure", 1);
writer.addCellData(state.concentration(), "concentration", 1);
writer.addCellData(state.maxconcentration(), "maxconcentration", 1);
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
AutoDiffGrid::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(grid),
state.faceflux(), cell_velocity);
writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension);
writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii);
}
#endif
}

96
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp
View File

@ -1,96 +0,0 @@
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <boost/filesystem.hpp>
#ifdef HAVE_DUNE_CORNERPOINT
#include <dune/grid/CpGrid.hpp>
#endif
namespace Opm
{
void outputStateVtk(const UnstructuredGrid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir);
void outputStateMatlab(const UnstructuredGrid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir);
void outputWellStateMatlab(const Opm::WellState& well_state,
const int step,
const std::string& output_dir);
#ifdef HAVE_DUNE_CORNERPOINT
void outputStateVtk(const Dune::CpGrid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir);
#endif
template<class Grid>
void outputStateMatlab(const Grid& grid,
const PolymerBlackoilState& state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["surfvolume"] = &state.surfacevol();
dm["rs"] = &state.gasoilratio();
dm["rv"] = &state.rv();
dm["concentration"] = &state.concentration();
dm["maxconcentration"] = &state.maxconcentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
AutoDiffGrid::numFaces(grid),
AutoDiffGrid::beginFaceCentroids(grid),
UgGridHelpers::faceCells(grid),
AutoDiffGrid::beginCellCentroids(grid),
AutoDiffGrid::beginCellVolumes(grid),
AutoDiffGrid::dimensions(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"));
}
}
}
#endif

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

@ -1,5 +1,6 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 IRIS AS
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
@ -18,11 +19,13 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp>
//#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
@ -36,9 +39,9 @@
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#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>
@ -46,6 +49,7 @@
#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>
@ -80,7 +84,7 @@ namespace Opm
const Grid& grid,
const DerivedGeology& geo,
BlackoilPropsAdInterface& props,
const PolymerPropsAd& polymer_props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
@ -88,7 +92,7 @@ namespace Opm
bool has_vapoil,
bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face);
@ -103,11 +107,6 @@ namespace Opm
const parameter::ParameterGroup param_;
// Parameters for output.
bool output_;
bool output_vtk_;
std::string output_dir_;
int output_interval_;
// Observed objects.
const Grid& grid_;
BlackoilPropsAdInterface& props_;
@ -122,10 +121,11 @@ namespace Opm
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
EclipseWriter& output_writer_;
BlackoilOutputWriter& output_writer_;
Opm::DeckConstPtr& deck_;
RateConverterType rateConverter_;
// Threshold pressures.
@ -134,7 +134,7 @@ namespace Opm
void
computeRESV(const std::size_t step,
const Wells* wells,
const PolymerBlackoilState& x,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw);
};
@ -154,13 +154,13 @@ namespace Opm
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
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));
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));
}
@ -175,64 +175,6 @@ namespace Opm
}
static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
const int step,
const std::string& output_dir)
{
Opm::DataMap dm;
dm["bhp"] = &well_state.bhp();
dm["wellrates"] = &well_state.wellRates();
// 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);
}
static void outputWellReport(const Opm::WellReport& wellreport,
const std::string& output_dir)
{
// Write well report.
std::string fname = output_dir + "/wellreport.txt";
std::ofstream os(fname.c_str());
if (!os) {
OPM_THROW(std::runtime_error, "Failed to open " << fname);
}
wellreport.write(os);
}
#endif
// \TODO: Treat bcs.
template<class T>
SimulatorFullyImplicitBlackoilPolymer<T>::Impl::Impl(const parameter::ParameterGroup& param,
@ -247,7 +189,7 @@ namespace Opm
const bool has_vapoil,
const bool has_polymer,
std::shared_ptr<EclipseState> eclipse_state,
EclipseWriter& output_writer,
BlackoilOutputWriter& output_writer,
Opm::DeckConstPtr& deck,
const std::vector<double>& threshold_pressures_by_face)
: param_(param),
@ -261,34 +203,30 @@ namespace Opm
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)
{
// 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);
}
// 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
}
@ -306,14 +244,40 @@ namespace Opm
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
std::string tstep_filename = output_dir_ + "/step_timing.txt";
std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
std::ofstream tstep_os(tstep_filename.c_str());
typename FullyImplicitBlackoilPolymerSolver<T>::SolverParameter solverParam( param_ );
//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();
timer.report(std::cout);
if ( terminal_output_ )
{
timer.report(std::cout);
}
// Create wells and well state.
WellsManager wells_manager(eclipse_state_,
@ -328,6 +292,7 @@ namespace Opm
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")) {
@ -344,42 +309,54 @@ namespace Opm
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
// Output state at start of time step.
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
}
if (output_) {
if (timer.currentStepNum() == 0) {
output_writer_.writeInit(timer);
}
output_writer_.writeTimeStep(timer, state.blackoilState(), well_state);
}
// 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, well_state);
computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state);
// Run a single step of the solver.
// Run a multiple steps of the solver depending on the time step control.
solver_timer.start();
FullyImplicitBlackoilPolymerSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells, solver_, has_disgas_, has_vapoil_, has_polymer_);
FullyImplicitBlackoilPolymerSolver<T> solver(solverParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_);
if (!threshold_pressures_by_face_.empty()) {
solver.setThresholdPressures(threshold_pressures_by_face_);
}
// 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, polymer_inflow_c);
// }
// 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();
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
if ( terminal_output_ )
{
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
}
stime += st;
if (output_) {
if ( output_writer_.output() ) {
SimulatorReport step_report;
step_report.pressure_time = st;
step_report.total_time = step_timer.secsSinceStart();
@ -392,14 +369,7 @@ namespace Opm
}
// Write final simulation state.
if (output_) {
if (output_vtk_) {
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_);
output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state);
}
output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state );
// Stop timer and create timing report
total_timer.stop();
@ -407,6 +377,8 @@ namespace Opm
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;
}
@ -444,17 +416,16 @@ namespace Opm
}
inline bool
is_resv_prod(const Wells& wells,
const int w)
is_resv(const Wells& wells,
const int w)
{
return ((wells.type[w] == PRODUCER) &&
(0 <= resv_control(wells.ctrls[w])));
return (0 <= resv_control(wells.ctrls[w]));
}
inline bool
is_resv_prod(const WellMap& wmap,
const std::string& name,
const std::size_t step)
is_resv(const WellMap& wmap,
const std::string& name,
const std::size_t step)
{
bool match = false;
@ -465,29 +436,34 @@ namespace Opm
match = (wp->isProducer(step) &&
wp->getProductionProperties(step)
.hasProductionControl(WellProducer::RESV));
.hasProductionControl(WellProducer::RESV))
|| (wp->isInjector(step) &&
wp->getInjectionProperties(step)
.hasInjectionControl(WellInjector::RESV));
}
return match;
}
inline std::vector<int>
resvProducers(const Wells& wells,
const std::size_t step,
const WellMap& wmap)
resvWells(const Wells* wells,
const std::size_t step,
const WellMap& wmap)
{
std::vector<int> resv_prod;
for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) {
if (is_resv_prod(wells, w) ||
((wells.name[w] != 0) &&
is_resv_prod(wmap, wells.name[w], step)))
{
resv_prod.push_back(w);
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_prod;
return resv_wells;
}
inline void
@ -527,7 +503,7 @@ namespace Opm
SimulatorFullyImplicitBlackoilPolymer<T>::
Impl::computeRESV(const std::size_t step,
const Wells* wells,
const PolymerBlackoilState& x,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw)
{
typedef SimFIBODetails::WellMap WellMap;
@ -535,24 +511,24 @@ namespace Opm
const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
const std::vector<int>& resv_prod =
SimFIBODetails::resvProducers(*wells, step, wmap);
const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
if (! resv_prod.empty()) {
if (! resv_wells.empty()) {
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x.blackoilState());
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_prod.begin(), e = resv_prod.end();
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
{
@ -561,11 +537,17 @@ namespace Opm
if (0 <= rctrl) {
const std::vector<double>::size_type off = (*rp) * np;
// 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>());
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);
@ -576,7 +558,7 @@ namespace Opm
// RESV control, WCONHIST wells. A bit of duplicate
// work, regrettably.
if (wells->name[*rp] != 0) {
if (is_producer && wells->name[*rp] != 0) {
WellMap::const_iterator i = wmap.find(wells->name[*rp]);
if (i != wmap.end()) {
@ -603,11 +585,18 @@ namespace Opm
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
const int ok =
const int ok_resv =
well_controls_add_new(RESERVOIR_RATE, target,
& distr[0], ctrl);
if (ok != 0) {
// 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);
}