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 remote-tracking branch 'upstream/master' into master-refactor-for-cpgrid-support

Browse Source 
Manually resolved conflicts:
	opm/core/io/eclipse/EclipseWriter.cpp
	opm/core/io/eclipse/EclipseWriter.hpp
	opm/core/props/BlackoilPropertiesFromDeck.cpp
	opm/core/simulator/initState_impl.hpp
	opm/core/wells/WellsManager.cpp
	opm/core/wells/WellsManager.hpp
This commit is contained in:
Markus Blatt 2014-03-27 15:27:11 +01:00
commit c282949400
11 changed files with 396 additions and 838 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -30,6 +30,7 @@
// TODO: clean up includes.
#include <dune/common/deprecated.hh>
#include <dune/common/version.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
@ -39,6 +40,10 @@
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
#include <dune/istl/paamg/fastamg.hh>
#endif
#include <stdexcept>
#include <iostream>
@ -61,7 +66,7 @@ namespace Opm
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
@ -173,7 +178,7 @@ namespace Opm
linsolver_prolongate_factor_, linsolver_smooth_steps_);
break;
case KAMG:
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_);
#else
@ -181,7 +186,7 @@ namespace Opm
#endif
break;
case FastAMG:
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_);
#else
@ -311,7 +316,7 @@ namespace Opm
}
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps)

2
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@ -240,7 +240,7 @@ namespace Opm
"SaturationPropsFromDeck::init() -- ENDSCALE: "
"Currently only 'NODIR' accepted.");
}
if (endscale.isReversible()) {
if (!endscale.isReversible()) {
OPM_THROW(std::runtime_error,
"SaturationPropsFromDeck::init() -- ENDSCALE: "
"Currently only 'REVERS' accepted.");

58
opm/core/simulator/WellState.hpp
View File

@ -23,6 +23,7 @@
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <vector>
#include <cassert>
namespace Opm
{
@ -45,27 +46,56 @@ namespace Opm
bhp_.resize(nw);
wellrates_.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
// Initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a little
// above or below (depending on if the well is an
// injector or producer) pressure in first perforation
// cell.
if (well_controls_well_is_shut(ctrl) || (well_controls_get_current_type(ctrl) != BHP)) {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*state.pressure()[first_cell];
} else {
bhp_[w] = well_controls_get_current_target( ctrl );
if (well_controls_well_is_shut(ctrl)) {
// Shut well:
// 1. Assign zero well rates.
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = 0.0;
}
// Initialize well rates to match controls if type is SURFACE_RATE
if (well_controls_well_is_open( ctrl ) || (well_controls_get_current_type(ctrl) == SURFACE_RATE)) {
// 2. Assign bhp equal to bhp control, if
// applicable, otherwise assign equal to
// first perforation cell pressure.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[first_cell];
}
} else {
// Open well:
// 1. Initialize well rates to match controls
// if type is SURFACE_RATE. Otherwise, we
// cannot set the correct value here, so we
// assign a small rate with the correct
// sign so that any logic depending on that
// sign will work as expected.
if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
const double rate_target = well_controls_get_current_target(ctrl);
const double * distr = well_controls_get_current_distr( ctrl );
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = rate_target * distr[p];
}
} else {
const double small_rate = 1e-14;
const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0;
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = small_rate * sign;
}
}
// 2. Initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*state.pressure()[first_cell];
}
}
}
// The perforation rates and perforation pressures are

9
opm/core/simulator/initState_impl.hpp
View File

@ -584,6 +584,11 @@ namespace Opm
"found " << pu.num_phases << " phases in deck.");
}
state.init(number_of_cells, number_of_faces, num_phases);
if (newParserDeck->hasKeyword("EQUIL") && newParserDeck->hasKeyword("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
if (newParserDeck->hasKeyword("EQUIL")) {
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
@ -703,6 +708,10 @@ namespace Opm
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
state.init(number_of_cells, number_of_faces, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {

22
opm/core/wells/WellsGroup.cpp
View File

@ -1171,19 +1171,21 @@ namespace Opm
InjectionSpecification injection_specification;
ProductionSpecification production_specification;
if (well->isInjector(timeStep)) {
injection_specification.BHP_limit_ = well->getBHPLimit(timeStep);
injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep)));
injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep)));
injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep);
injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep);
const WellInjectionProperties& properties = well->getInjectionProperties(timeStep);
injection_specification.BHP_limit_ = properties.BHPLimit;
injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(properties.injectorType));
injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(properties.controlMode));
injection_specification.surface_flow_max_rate_ = properties.surfaceInjectionRate;
injection_specification.reservoir_flow_max_rate_ = properties.reservoirInjectionRate;
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
else if (well->isProducer(timeStep)) {
production_specification.BHP_limit_ = well->getBHPLimit(timeStep);
production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep);
production_specification.oil_max_rate_ = well->getOilRate(timeStep);
production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep)));
production_specification.water_max_rate_ = well->getWaterRate(timeStep);
const WellProductionProperties& properties = well->getProductionProperties(timeStep);
production_specification.BHP_limit_ = properties.BHPLimit;
production_specification.reservoir_flow_max_rate_ = properties.ResVRate;
production_specification.oil_max_rate_ = properties.OilRate;
production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(properties.controlMode));
production_specification.water_max_rate_ = properties.WaterRate;
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
std::shared_ptr<WellsGroupInterface> wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage));

