20#ifndef OPM_BLACKOILWELLS_HEADER_INCLUDED
21#define OPM_BLACKOILWELLS_HEADER_INCLUDED
24#include <dune/grid/CpGrid.hpp>
25#include <opm/common/ErrorMacros.hpp>
26#include <opm/core/utility/SparseTable.hpp>
28#include <opm/parser/eclipse/Deck/Deck.hpp>
29#include <dune/common/fvector.hh>
41 void init(Opm::DeckConstPtr deck,
42 const Dune::CpGrid& grid,
51 double target(
int wellnum)
const;
54 int wellCell(
int wellnum,
int perfnum)
const;
55 double wellIndex(
int wellnum,
int perfnum)
const;
59 const std::vector<double>& well_perf_pressures,
60 const std::vector<double>& well_perf_fluxes);
73 if (well_report_ == 0x0)
96 double computeWellIndex(
double radius,
const Dune::FieldVector<double, 3>& cubical,
100 std::vector<WellData> well_data_;
101 struct PerfData {
int cell;
double well_index; };
102 Opm::SparseTable<PerfData> perf_data_;
103 std::vector<double> well_cell_flux_;
104 std::vector<double> well_cell_pressure_;
105 Dune::FieldVector<double, 3> injection_mixture_;
106 std::vector<std::string> well_names_;
109 BlackoilWells::WellReport* BlackoilWells::WellReport::well_report_ = 0x0;
117 int prod_control_mode(
const std::string& control);
118 int inje_control_mode(
const std::string& control);
119 template<
class gr
id_t>
120 const Dune::FieldVector<double,3> getCubeDim(
const grid_t& grid,
int cell);
125 const Dune::CpGrid& grid,
128 if (!deck->hasKeyword(
"WELSPECS")) {
129 OPM_MESSAGE(
"Missing keyword \"WELSPECS\" in deck, initializing no wells.");
132 if (!deck->hasKeyword(
"COMPDAT")) {
133 OPM_MESSAGE(
"Missing keyword \"COMPDAT\" in deck, initializing no wells.");
136 if (!(deck->hasKeyword(
"WCONINJE") || deck->hasKeyword(
"WCONPROD")) ) {
137 OPM_THROW(std::runtime_error,
"Needed field is missing in file");
142 const auto& welspecsKeyword = deck->getKeyword(
"WELSPECS");
143 const int num_welspecs = welspecsKeyword.size();
144 well_names_.reserve(num_welspecs);
145 well_data_.reserve(num_welspecs);
146 for (
int i=0; i<num_welspecs; ++i) {
147 well_names_.push_back(welspecsKeyword.getRecord(i).getItem(
"WELL").get< std::string >(0));
149 well_data_.push_back(wd);
150 well_data_.back().reference_bhp_depth =
151 welspecsKeyword.getRecord(i).getItem(
"REF_DEPTH").getSIDouble(0);
155 const auto& compdatKeyword = deck->getKeyword(
"COMPDAT");
156 const int num_compdats = compdatKeyword.size();
157 std::vector<std::vector<PerfData> > wellperf_data(num_welspecs);
158 for (
int kw=0; kw<num_compdats; ++kw) {
159 const auto& compdatRecord = compdatKeyword.getRecord(kw);
160 std::string name = compdatRecord.getItem(
"WELL").get< std::string >(0);
161 std::string::size_type len = name.find(
'*');
162 if (len != std::string::npos) {
163 name = name.substr(0, len);
167 const std::vector<int>& global_cell = grid.globalCell();
168 const std::array<int, 3>& cpgdim = grid.logicalCartesianSize();
169 std::map<int,int> cartesian_to_compressed;
170 for (
int i=0; i<int(global_cell.size()); ++i) {
171 cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
174 for (
int wix=0; wix<num_welspecs; ++wix) {
175 if (well_names_[wix].compare(0,len, name) == 0) {
176 int ix = compdatRecord.getItem(
"I").get<
int >(0) - 1;
177 int jy = compdatRecord.getItem(
"J").get<
int >(0) - 1;
178 int kz1 = compdatRecord.getItem(
"K1").get<
int >(0) - 1;
179 int kz2 = compdatRecord.getItem(
"K2").get<
int >(0) - 1;
180 for (
int kz = kz1; kz <= kz2; ++kz) {
181 int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz);
182 std::map<int, int>::const_iterator cgit =
183 cartesian_to_compressed.find(cart_grid_indx);
184 if (cgit == cartesian_to_compressed.end()) {
185 OPM_THROW(std::runtime_error,
"Cell with i,j,k indices " << ix <<
' ' << jy <<
' '
186 << kz <<
" not found!");
188 int cell = cgit->second;
191 if (compdatRecord.getItem(
"CF").getSIDouble(0) > 0.0) {
192 pd.well_index = compdatRecord.getItem(
"CF").getSIDouble(0);
194 double radius = 0.5*compdatRecord.getItem(
"DIAMETER").getSIDouble(0);
196 radius = 0.5*unit::feet;
197 OPM_MESSAGE(
"Warning: Well bore internal radius set to " << radius);
199 Dune::FieldVector<double, 3> cubical = getCubeDim(grid, cell);
201 pd.well_index = computeWellIndex(radius, cubical, permeability,
202 compdatRecord.getItem(
"SKIN").getSIDouble(0));
204 wellperf_data[wix].push_back(pd);
211 OPM_THROW(std::runtime_error,
"Undefined well name: " << compdatRecord.getItem(
"WELL").get< std::string >(0)
215 for (
int w = 0; w < num_welspecs; ++w) {
216 perf_data_.appendRow(wellperf_data[w].begin(), wellperf_data[w].end());
217 if (well_data_[w].reference_bhp_depth == -1e100) {
219 double min_depth = 1e100;
220 int num_wperfs = wellperf_data[w].size();
221 for (
int perf = 0; perf < num_wperfs; ++perf) {
222 double depth = grid.cellCentroid(wellperf_data[w][perf].cell)[2];
223 min_depth = std::min(min_depth, depth);
225 well_data_[w].reference_bhp_depth = min_depth;
230 injection_mixture_ = 0.0;
231 if (deck->hasKeyword(
"WCONINJE")) {
232 const auto& wconinjeKeyword = deck->getKeyword(
"WCONINJE");
233 const int num_wconinjes = wconinjeKeyword.size();
234 int injector_component = -1;
235 for (
int kw=0; kw<num_wconinjes; ++kw) {
236 const auto& wconinjeRecord = wconinjeKeyword.getRecord(kw);
237 std::string name = wconinjeRecord.getItem(
"WELL").get< std::string >(0);
238 std::string::size_type len = name.find(
'*');
239 if (len != std::string::npos) {
240 name = name.substr(0, len);
243 bool well_found =
false;
244 for (
int wix=0; wix<num_welspecs; ++wix) {
245 if (well_names_[wix].compare(0,len, name) == 0) {
248 int m = inje_control_mode(wconinjeRecord.getItem(
"CMODE").get< std::string >(0));
251 well_data_[wix].control =
Rate;
253 well_data_[wix].target = wconinjeRecord.getItem(
"RATE").get<
double >(0);
256 well_data_[wix].control =
Rate;
258 well_data_[wix].target = wconinjeRecord.getItem(
"RESV").get<
double >(0);
262 well_data_[wix].target = wconinjeRecord.getItem(
"BHP").getSIDouble(0);
266 well_data_[wix].target = wconinjeRecord.getItem(
"THP").getSIDouble(0);
269 OPM_THROW(std::runtime_error,
"Unknown well control mode; WCONIJE = "
270 << wconinjeRecord.getItem(
"CMODE").get< std::string >(0)
271 <<
" in input file");
274 if (wconinjeRecord.getItem(
"TYPE").get< std::string >(0) ==
"WATER") {
276 }
else if (wconinjeRecord.getItem(
"TYPE").get< std::string >(0) ==
"OIL") {
278 }
else if (wconinjeRecord.getItem(
"TYPE").get< std::string >(0) ==
"GAS") {
281 if (itp == -1 || (injector_component != -1 && itp != injector_component)) {
283 OPM_THROW(std::runtime_error,
"Error in injector specification, found no known fluid type.");
285 OPM_THROW(std::runtime_error,
"Error in injector specification, we can only handle a single injection fluid.");
288 injector_component = itp;
293 OPM_THROW(std::runtime_error,
"Undefined well name: " << wconinjeRecord.getItem(
"WELL").get< std::string >(0)
297 if (injector_component != -1) {
298 injection_mixture_[injector_component] = 1.0;
307 injection_mixture_[
Oil] = 1.0;
311 if (deck->hasKeyword(
"WCONPROD")) {
312 const auto& wconprodKeyword = deck->getKeyword(
"WCONPROD");
313 const int num_wconprods = wconprodKeyword.size();
314 for (
int kw=0; kw<num_wconprods; ++kw) {
315 const auto& wconprodRecord = wconprodKeyword.getRecord(kw);
316 std::string name = wconprodRecord.getItem(
"WELL").get< std::string >(0);
317 std::string::size_type len = name.find(
'*');
318 if (len != std::string::npos) {
319 name = name.substr(0, len);
322 bool well_found =
false;
323 for (
int wix=0; wix<num_welspecs; ++wix) {
324 if (well_names_[wix].compare(0,len, name) == 0) {
327 int m = prod_control_mode(wconprodRecord.getItem(
"CMODE").get< std::string >(0));
330 well_data_[wix].control =
Rate;
331 well_data_[wix].target = wconprodRecord.getItem(
"ORAT").getSIDouble(0);
334 well_data_[wix].control =
Rate;
335 well_data_[wix].target = wconprodRecord.getItem(
"WRAT").getSIDouble(0);
338 well_data_[wix].control =
Rate;
339 well_data_[wix].target = wconprodRecord.getItem(
"GRAT").getSIDouble(0);
342 well_data_[wix].control =
Rate;
343 well_data_[wix].target = wconprodRecord.getItem(
"LRAT").getSIDouble(0);
346 well_data_[wix].control =
Rate;
347 well_data_[wix].target = wconprodRecord.getItem(
"RESV").getSIDouble(0);
351 well_data_[wix].target = wconprodRecord.getItem(
"BHP").getSIDouble(0);
355 well_data_[wix].target = wconprodRecord.getItem(
"THP").getSIDouble(0);
358 OPM_THROW(std::runtime_error,
"Unknown well control mode; WCONPROD = "
359 << wconprodRecord.getItem(
"CMODE").get< std::string >(0)
360 <<
" in input file");
365 OPM_THROW(std::runtime_error,
"Undefined well name: " << wconprodRecord.getItem(
"WELL").get< std::string >(0)
372 if (deck->hasKeyword(
"WELTARG")) {
373 const auto& weltargKeyword = deck->getKeyword(
"WELTARG");
374 const int num_weltargs = weltargKeyword.size();
375 for (
int kw=0; kw<num_weltargs; ++kw) {
376 const auto& weltargRecord = weltargKeyword.getRecord(kw);
377 std::string name = weltargRecord.getItem(
"WELL").get< std::string >(0);
378 std::string::size_type len = name.find(
'*');
379 if (len != std::string::npos) {
380 name = name.substr(0, len);
382 bool well_found =
false;
383 for (
int wix=0; wix<num_welspecs; ++wix) {
384 if (well_names_[wix].compare(0,len, name) == 0) {
387 well_data_[wix].target = weltargRecord.getItem(
"NEW_VALUE").get<
double >(0);
392 OPM_THROW(std::runtime_error,
"Undefined well name: " << weltargRecord.getItem(
"WELL").get< std::string >(0)
399 std::cout <<
"\t WELL DATA" << std::endl;
400 for(
int i=0; i< int(well_data_.size()); ++i) {
401 std::cout << i <<
": " << well_data_[i].type <<
" "
402 << well_data_[i].control <<
" " << well_data_[i].target
406 std::cout <<
"\n\t PERF DATA" << std::endl;
407 for(
int i=0; i< int(perf_data_.size()); ++i) {
408 for(
int j=0; j< int(perf_data_[i].size()); ++j) {
409 std::cout << i <<
": " << perf_data_[i][j].cell <<
" "
410 << perf_data_[i][j].well_index << std::endl;
415 well_cell_pressure_.resize(grid.numCells(), -1e100);
416 well_cell_flux_.resize(grid.numCells(), 0.0);
421 return well_data_.size();
426 return well_data_[wellnum].type;
431 return well_data_[wellnum].control;
436 return well_data_[wellnum].target;
441 return well_data_[wellnum].reference_bhp_depth;
446 return perf_data_[wellnum].size();
451 return perf_data_[wellnum][perfnum].cell;
456 return perf_data_[wellnum][perfnum].well_index;
460 const std::vector<double>& well_perf_pressures,
461 const std::vector<double>& well_perf_fluxes)
464 assert(perf_data_.dataSize() ==
int(well_perf_pressures.size()));
465 well_cell_pressure_.resize(num_cells, -1e100);
466 well_cell_flux_.resize(num_cells, 0.0);
468 for (
int w = 0; w <
numWells(); ++w) {
471 well_cell_pressure_[cell] = well_perf_pressures[pcount];
472 well_cell_flux_[cell] = well_perf_fluxes[pcount];
476 assert(pcount == perf_data_.dataSize());
481 return well_cell_flux_[cell];
486 return well_cell_pressure_[cell];
491 return injection_mixture_;
494 inline double BlackoilWells::computeWellIndex(
double radius,
495 const Dune::FieldVector<double, 3>& cubical,
497 double skin_factor)
const
509 double effective_perm = sqrt(permeability(0,0) * permeability(1,1));
513 assert(permeability(0,0) > 0.0);
514 assert(permeability(1,1) > 0.0);
515 double kxoy = permeability(0,0) / permeability(1,1);
516 double kyox = permeability(1,1) / permeability(0,0);
517 double r0_denominator = pow(kyox, 0.25) + pow(kxoy, 0.25);
518 double r0_numerator = sqrt((sqrt(kyox)*cubical[0]*cubical[0]) +
519 (sqrt(kxoy)*cubical[1]*cubical[1]));
520 assert(r0_denominator > 0.0);
521 double r0 = 0.28 * r0_numerator / r0_denominator;
522 assert(radius > 0.0);
525 std::cout <<
"ERROR: Too big well radius detected.";
526 std::cout <<
"Specified well radius is " << radius
527 <<
" while r0 is " << r0 <<
".\n";
530 const double two_pi = 6.2831853071795864769252867665590057683943387987502116419498;
532 double wi_denominator = log(r0 / radius) + skin_factor;
533 double wi_numerator = two_pi * cubical[2];
534 assert(wi_denominator > 0.0);
535 double wi = effective_perm * wi_numerator / wi_denominator;
547 int prod_control_mode(
const std::string& control){
548 const int num_prod_control_modes = 8;
549 static std::string prod_control_modes[num_prod_control_modes] =
550 {std::string(
"ORAT"), std::string(
"WRAT"), std::string(
"GRAT"),
551 std::string(
"LRAT"), std::string(
"RESV"), std::string(
"BHP"),
552 std::string(
"THP"), std::string(
"GRUP") };
554 for (
int i=0; i<num_prod_control_modes; ++i) {
555 if (control == prod_control_modes[i]) {
563 OPM_THROW(std::runtime_error,
"Unknown well control mode = " << control <<
" in input file");
567 int inje_control_mode(
const std::string& control)
569 const int num_inje_control_modes = 5;
570 static std::string inje_control_modes[num_inje_control_modes] =
571 {std::string(
"RATE"), std::string(
"RESV"), std::string(
"BHP"),
572 std::string(
"THP"), std::string(
"GRUP") };
574 for (
int i=0; i<num_inje_control_modes; ++i) {
575 if (control == inje_control_modes[i]) {
584 OPM_THROW(std::runtime_error,
"Unknown well control mode = " << control <<
" in input file");
589 template<
class gr
id_t>
590 const Dune::FieldVector<double,3> getCubeDim(
const grid_t& grid,
int cell)
592 Dune::FieldVector<double, 3> cube;
593 int num_local_faces = grid.numCellFaces(cell);
594 std::vector<double> x(num_local_faces);
595 std::vector<double> y(num_local_faces);
596 std::vector<double> z(num_local_faces);
597 for (
int lf=0; lf<num_local_faces; ++ lf) {
598 int face = grid.cellFace(cell,lf);
599 const Dune::FieldVector<double,3>& centroid =
600 grid.faceCentroid(face);
605 cube[0] = *max_element(x.begin(), x.end()) -
606 *min_element(x.begin(), x.end());
607 cube[1] = *max_element(y.begin(), y.end()) -
608 *min_element(y.begin(), y.end());
609 cube[2] = *max_element(z.begin(), z.end()) -
610 *min_element(z.begin(), z.end());
Definition: BlackoilDefs.hpp:31
Dune::FieldVector< Scalar, numComponents > CompVec
Definition: BlackoilDefs.hpp:40
@ Gas
Definition: BlackoilDefs.hpp:36
@ Oil
Definition: BlackoilDefs.hpp:36
@ Water
Definition: BlackoilDefs.hpp:36
Definition: BlackoilWells.hpp:69
static WellReport * report()
Definition: BlackoilWells.hpp:71
std::vector< double > cellPressure
Definition: BlackoilWells.hpp:85
std::vector< BlackoilDefs::CompVec > massRate
Definition: BlackoilWells.hpp:86
void clearAll()
Definition: BlackoilWells.hpp:77
WellReport()
Definition: BlackoilWells.hpp:89
std::vector< double > perfPressure
Definition: BlackoilWells.hpp:84
std::vector< int > cellId
Definition: BlackoilWells.hpp:87
Definition: BlackoilWells.hpp:39
double target(int wellnum) const
Definition: BlackoilWells.hpp:434
double wellToReservoirFlux(int cell) const
Definition: BlackoilWells.hpp:479
WellType
Definition: BlackoilWells.hpp:47
@ Injector
Definition: BlackoilWells.hpp:47
@ Producer
Definition: BlackoilWells.hpp:47
int numWells() const
Definition: BlackoilWells.hpp:419
double referenceDepth(int wellnum) const
Definition: BlackoilWells.hpp:439
double perforationPressure(int cell) const
Definition: BlackoilWells.hpp:484
double wellIndex(int wellnum, int perfnum) const
Definition: BlackoilWells.hpp:454
void init(Opm::DeckConstPtr deck, const Dune::CpGrid &grid, const Opm::Rock< 3 > &rock)
Definition: BlackoilWells.hpp:124
WellType type(int wellnum) const
Definition: BlackoilWells.hpp:424
WellControl control(int wellnum) const
Definition: BlackoilWells.hpp:429
void update(int num_cells, const std::vector< double > &well_perf_pressures, const std::vector< double > &well_perf_fluxes)
Definition: BlackoilWells.hpp:459
int numPerforations(int wellnum) const
Definition: BlackoilWells.hpp:444
WellControl
Definition: BlackoilWells.hpp:49
@ Rate
Definition: BlackoilWells.hpp:49
@ Pressure
Definition: BlackoilWells.hpp:49
CompVec injectionMixture(int cell) const
Definition: BlackoilWells.hpp:489
int wellCell(int wellnum, int perfnum) const
Definition: BlackoilWells.hpp:449
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:80
PermTensor permeability(int cell_index) const
Read-access to permeability.
Definition: Rock_impl.hpp:84
ImmutableCMatrix PermTensor
Tensor type for read-only access to permeability.
Definition: Rock.hpp:39
Definition: BlackoilFluid.hpp:32