adding computeRepRadiusPerfLength() to WellInterface.

This commit is contained in:
Kai Bao 2017-07-31 16:42:26 +02:00
parent d1727d0183
commit d4fa8c06f1
3 changed files with 98 additions and 95 deletions

View File

@ -964,106 +964,13 @@ namespace Opm {
// to be fixed later.
int number_of_cells = Opm::UgGridHelpers::numCells(grid);
const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
if (wells_ecl_.size() == 0) {
OPM_MESSAGE("No wells specified in Schedule section, "
"initializing no wells");
return;
}
const int nw = number_of_wells_;
const int nperf = wells().well_connpos[nw];
const size_t timeStep = current_timeIdx_;
wells_rep_radius_.clear();
wells_perf_length_.clear();
wells_bore_diameter_.clear();
wells_rep_radius_.reserve(nperf);
wells_perf_length_.reserve(nperf);
wells_bore_diameter_.reserve(nperf);
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(global_cell, number_of_cells,
cartesian_to_compressed);
int well_index = 0;
for (auto wellIter= wells_ecl_.begin(); wellIter != wells_ecl_.end(); ++wellIter) {
const auto* well = (*wellIter);
if (well->getStatus(timeStep) == WellCommon::SHUT) {
continue;
}
{ // COMPDAT handling
const auto& completionSet = well->getCompletions(timeStep);
for (size_t c=0; c<completionSet.size(); c++) {
const auto& completion = completionSet.get(c);
if (completion.getState() == WellCompletion::OPEN) {
int i = completion.getI();
int j = completion.getJ();
int k = completion.getK();
const int* cpgdim = cart_dims;
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()) {
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;
{
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);
}
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
WellCompletion::DirectionEnum direction = completion.getDirection();
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
break;
case Opm::WellCompletion::DirectionEnum::Y:
re = std::sqrt(cubical[0] * cubical[2] / M_PI);
perf_length = cubical[1];
break;
case Opm::WellCompletion::DirectionEnum::Z:
re = std::sqrt(cubical[0] * cubical[1] / M_PI);
perf_length = cubical[2];
break;
default:
OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
}
double repR = std::sqrt(re * radius);
wells_rep_radius_.push_back(repR);
wells_perf_length_.push_back(perf_length);
wells_bore_diameter_.push_back(2. * radius);
}
} else {
if (completion.getState() != WellCompletion::SHUT) {
OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
}
}
}
}
well_index++;
for (const auto& well : well_container_) {
well->computeRepRadiusPerfLength(grid, cartesian_to_compressed);
}
}

View File

@ -63,6 +63,7 @@ namespace Opm
using WellState = WellStateFullyImplicitBlackoilDense;
typedef BlackoilModelParameters ModelParameters;
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
@ -213,6 +214,9 @@ namespace Opm
void updateListEconLimited(const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
void computeRepRadiusPerfLength(const Grid& grid, const std::map<int, int>& cartesian_to_compressed);
protected:
// to indicate a invalid connection
@ -272,6 +276,15 @@ namespace Opm
// saturation table nubmer for each well perforation
std::vector<int> saturation_table_number_;
// representative radius of the perforations, used in shear calculation
std::vector<double> perf_rep_radius_;
// length of the perforations, use in shear calculation
std::vector<double> perf_length_;
// well bore diameter
std::vector<double> bore_diameters_;
const PhaseUsage* phase_usage_;
const std::vector<bool>* active_;

View File

@ -779,4 +779,87 @@ namespace Opm
}
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
computeRepRadiusPerfLength(const Grid& grid,
const std::map<int, int>& cartesian_to_compressed)
{
const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
const int nperf = numberOfPerforations();
perf_rep_radius_.clear();
perf_length_.clear();
bore_diameters_.clear();
perf_rep_radius_.resize(nperf);
perf_length_.resize(nperf);
bore_diameters_.resize(nperf);
// COMPDAT handling
const auto& completionSet = well_ecl_->getCompletions(current_step_);
for (size_t c=0; c<completionSet.size(); c++) {
const auto& completion = completionSet.get(c);
if (completion.getState() == WellCompletion::OPEN) {
const int i = completion.getI();
const int j = completion.getJ();
const int k = completion.getK();
const int* cpgdim = cart_dims;
const int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
const std::map<int, int>::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 = " << name() << ')');
}
const int cell = cgit->second;
{
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);
}
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
WellCompletion::DirectionEnum direction = completion.getDirection();
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
break;
case Opm::WellCompletion::DirectionEnum::Y:
re = std::sqrt(cubical[0] * cubical[2] / M_PI);
perf_length = cubical[1];
break;
case Opm::WellCompletion::DirectionEnum::Z:
re = std::sqrt(cubical[0] * cubical[1] / M_PI);
perf_length = cubical[2];
break;
default:
OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
}
const double repR = std::sqrt(re * radius);
perf_rep_radius_.push_back(repR);
perf_length_.push_back(perf_length);
bore_diameters_.push_back(2. * radius);
}
}
}
}
}