180
opm/core/wells/WellsManager.cpp
View File

@ -231,113 +231,20 @@ namespace Opm
{
}
/// Construct from existing wells object.
WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
: w_(clone_wells(W)), checkCellExistence_(checkCellExistence)
{
}
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence)
: w_(0), checkCellExistence_(checkCellExistence)
const double* permeability)
: w_(0)
{
if (UgGridHelpers::dimensions(grid) != 3) {
OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
}
if (eclipseState->getSchedule()->numWells() == 0) {
OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
return;
}
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(UgGridHelpers::globalCell(grid),
UgGridHelpers::numCells(grid), cartesian_to_compressed);
// Obtain phase usage data.
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
// These data structures will be filled in this constructor,
// then used to initialize the Wells struct.
std::vector<std::string> well_names;
std::vector<WellData> well_data;
// For easy lookup:
std::map<std::string, int> well_names_to_index;
ScheduleConstPtr schedule = eclipseState->getSchedule();
std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
well_names.reserve(wells.size());
well_data.reserve(wells.size());
createWellsFromSpecs(wells, timeStep, UgGridHelpers::cell2Faces(grid),
UgGridHelpers::cartDims(grid),
UgGridHelpers::beginFaceCentroids(grid),
UgGridHelpers::beginCellCentroids(grid),
UgGridHelpers::dimensions(grid),
well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
setupWellControls(wells, timeStep, well_names, pu);
{
GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
well_collection_.addField(fieldGroup, timeStep, pu);
addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
}
for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
well_collection_.addWell((*wellIter), timeStep, pu);
}
well_collection_.setWellsPointer(w_);
well_collection_.applyGroupControls();
setupGuideRates(wells, timeStep, well_data, well_names_to_index);
// Debug output.
#define EXTRA_OUTPUT
#ifdef EXTRA_OUTPUT
/*
std::cout << "\t WELL DATA" << std::endl;
for(int i = 0; i< num_wells; ++i) {
std::cout << i << ": " << well_data[i].type << " "
<< well_data[i].control << " " << well_data[i].target
<< std::endl;
}
std::cout << "\n\t PERF DATA" << std::endl;
for(int i=0; i< int(wellperf_data.size()); ++i) {
for(int j=0; j< int(wellperf_data[i].size()); ++j) {
std::cout << i << ": " << wellperf_data[i][j].cell << " "
<< wellperf_data[i][j].well_index << std::endl;
}
}
*/
#endif
}
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence)
: w_(0), checkCellExistence_(checkCellExistence)
{
init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.dimensions,
grid.cell_centroids, UgGridHelpers::cell2Faces(grid), grid.face_centroids,
init(eclipseState, timeStep, UgGridHelpers::numCells(grid),
UgGridHelpers::globalCell(grid), UgGridHelpers::cartDims(grid),
UgGridHelpers::dimensions(grid), UgGridHelpers::beginCellCentroids(grid),
UgGridHelpers::cell2Faces(grid), UgGridHelpers::beginFaceCentroids(grid),
permeability);
}
}
/// Destructor.
WellsManager::~WellsManager()
@ -423,14 +330,15 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
}
if (well->isInjector(timeStep)) {
clear_well_controls(well_index, w_);
const WellInjectionProperties& injectionProperties = well->getInjectionProperties(timeStep);
int ok = 1;
int control_pos[5] = { -1, -1, -1, -1, -1 };
if (well->hasInjectionControl(timeStep , WellInjector::RATE)) {
clear_well_controls(well_index, w_);
if (injectionProperties.hasInjectionControl(WellInjector::RATE)) {
control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]);
double distr[3] = { 0.0, 0.0, 0.0 };
WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
if (injectorType == WellInjector::TypeEnum::WATER) {
distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
@ -441,16 +349,16 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
}
ok = append_well_controls(SURFACE_RATE,
well->getSurfaceInjectionRate( timeStep ) ,
injectionProperties.surfaceInjectionRate,
distr,
well_index,
w_);
}
if (ok && well->hasInjectionControl(timeStep , WellInjector::RESV)) {
if (ok && injectionProperties.hasInjectionControl(WellInjector::RESV)) {
control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
double distr[3] = { 0.0, 0.0, 0.0 };
WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep);
WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
if (injectorType == WellInjector::TypeEnum::WATER) {
distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
@ -461,22 +369,23 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
}
ok = append_well_controls(RESERVOIR_RATE,
well->getReservoirInjectionRate( timeStep ),
injectionProperties.reservoirInjectionRate,
distr,
well_index,
w_);
}
if (ok && well->hasInjectionControl(timeStep , WellInjector::BHP)) {
if (ok && injectionProperties.hasInjectionControl(WellInjector::BHP)) {
control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
ok = append_well_controls(BHP,
well->getBHPLimit(timeStep),
injectionProperties.BHPLimit,
NULL,
well_index,
w_);
}
if (ok && well->hasInjectionControl(timeStep , WellInjector::THP)) {
if (ok && injectionProperties.hasInjectionControl(WellInjector::THP)) {
OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
}
@ -487,7 +396,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
{
WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( well->getInjectorControlMode(timeStep) );
WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( injectionProperties.controlMode );
int cpos = control_pos[mode];
if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) {
OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]);
@ -495,7 +404,7 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
// We need to check if the well is shut or not
if (well->getStatus( timeStep ) == WellCommon::SHUT) {
cpos = ~cpos;
well_controls_shut_well( w_->ctrls[well_index] );
}
set_current_control(well_index, cpos, w_);
}
@ -503,34 +412,39 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
// Set well component fraction.
double cf[3] = { 0.0, 0.0, 0.0 };
if (well->getInjectorType(timeStep) == WellInjector::WATER) {
{
WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
if (injectorType == WellInjector::WATER) {
if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
} else if (well->getInjectorType(timeStep) == WellInjector::OIL) {
} else if (injectorType == WellInjector::OIL) {
if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (well->getInjectorType(timeStep) == WellInjector::GAS) {
} else if (injectorType == WellInjector::GAS) {
if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
}
std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
}
}
if (well->isProducer(timeStep)) {
// Add all controls that are present in well.
// First we must clear existing controls, in case the
// current WCONPROD line is modifying earlier controls.
clear_well_controls(well_index, w_);
const WellProductionProperties& productionProperties = well->getProductionProperties(timeStep);
int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
int ok = 1;
if (ok && well->hasProductionControl(timeStep , WellProducer::ORAT)) {
clear_well_controls(well_index, w_);
if (ok && productionProperties.hasProductionControl(WellProducer::ORAT)) {
if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
}
@ -539,13 +453,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
double distr[3] = { 0.0, 0.0, 0.0 };
distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
ok = append_well_controls(SURFACE_RATE,
-well->getOilRate( timeStep ),
-productionProperties.OilRate,
distr,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::WRAT)) {
if (ok && productionProperties.hasProductionControl(WellProducer::WRAT)) {
if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
}
@ -553,13 +467,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
double distr[3] = { 0.0, 0.0, 0.0 };
distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
ok = append_well_controls(SURFACE_RATE,
-well->getWaterRate(timeStep),
-productionProperties.WaterRate,
distr,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::GRAT)) {
if (ok && productionProperties.hasProductionControl(WellProducer::GRAT)) {
if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
}
@ -567,13 +481,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
double distr[3] = { 0.0, 0.0, 0.0 };
distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
ok = append_well_controls(SURFACE_RATE,
-well->getGasRate( timeStep ),
-productionProperties.GasRate,
distr,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::LRAT)) {
if (ok && productionProperties.hasProductionControl(WellProducer::LRAT)) {
if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified.");
}
@ -585,32 +499,32 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
ok = append_well_controls(SURFACE_RATE,
-well->getLiquidRate(timeStep),
-productionProperties.LiquidRate ,
distr,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::RESV)) {
if (ok && productionProperties.hasProductionControl(WellProducer::RESV)) {
control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
double distr[3] = { 1.0, 1.0, 1.0 };
ok = append_well_controls(RESERVOIR_RATE,
-well->getResVRate(timeStep),
-productionProperties.ResVRate ,
distr,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::BHP)) {
if (ok && productionProperties.hasProductionControl(WellProducer::BHP)) {
control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
ok = append_well_controls(BHP,
well->getBHPLimit( timeStep ) ,
productionProperties.BHPLimit ,
NULL,
well_index,
w_);
}
if (ok && well->hasProductionControl(timeStep , WellProducer::THP)) {
if (ok && productionProperties.hasProductionControl(WellProducer::THP)) {
OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]);
}
@ -618,14 +532,16 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]);
}
WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(well->getProducerControlMode(timeStep));
WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(productionProperties.controlMode);
int cpos = control_pos[mode];
if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
}
// If it's shut, we complement the cpos
if (well->getStatus(timeStep) == WellCommon::SHUT) {
cpos = ~cpos; // So we can easily retrieve the cpos later
well_controls_shut_well( w_->ctrls[well_index] );
} else if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
}
set_current_control(well_index, cpos, w_);
}

