diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 77aaf4ca7..c32ad02ee 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -30,6 +30,7 @@ // TODO: clean up includes. #include +#include #include #include #include @@ -39,6 +40,10 @@ #include #include +#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3) +#include +#endif + #include #include @@ -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) diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index 47108ef23..fb6f8f94f 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -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."); diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 39fbdbd3e..090676ee5 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -23,6 +23,7 @@ #include #include #include +#include namespace Opm { @@ -45,26 +46,55 @@ 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 ); - } - - // 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)) { - const double rate_target = well_controls_get_current_target(ctrl); - const double * distr = well_controls_get_current_distr( 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] = rate_target * distr[p]; + wellrates_[np*w + p] = 0.0; + } + // 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]; } } } diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index c512efaee..886532cea 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -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) { diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index ee45b6ed5..050e18712 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -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 wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage)); diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 24ad6c288..3beab7567 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -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 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 well_names; - std::vector well_data; - - - // For easy lookup: - std::map well_names_to_index; - - ScheduleConstPtr schedule = eclipseState->getSchedule(); - std::vector 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) { - 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) { - 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) { - 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); + { + 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 (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 (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_); } diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index c92ba3a2c..41f567573 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -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 - 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 - void init(const Opm::EclipseGridParser& deck, + template + 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& cartesian_to_compressed ); void setupWellControls(std::vector& wells, size_t timeStep, std::vector& well_names, const PhaseUsage& phaseUsage); template void createWellsFromSpecs( std::vector& 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 diff --git a/opm/core/wells/WellsManager_impl.hpp b/opm/core/wells/WellsManager_impl.hpp index 203b10c25..058f32ea1 100644 --- a/opm/core/wells/WellsManager_impl.hpp +++ b/opm/core/wells/WellsManager_impl.hpp @@ -127,10 +127,8 @@ void WellsManager::createWellsFromSpecs(std::vector& wells, size_t int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); std::map::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& wells, size_t } template -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. +/// Construct wells from deck. template -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 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 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 well_names; std::vector well_data; - std::vector > wellperf_data; + // For easy lookup: std::map well_names_to_index; - typedef std::map::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 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 cartesian_to_compressed; + createWellsFromSpecs(wells, timeStep, cell_to_faces, + cart_dims, + begin_face_centroids, + begin_cell_centroids, + dimensions, + well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability); + setupWellControls(wells, timeStep, well_names, pu); - 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)); - } + { + 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); } - // 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::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 cubical = WellsManagerDetail::getCubeDim(c2f, - begin_face_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"); - } + for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) { + well_collection_.addWell((*wellIter), 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; - } - } + well_collection_.setWellsPointer(w_); + well_collection_.applyGroupControls(); - // 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."); - } - // Classify wells - if (deck.hasField("WCONINJE")) { - const std::vector& lines = deck.getWCONINJE().wconinje; - for (size_t i = 0 ; i < lines.size(); ++i) { - const std::map::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& lines = deck.getWCONPROD().wconprod; - for (size_t i = 0; i < lines.size(); ++i) { - const std::map::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 perf_cells(w_num_perf); - std::vector 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::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::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& 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 diff --git a/tests/test_linearsolver.cpp b/tests/test_linearsolver.cpp new file mode 100644 index 000000000..184b7765d --- /dev/null +++ b/tests/test_linearsolver.cpp @@ -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 . +*/ + +#include + +#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 + +#include +#include + +#include +#include +#include +#include + +struct MyMatrix +{ + MyMatrix(int rows, int nnz) + : data(nnz, 0.0), rowStart(rows+1, -1), + colIndex(nnz, -1) + {} + MyMatrix() + : data(), rowStart(), colIndex() + {} + + std::vector data; + std::vector rowStart; + std::vector colIndex; +}; + +std::shared_ptr createLaplacian(int N) +{ + MyMatrix* mm=new MyMatrix(N*N, N*N*5); + int nnz=0; + mm->rowStart[0]=0; + for(int row=0; row0) + { + 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(xcolIndex[nnz]=row+1; + mm->data[nnz++]=-1; + } + if(ycolIndex[nnz]=row+N; + mm->data[nnz++]=-1; + } + mm->rowStart[row+1]=nnz; + } + mm->data.resize(nnz); + mm->colIndex.resize(nnz); + return std::shared_ptr(mm); +} + +void createRandomVectors(int NN, std::vector& x, std::vector& 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 x(N*N), b(N*N); + createRandomVectors(100*100, x, b, *mat); + std::vector 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 diff --git a/tests/test_wellsgroup.cpp b/tests/test_wellsgroup.cpp index 9ba9e0ef3..46fd8eb0c 100644 --- a/tests/test_wellsgroup.cpp +++ b/tests/test_wellsgroup.cpp @@ -59,16 +59,18 @@ BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) { std::shared_ptr 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_); } } diff --git a/tests/test_wellsmanager.cpp b/tests/test_wellsmanager.cpp index 8cb097e2e..2141e4bf8 100644 --- a/tests/test_wellsmanager.cpp +++ b/tests/test_wellsmanager.cpp @@ -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));