diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 81a3f1c66..a03e3cc53 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -253,6 +253,24 @@ namespace Opm return; } + // global_cell is a map from compressed cells to Cartesian grid cells. + // We must make the inverse lookup. + const int* global_cell = grid.global_cell; + const int* cpgdim = grid.cartdims; + std::map cartesian_to_compressed; + + if (global_cell) { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + // Obtain phase usage data. PhaseUsage pu = phaseUsageFromDeck(eclipseState); @@ -273,123 +291,146 @@ namespace Opm std::vector wells = schedule->getWells(timeStep); well_names.reserve(wells.size()); well_data.reserve(wells.size()); + wellperf_data.resize(wells.size()); int well_index = 0; for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { WellConstPtr well = (*wellIter); - well_names_to_index[well->name()] = well_index; - well_names.push_back(well->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 = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); - wd.welspecsline = -1; - well_data.push_back(wd); + { // WELSPECS handling + well_names_to_index[well->name()] = well_index; + well_names.push_back(well->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 = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); + wd.welspecsline = -1; + well_data.push_back(wd); + } + } + + { // COMPDAT handling + CompletionSetConstPtr completionSet = well->getCompletions(timeStep); + for (size_t c=0; csize(); c++) { + CompletionConstPtr completion = completionSet->get(c); + int i = completion->getI() - 1; + int j = completion->getJ() - 1; + int k = completion->getK() - 1; + 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()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + } + int cell = cgit->second; + PerfData pd; + pd.cell = cell; + if (completion->getCF() > 0.0) { + pd.well_index = completion->getCF(); + } else { + double radius = 0.5*completion->getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + std::array cubical = getCubeDim(grid, cell); + const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; + pd.well_index = computeWellIndex(radius, cubical, cell_perm, + completion->getDiameter()); + } + wellperf_data[well_index].push_back(pd); + } } well_index++; } - const int num_wells = well_data.size(); - wellperf_data.resize(num_wells); - // global_cell is a map from compressed cells to Cartesian grid cells. - // We must make the inverse lookup. - const int* global_cell = grid.global_cell; - const int* cpgdim = grid.cartdims; - std::map cartesian_to_compressed; - if (global_cell) { - for (int i = 0; i < grid.number_of_cells; ++i) { - cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); - } - } - else { - for (int i = 0; i < grid.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 - // 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; - WellConstPtr well = schedule->getWell(well_names[wix]); - if (ix < 0) { - // Defaulted I location. Extract from WELSPECS. - ix = well->getHeadI() - 1; - } - if (jy < 0) { - // Defaulted J location. Extract from WELSPECS. - jy = well->getHeadJ() - 1; - } - if (kz1 < 0) { - // Defaulted KZ1. Use top layer. - kz1 = 0; - } - if (kz2 < 0) { - // Defaulted KZ2. Use bottom layer. - kz2 = cpgdim[2] - 1; - } +// // 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 +// // 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; - 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()) { - OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' ' - << kz << " not found in grid (well = " << name << ')'); - } - 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 = getCubeDim(grid, cell); - const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; - pd.well_index = 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"); - } - } +// WellConstPtr well = schedule->getWell(well_names[wix]); +// if (ix < 0) { +// // Defaulted I location. Extract from WELSPECS. +// ix = well->getHeadI() - 1; +// } +// if (jy < 0) { +// // Defaulted J location. Extract from WELSPECS. +// jy = well->getHeadJ() - 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()) { +// OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' ' +// << kz << " not found in grid (well = " << name << ')'); +// } +// 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 = getCubeDim(grid, cell); +// const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; +// pd.well_index = 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"); +// } +// } // Set up reference depths that were defaulted. Count perfs. + + const int num_wells = well_data.size(); + int num_perfs = 0; assert(grid.dimensions == 3); for (int w = 0; w < num_wells; ++w) { diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index f14be7946..145fb2e25 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -124,7 +124,7 @@ namespace Opm // Data Wells* w_; - WellCollection well_collection_; + WellCollection well_collection_; diff --git a/tests/wells_manager_data_expanded.data b/tests/wells_manager_data_expanded.data index dfcd0b48e..3d97003ae 100755 --- a/tests/wells_manager_data_expanded.data +++ b/tests/wells_manager_data_expanded.data @@ -26,7 +26,7 @@ WELSPECS COMPDAT 'INJ1' 1 1 1 2 'OPEN' 1 10.6092 0.5 / 'INJ1' 1 1 3 5 'OPEN' 1 12.6092 0.5 / - 'INJ1' 1 1 2* 'OPEN' 1 14.6092 0.5 / + 'INJ1' 1 1 1 1 'OPEN' 1 14.6092 0.5 / 'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 / /