30
opm/core/wells/WellsManager.hpp
View File

@ -66,20 +66,15 @@ namespace Opm
/// manage control switching does not exist.
///
/// @param[in] W Existing wells object.
WellsManager(struct Wells* W, bool checkCellExistence=true);
explicit WellsManager(struct Wells* W, bool checkCellExistence=true);
/// Construct from input deck and grid.
/// The permeability argument may be zero if the input contain
/// well productivity indices, otherwise it must be given in
/// order to approximate these by the Peaceman formula.
WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence=true);
template<class CC, class F2C, class FC>
WellsManager(const Opm::EclipseGridParser& deck,
WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
int num_cells,
const int* global_cell,
const int* cart_dims,
@ -87,15 +82,12 @@ namespace Opm
CC begin_cell_centroids,
const F2C& f2c,
FC begin_face_centroids,
const double* permeability,
bool checkCellExistence=true);
const double* permeability);
WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence=true);
const double* permeability);
/// Destructor.
~WellsManager();
@ -148,25 +140,26 @@ namespace Opm
private:
template<class CC, class F2C, class FC>
void init(const Opm::EclipseGridParser& deck,
template<class CC, class C2F, class FC>
void init(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
int num_cells,
const int* global_cell,
const int* cart_dims,
int dimensions,
CC begin_cell_centroids,
const F2C& f2c,
const C2F& cell_to_faces,
FC begin_face_centroids,
const double* permeability);
// Disable copying and assignment.
WellsManager(const WellsManager& other, bool checkCellExistence=true);
WellsManager(const WellsManager& other);
WellsManager& operator=(const WellsManager& other);
static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
template<class C2F, class CC, class FC>
void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
const C2F& c2f,
const C2F& cell_to_faces,
const int* cart_dims,
FC begin_face_centroids,
CC begin_cell_centroids,
@ -185,7 +178,6 @@ namespace Opm
// Data
Wells* w_;
WellCollection well_collection_;
bool checkCellExistence_;
};
} // namespace Opm

583
opm/core/wells/WellsManager_impl.hpp
View File

@ -127,10 +127,8 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
if (checkCellExistence_)
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << well->name() << ')');
continue;
}
int cell = cgit->second;
PerfData pd;
@ -206,522 +204,90 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
}
template<class CC, class C2F, class FC>
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
int dimensions,
CC begin_cell_centroids,
const C2F& c2f,
const C2F& cell_to_faces,
FC begin_face_centroids,
const double* permeability,
bool checkCellExistence)
: w_(0), checkCellExistence_(checkCellExistence)
const double* permeability)
: w_(0)
{
init(deck, number_of_cells, global_cell, cart_dims, dimensions,
begin_cell_centroids, c2f, begin_face_centroids, permeability);
init(eclipseState, timeStep, number_of_cells, global_cell, cart_dims, dimensions,
begin_cell_centroids, cell_to_faces, begin_face_centroids, permeability);
}
/// Construct wells from deck.
template<class CC, class C2F, class FC>
void WellsManager::init(const Opm::EclipseGridParser& deck,
void WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
int dimensions,
CC begin_cell_centroids,
const C2F& c2f,
const C2F& cell_to_faces,
FC begin_face_centroids,
const double* permeability)
{
if (dimensions != 3) {
OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
}
// NOTE: Implementation copied and modified from dune-porsol's class BlackoilWells.
std::vector<std::string> keywords;
keywords.push_back("WELSPECS");
keywords.push_back("COMPDAT");
// keywords.push_back("WELTARG");
if (!deck.hasFields(keywords)) {
OPM_MESSAGE("Missing well keywords in deck, initializing no wells.");
if (eclipseState->getSchedule()->numWells() == 0) {
OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
return;
}
if (!(deck.hasField("WCONINJE") || deck.hasField("WCONPROD")) ) {
OPM_THROW(std::runtime_error, "Needed field is missing in file");
}
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(global_cell,
number_of_cells, cartesian_to_compressed);
// Obtain phase usage data.
PhaseUsage pu = phaseUsageFromDeck(deck);
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
// These data structures will be filled in this constructor,
// then used to initialize the Wells struct.
std::vector<std::string> well_names;
std::vector<WellData> well_data;
std::vector<std::vector<PerfData> > wellperf_data;
// For easy lookup:
std::map<std::string, int> well_names_to_index;
typedef std::map<std::string, int>::const_iterator WNameIt;
// Get WELSPECS data.
// It is allowed to have multiple lines corresponding to
// the same well, in which case the last one encountered
// is the one used.
const WELSPECS& welspecs = deck.getWELSPECS();
const int num_welspecs = welspecs.welspecs.size();
well_names.reserve(num_welspecs);
well_data.reserve(num_welspecs);
for (int w = 0; w < num_welspecs; ++w) {
// First check if this well has already been encountered.
// If so, we modify it's data instead of appending a new well
// to the well_data and well_names vectors.
const std::string& name = welspecs.welspecs[w].name_;
const double refdepth = welspecs.welspecs[w].datum_depth_BHP_;
WNameIt wit = well_names_to_index.find(name);
if (wit == well_names_to_index.end()) {
// New well, append data.
well_names_to_index[welspecs.welspecs[w].name_] = well_data.size();
well_names.push_back(name);
WellData wd;
// If negative (defaulted), set refdepth to a marker
// value, will be changed after getting perforation
// data to the centroid of the cell of the top well
// perforation.
wd.reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth;
wd.welspecsline = w;
well_data.push_back(wd);
} else {
// Existing well, change data.
const int wix = wit->second;
well_data[wix].reference_bhp_depth = (refdepth < 0.0) ? -1e100 : refdepth;
well_data[wix].welspecsline = w;
}
}
const int num_wells = well_data.size();
wellperf_data.resize(num_wells);
ScheduleConstPtr schedule = eclipseState->getSchedule();
std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
well_names.reserve(wells.size());
well_data.reserve(wells.size());
// global_cell is a map from compressed cells to Cartesian grid cells.
// We must make the inverse lookup.
const int* cpgdim = cart_dims;
std::map<int,int> cartesian_to_compressed;
if (global_cell) {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
}
else {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(i, i));
}
}
// Get COMPDAT data
// It is *not* allowed to have multiple lines corresponding to
// the same perforation!
const COMPDAT& compdat = deck.getCOMPDAT();
const int num_compdat = compdat.compdat.size();
for (int kw = 0; kw < num_compdat; ++kw) {
// Extract well name, or the part of the well name that
// comes before the '*'.
std::string name = compdat.compdat[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
// Look for well with matching name.
bool found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { // equal
// Extract corresponding WELSPECS defintion for
// purpose of default location specification.
const WelspecsLine& wspec = welspecs.welspecs[well_data[wix].welspecsline];
// We have a matching name.
int ix = compdat.compdat[kw].grid_ind_[0] - 1;
int jy = compdat.compdat[kw].grid_ind_[1] - 1;
int kz1 = compdat.compdat[kw].grid_ind_[2] - 1;
int kz2 = compdat.compdat[kw].grid_ind_[3] - 1;
if (ix < 0) {
// Defaulted I location. Extract from WELSPECS.
ix = wspec.I_ - 1;
}
if (jy < 0) {
// Defaulted J location. Extract from WELSPECS.
jy = wspec.J_ - 1;
}
if (kz1 < 0) {
// Defaulted KZ1. Use top layer.
kz1 = 0;
}
if (kz2 < 0) {
// Defaulted KZ2. Use bottom layer.
kz2 = cpgdim[2] - 1;
}
for (int kz = kz1; kz <= kz2; ++kz) {
int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz);
std::map<int, int>::const_iterator cgit =
cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
if(checkCellExistence_)
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' '
<< kz << " not found in grid (well = " << name << ')');
continue;
}
int cell = cgit->second;
PerfData pd;
pd.cell = cell;
if (compdat.compdat[kw].connect_trans_fac_ > 0.0) {
pd.well_index = compdat.compdat[kw].connect_trans_fac_;
} else {
double radius = 0.5*compdat.compdat[kw].diameter_;
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
std::array<double, 3> cubical = WellsManagerDetail::getCubeDim(c2f,
createWellsFromSpecs(wells, timeStep, cell_to_faces,
cart_dims,
begin_face_centroids,
begin_cell_centroids,
dimensions,
cell);
const double* cell_perm = &permeability[dimensions*dimensions*cell];
pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
compdat.compdat[kw].skin_factor_);
}
wellperf_data[wix].push_back(pd);
}
found = true;
break;
}
}
if (!found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << compdat.compdat[kw].well_
<< " in COMPDAT");
}
well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
setupWellControls(wells, timeStep, well_names, pu);
{
GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
well_collection_.addField(fieldGroup, timeStep, pu);
addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
}
// Set up reference depths that were defaulted. Count perfs.
int num_perfs = 0;
assert(dimensions == 3);
for (int w = 0; w < num_wells; ++w) {
num_perfs += wellperf_data[w].size();
if (well_data[w].reference_bhp_depth < 0.0) {
// It was defaulted. Set reference depth to minimum perforation depth.
double min_depth = 1e100;
int num_wperfs = wellperf_data[w].size();
for (int perf = 0; perf < num_wperfs; ++perf) {
double depth = UgGridHelpers::
getCoordinate(UgGridHelpers::increment(begin_cell_centroids,
wellperf_data[w][perf].cell,
dimensions), 2);
min_depth = std::min(min_depth, depth);
}
well_data[w].reference_bhp_depth = min_depth;
}
for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
well_collection_.addWell((*wellIter), timeStep, pu);
}
// Create the well data structures.
w_ = create_wells(pu.num_phases, num_wells, num_perfs);
if (!w_) {
OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
}
well_collection_.setWellsPointer(w_);
well_collection_.applyGroupControls();
// Classify wells
if (deck.hasField("WCONINJE")) {
const std::vector<WconinjeLine>& lines = deck.getWCONINJE().wconinje;
for (size_t i = 0 ; i < lines.size(); ++i) {
const std::map<std::string, int>::const_iterator it = well_names_to_index.find(lines[i].well_);
if (it != well_names_to_index.end()) {
const int well_index = it->second;
well_data[well_index].type = INJECTOR;
} else {
OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONINJE");
}
}
}
if (deck.hasField("WCONPROD")) {
const std::vector<WconprodLine>& lines = deck.getWCONPROD().wconprod;
for (size_t i = 0; i < lines.size(); ++i) {
const std::map<std::string, int>::const_iterator it = well_names_to_index.find(lines[i].well_);
if (it != well_names_to_index.end()) {
const int well_index = it->second;
well_data[well_index].type = PRODUCER;
} else {
OPM_THROW(std::runtime_error, "Unseen well name: " << lines[i].well_ << " first seen in WCONPROD");
}
}
}
// Add wells.
for (int w = 0; w < num_wells; ++w) {
const int w_num_perf = wellperf_data[w].size();
std::vector<int> perf_cells(w_num_perf);
std::vector<double> perf_prodind(w_num_perf);
for (int perf = 0; perf < w_num_perf; ++perf) {
perf_cells[perf] = wellperf_data[w][perf].cell;
perf_prodind[perf] = wellperf_data[w][perf].well_index;
}
const double* comp_frac = NULL;
// We initialize all wells with a null component fraction,
// and must (for injection wells) overwrite it later.
int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf,
comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_);
if (!ok) {
OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure.");
}
}
// Get WCONINJE data, add injection controls to wells.
// It is allowed to have multiple lines corresponding to
// the same well, in which case the last one encountered
// is the one used.
if (deck.hasField("WCONINJE")) {
const WCONINJE& wconinjes = deck.getWCONINJE();
const int num_wconinjes = wconinjes.wconinje.size();
for (int kw = 0; kw < num_wconinjes; ++kw) {
const WconinjeLine& wci_line = wconinjes.wconinje[kw];
// Extract well name, or the part of the well name that
// comes before the '*'.
std::string name = wci_line.well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
assert(well_data[wix].type == w_->type[wix]);
if (well_data[wix].type != INJECTOR) {
OPM_THROW(std::runtime_error, "Found WCONINJE entry for a non-injector well: " << well_names[wix]);
}
// Add all controls that are present in well.
// First we must clear existing controls, in case the
// current WCONINJE line is modifying earlier controls.
clear_well_controls(wix, w_);
int ok = 1;
int control_pos[5] = { -1, -1, -1, -1, -1 };
if (ok && wci_line.surface_flow_max_rate_ >= 0.0) {
control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
if (wci_line.injector_type_ == "WATER") {
distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
} else if (wci_line.injector_type_ == "OIL") {
distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (wci_line.injector_type_ == "GAS") {
distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
} else {
OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported."
"WellsManager only supports WATER, OIL and GAS injector types.");
}
ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_,
distr, wix, w_);
}
if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) {
control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
if (wci_line.injector_type_ == "WATER") {
distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
} else if (wci_line.injector_type_ == "OIL") {
distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (wci_line.injector_type_ == "GAS") {
distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
} else {
OPM_THROW(std::runtime_error, "Injector type " << wci_line.injector_type_ << " not supported."
"WellsManager only supports WATER, OIL and GAS injector types.");
}
ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_,
distr, wix, w_);
}
if (ok && wci_line.BHP_limit_ > 0.0) {
control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[wix]);
ok = append_well_controls(BHP, wci_line.BHP_limit_,
NULL, wix, w_);
}
if (ok && wci_line.THP_limit_ > 0.0) {
OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]);
}
if (!ok) {
OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]);
}
WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode(wci_line.control_mode_);
int cpos = control_pos[mode];
if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) {
OPM_THROW(std::runtime_error, "Control for " << wci_line.control_mode_ << " not specified in well " << well_names[wix]);
}
// We need to check if the well is shut or not
if (wci_line.open_shut_flag_ == "SHUT") {
cpos = ~cpos;
}
set_current_control(wix, cpos, w_);
// Set well component fraction.
double cf[3] = { 0.0, 0.0, 0.0 };
if (wci_line.injector_type_[0] == 'W') {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well.");
}
cf[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
} else if (wci_line.injector_type_[0] == 'O') {
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well.");
}
cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
} else if (wci_line.injector_type_[0] == 'G') {
if (!pu.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well.");
}
cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
}
std::copy(cf, cf + pu.num_phases, w_->comp_frac + wix*pu.num_phases);
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << wci_line.well_
<< " in WCONINJE");
}
}
}
// Get WCONPROD data
// It is allowed to have multiple lines corresponding to
// the same well, in which case the last one encountered
// is the one used.
if (deck.hasField("WCONPROD")) {
const WCONPROD& wconprods = deck.getWCONPROD();
const int num_wconprods = wconprods.wconprod.size();
for (int kw = 0; kw < num_wconprods; ++kw) {
const WconprodLine& wcp_line = wconprods.wconprod[kw];
std::string name = wcp_line.well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
assert(well_data[wix].type == w_->type[wix]);
if (well_data[wix].type != PRODUCER) {
OPM_THROW(std::runtime_error, "Found WCONPROD entry for a non-producer well: " << well_names[wix]);
}
// Add all controls that are present in well.
// First we must clear existing controls, in case the
// current WCONPROD line is modifying earlier controls.
clear_well_controls(wix, w_);
int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
int ok = 1;
if (ok && wcp_line.oil_max_rate_ >= 0.0) {
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
}
control_pos[WellsManagerDetail::ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_,
distr, wix, w_);
}
if (ok && wcp_line.water_max_rate_ >= 0.0) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
}
control_pos[WellsManagerDetail::ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_,
distr, wix, w_);
}
if (ok && wcp_line.gas_max_rate_ >= 0.0) {
if (!pu.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
}
control_pos[WellsManagerDetail::ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_,
distr, wix, w_);
}
if (ok && wcp_line.liquid_max_rate_ >= 0.0) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified.");
}
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified.");
}
control_pos[WellsManagerDetail::ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 0.0, 0.0, 0.0 };
distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_,
distr, wix, w_);
}
if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) {
control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[wix]);
double distr[3] = { 1.0, 1.0, 1.0 };
ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_,
distr, wix, w_);
}
if (ok && wcp_line.BHP_limit_ > 0.0) {
control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[wix]);
ok = append_well_controls(BHP, wcp_line.BHP_limit_,
NULL, wix, w_);
}
if (ok && wcp_line.THP_limit_ > 0.0) {
OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[wix]);
}
if (!ok) {
OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[wix]);
}
WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(wcp_line.control_mode_);
int cpos = control_pos[mode];
if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[wix]);
}
// If it's shut, we complement the cpos
if (wcp_line.open_shut_flag_ == "SHUT") {
cpos = ~cpos; // So we can easily retrieve the cpos later
}
set_current_control(wix, cpos, w_);
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << wcp_line.well_
<< " in WCONPROD");
}
}
}
// Get WELTARG data
if (deck.hasField("WELTARG")) {
OPM_THROW(std::runtime_error, "We currently do not handle WELTARG.");
/*
const WELTARG& weltargs = deck.getWELTARG();
const int num_weltargs = weltargs.weltarg.size();
for (int kw = 0; kw < num_weltargs; ++kw) {
std::string name = weltargs.weltarg[kw].well_;
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix = 0; wix < num_wells; ++wix) {
if (well_names[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data[wix].target = weltargs.weltarg[kw].new_value_;
break;
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << weltargs.weltarg[kw].well_
<< " in WELTARG");
}
}
*/
}
setupGuideRates(wells, timeStep, well_data, well_names_to_index);
// Debug output.
#define EXTRA_OUTPUT
@ -743,77 +309,6 @@ void WellsManager::init(const Opm::EclipseGridParser& deck,
}
*/
#endif
if (deck.hasField("WELOPEN")) {
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
OPM_THROW(std::runtime_error, "Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
const int index = it->second;
if (line.openshutflag_ == "SHUT") {
well_controls_shut_well( w_->ctrls[index] );
} else if (line.openshutflag_ == "OPEN") {
well_controls_open_well( w_->ctrls[index] );
} else {
OPM_THROW(std::runtime_error, "Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
}
}
// Build the well_collection_ well group hierarchy.
if (deck.hasField("GRUPTREE")) {
std::cout << "Found gruptree" << std::endl;
const GRUPTREE& gruptree = deck.getGRUPTREE();
std::map<std::string, std::string>::const_iterator it = gruptree.tree.begin();
for( ; it != gruptree.tree.end(); ++it) {
well_collection_.addChild(it->first, it->second, deck);
}
}
for (size_t i = 0; i < welspecs.welspecs.size(); ++i) {
WelspecsLine line = welspecs.welspecs[i];
well_collection_.addChild(line.name_, line.group_, deck);
}
// Set the guide rates:
if (deck.hasField("WGRUPCON")) {
std::cout << "Found Wgrupcon" << std::endl;
WGRUPCON wgrupcon = deck.getWGRUPCON();
const std::vector<WgrupconLine>& lines = wgrupcon.wgrupcon;
std::cout << well_collection_.getLeafNodes().size() << std::endl;
for (size_t i = 0; i < lines.size(); i++) {
std::string name = lines[i].well_;
const int wix = well_names_to_index[name];
WellNode& wellnode = *well_collection_.getLeafNodes()[wix];
assert(wellnode.name() == name);
if (well_data[wix].type == PRODUCER) {
wellnode.prodSpec().guide_rate_ = lines[i].guide_rate_;
if (lines[i].phase_ == "OIL") {
wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL;
} else {
OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for producer "
<< name << " in WGRUPCON, cannot handle.");
}
} else if (well_data[wix].type == INJECTOR) {
wellnode.injSpec().guide_rate_ = lines[i].guide_rate_;
if (lines[i].phase_ == "RAT") {
wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT;
} else {
OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for injector "
<< name << " in WGRUPCON, cannot handle.");
}
} else {
OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << name);
}
}
}
well_collection_.setWellsPointer(w_);
well_collection_.applyGroupControls();
}
} // end namespace Opm

186
tests/test_linearsolver.cpp Normal file
View File

@ -0,0 +1,186 @@
/*
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
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/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE OPM-IterativeSolverTest
#include <boost/test/unit_test.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <dune/common/version.hh>
#include <memory>
#include <cstdlib>
#include <string>
struct MyMatrix
{
MyMatrix(int rows, int nnz)
: data(nnz, 0.0), rowStart(rows+1, -1),
colIndex(nnz, -1)
{}
MyMatrix()
: data(), rowStart(), colIndex()
{}
std::vector<double> data;
std::vector<int> rowStart;
std::vector<int> colIndex;
};
std::shared_ptr<MyMatrix> createLaplacian(int N)
{
MyMatrix* mm=new MyMatrix(N*N, N*N*5);
int nnz=0;
mm->rowStart[0]=0;
for(int row=0; row<N*N; row++)
{
int x=row%N;
int y=row/N;
if(y>0)
{
mm->colIndex[nnz]=row-N;
mm->data[nnz++]=-1;
}
if(x>0)
{
mm->colIndex[nnz]=row-1;
mm->data[nnz++]=-1;
}
mm->colIndex[nnz]=row;
mm->data[nnz++]=4;
if(x<N-1)
{
mm->colIndex[nnz]=row+1;
mm->data[nnz++]=-1;
}
if(y<N-1)
{
mm->colIndex[nnz]=row+N;
mm->data[nnz++]=-1;
}
mm->rowStart[row+1]=nnz;
}
mm->data.resize(nnz);
mm->colIndex.resize(nnz);
return std::shared_ptr<MyMatrix>(mm);
}
void createRandomVectors(int NN, std::vector<double>& x, std::vector<double>& b,
const MyMatrix& mat)
{
x.resize(NN);
for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry)
*entry=((double) (rand()%100))/10.0;
b.resize(NN);
std::fill(b.begin(), b.end(), 0.0);
// Construct the right hand side as b=A*x
for(std::size_t row=0; row<mat.rowStart.size()-1; ++row)
{
for(int i=mat.rowStart[row], end=mat.rowStart[row+1]; i!=end; ++i)
{
b[row]+= mat.data[i]*x[mat.colIndex[i]];
}
}
}
void run_test(const Opm::parameter::ParameterGroup& param)
{
int N=4;
auto mat = createLaplacian(N);
std::vector<double> x(N*N), b(N*N);
createRandomVectors(100*100, x, b, *mat);
std::vector<double> exact(x);
std::fill(x.begin(), x.end(), 0.0);
Opm::LinearSolverFactory ls(param);
ls.solve(N*N, mat->data.size(), &(mat->rowStart[0]),
&(mat->colIndex[0]), &(mat->data[0]), &(b[0]),
&(x[0]));
}
BOOST_AUTO_TEST_CASE(DefaultTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
#ifdef HAVE_DUNE_ISTL
BOOST_AUTO_TEST_CASE(CGAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("1"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(CGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("0"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(BiCGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("2"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
BOOST_AUTO_TEST_CASE(FastAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("3"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(KAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("4"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
run_test(param);
}
#endif
#endif

16
tests/test_wellsgroup.cpp
View File

@ -59,16 +59,18 @@ BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) {
std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
if (well->isInjector(2)) {
BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
const WellInjectionProperties& properties = well->getInjectionProperties(2);
BOOST_CHECK_EQUAL(properties.surfaceInjectionRate, wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->injSpec().BHP_limit_);
BOOST_CHECK_EQUAL(properties.reservoirInjectionRate, wellsGroup->injSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
}
if (well->isProducer(2)) {
BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_);
const WellProductionProperties& properties = well->getProductionProperties(2);
BOOST_CHECK_EQUAL(properties.ResVRate, wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->prodSpec().BHP_limit_);
BOOST_CHECK_EQUAL(properties.OilRate, wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(properties.WaterRate, wellsGroup->prodSpec().water_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_);
}
}

95
tests/test_wellsmanager.cpp
View File

@ -173,30 +173,6 @@ void check_controls_epoch1( struct WellControls ** ctrls) {
}
}
BOOST_AUTO_TEST_CASE(Constructor_Works) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
const Wells* wells = wellsManager.c_wells();
wells_static_check( wells );
check_controls_epoch0( wells->ctrls );
}
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
const Wells* wells = wellsManager.c_wells();
wells_static_check( wells );
check_controls_epoch1( wells->ctrls );
}
}
BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Opm::ParserPtr parser(new Opm::Parser());
@ -205,70 +181,16 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 0" << std::endl;
wells_static_check( wellsManager.c_wells() );
std::cout << "Checking old well structure, epoch 0" << std::endl;
wells_static_check( oldWellsManager.c_wells() );
check_controls_epoch0( wellsManager.c_wells()->ctrls );
BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells() , false));
}
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 1" << std::endl;
wells_static_check( wellsManager.c_wells() );
std::cout << "Checking old well structure, epoch 1" << std::endl;
wells_static_check( oldWellsManager.c_wells() );
check_controls_epoch1( wellsManager.c_wells()->ctrls );
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(),false));
}
}
BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_expanded.data")));
Opm::EclipseGridParser Deck("wells_manager_data_expanded.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false));
}
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
}
Deck.setCurrentEpoch(2);
{
Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
}
}
@ -277,12 +199,11 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
BOOST_AUTO_TEST_CASE(WellsEqual) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
Deck.setCurrentEpoch(0);
Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
Deck.setCurrentEpoch(1);
Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
BOOST_CHECK( wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) );
BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) );
@ -293,11 +214,11 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
Deck.setCurrentEpoch(1);
Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false));
BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false));