Merge pull request #643 from bska/ctf-horizontal

Introduce Support for Horizontal Completions and Net-To-Gross Ratios
This commit is contained in:
Atgeirr Flø Rasmussen 2014-09-01 10:04:56 +02:00
commit 18304eb54d
3 changed files with 273 additions and 118 deletions

View File

@ -168,50 +168,127 @@ namespace WellsManagerDetail
} // namespace InjectionControl
// Compute direction permutation corresponding to completion's
// direction. First two elements of return value are directions
// perpendicular to completion while last element is direction
// along completion.
inline std::array< std::array<double,3>::size_type, 3 >
directionIndices(const Opm::CompletionDirection::DirectionEnum direction)
{
typedef std::array<double,3>::size_type idx_t;
typedef std::array<idx_t,3> permutation;
switch (direction) {
case Opm::CompletionDirection::DirectionEnum::X:
return permutation { idx_t(1), idx_t(2), idx_t(0) };
case Opm::CompletionDirection::DirectionEnum::Y:
return permutation { idx_t(2), idx_t(0), idx_t(1) };
case Opm::CompletionDirection::DirectionEnum::Z:
return permutation { idx_t(0), idx_t(1), idx_t(2) };
}
}
// Permute (diagonal) permeability components according to
// completion's direction.
inline std::array<double,3>
permComponents(const Opm::CompletionDirection::DirectionEnum direction,
const double* perm)
{
const auto p = directionIndices(direction);
const std::array<double,3>::size_type d = 3;
std::array<double,3>
K = { perm[ p[0]*(d + 1) ] ,
perm[ p[1]*(d + 1) ] ,
perm[ p[2]*(d + 1) ] };
return K;
}
// Permute cell's geometric extent according to completion's
// direction. Honour net-to-gross ratio.
//
// Note: 'extent' is intentionally accepted by modifiable value
// rather than reference-to-const to support NTG manipulation.
inline std::array<double,3>
effectiveExtent(const Opm::CompletionDirection::DirectionEnum direction,
const double ntg,
std::array<double,3> extent)
{
// Vertical extent affected by net-to-gross ratio.
extent[2] *= ntg;
const auto p = directionIndices(direction);
const std::array<double,3>::size_type d = extent.size();
std::array<double,3>
D = { extent[ p[0]*(d + 1) ] ,
extent[ p[1]*(d + 1) ] ,
extent[ p[2]*(d + 1) ] };
return D;
}
// Compute Peaceman's effective radius of single completion.
inline double
effectiveRadius(const std::array<double,3>& K,
const std::array<double,3>& D)
{
const double K01 = K[0] / K[1];
const double K10 = K[1] / K[0];
const double D0_sq = D[0] * D[0];
const double D1_sq = D[1] * D[1];
const double num = std::sqrt(K10*D0_sq + K01*D1_sq);
const double den = std::pow(K01, 0.25) + std::pow(K10, 0.25);
// Note: Analytic constant 0.28 derived for infintely sized
// formation with repeating well placement.
return 0.28 * (num / den);
}
// Use the Peaceman well model to compute well indices.
// radius is the radius of the well.
// cubical contains [dx, dy, dz] of the cell.
// (Note that the well model asumes that each cell is a cuboid).
// cell_permeability is the permeability tensor of the given cell.
// returns the well index of the cell.
double computeWellIndex(const double radius,
const std::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor)
double
computeWellIndex(const double radius,
const std::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor,
const Opm::CompletionDirection::DirectionEnum direction,
const double ntg)
{
using namespace std;
// sse: Using the Peaceman model.
// NOTE: The formula is valid for cartesian grids, so the result can be a bit
// (in worst case: there is no upper bound for the error) off the mark.
const double permx = cell_permeability[0];
const double permy = cell_permeability[3*1 + 1];
double effective_perm = sqrt(permx*permy);
// sse: The formula for r_0 can be found on page 39 of
// "Well Models for Mimetic Finite Differerence Methods and Improved Representation
// of Wells in Multiscale Methods" by Ingeborg Skjelkvåle Ligaarden.
assert(permx > 0.0);
assert(permy > 0.0);
double kxoy = permx / permy;
double kyox = permy / permx;
double r0_denominator = pow(kyox, 0.25) + pow(kxoy, 0.25);
double r0_numerator = sqrt((sqrt(kyox)*cubical[0]*cubical[0]) +
(sqrt(kxoy)*cubical[1]*cubical[1]));
assert(r0_denominator > 0.0);
double r0 = 0.28 * r0_numerator / r0_denominator;
assert(radius > 0.0);
assert(r0 > 0.0);
if (r0 < radius) {
std::cout << "ERROR: Too big well radius detected.";
std::cout << "Specified well radius is " << radius
<< " while r0 is " << r0 << ".\n";
const std::array<double,3>& K =
permComponents(direction, cell_permeability);
const std::array<double,3>& D =
effectiveExtent(direction, ntg, cubical);
const double r0 = effectiveRadius(K, D);
const double Kh = std::sqrt(K[0] * K[1]) * D[2];
// Angle of completion exposed to flow. We assume centre
// placement so there's complete exposure (= 2\pi).
const double angle =
6.2831853071795864769252867665590057683943387987502116419498;
double rw = radius;
if (r0 < rw) {
std::cerr << "Completion radius exceeds effective radius\n"
<< "Reset to effective\n";
rw = r0;
}
const long double two_pi = 6.2831853071795864769252867665590057683943387987502116419498;
double wi_denominator = log(r0 / radius) + skin_factor;
double wi_numerator = two_pi * cubical[2];
assert(wi_denominator > 0.0);
double wi = effective_perm * wi_numerator / wi_denominator;
assert(wi > 0.0);
return wi;
// NOTE: The formula is originally derived and valid for
// Cartesian grids only.
return (angle * Kh) / (std::log(r0 / rw) + skin_factor);
}
} // anonymous namespace

View File

@ -27,6 +27,8 @@
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
#include <opm/core/utility/CompressedPropertyAccess.hpp>
struct Wells;
struct UnstructuredGrid;
@ -155,7 +157,8 @@ namespace Opm
static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
template<class C2F, class CC, class FC>
template<class C2F, class CC, class FC, class NTG>
void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
const C2F& cell_to_faces,
const int* cart_dims,
@ -166,8 +169,9 @@ namespace Opm
std::vector<WellData>& well_data,
std::map<std::string, int> & well_names_to_index,
const PhaseUsage& phaseUsage,
const std::map<int,int> cartesian_to_compressed,
const double* permeability);
const std::map<int,int>& cartesian_to_compressed,
const double* permeability,
const NTG& ntg);
void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage);
void setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index);

View File

@ -1,7 +1,13 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <array>
#include <cstddef>
#include <exception>
#include <iterator>
namespace WellsManagerDetail
{
@ -44,39 +50,57 @@ namespace WellsManagerDetail
double computeWellIndex(const double radius,
const std::array<double, 3>& cubical,
const double* cell_permeability,
const double skin_factor);
const double skin_factor,
const Opm::CompletionDirection::DirectionEnum direction,
const double ntg);
template<class C2F, class FC>
std::array<double, 3> getCubeDim(const C2F& c2f, FC begin_face_centroids, int dimensions,
int cell)
template <int dim, class C2F, class FC>
std::array<double, dim>
getCubeDim(const C2F& c2f,
FC begin_face_centroids,
int cell)
{
using namespace std;
std::array<double, 3> cube;
//typedef Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type Cell2Faces;
//Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid);
typedef typename C2F::row_type FaceRow;
FaceRow faces=c2f[cell];
typename FaceRow::const_iterator face=faces.begin();
int num_local_faces = faces.end()-face;
vector<double> x(num_local_faces);
vector<double> y(num_local_faces);
vector<double> z(num_local_faces);
for (int lf=0; lf<num_local_faces; ++ lf, ++face) {
FC centroid = Opm::UgGridHelpers::increment(begin_face_centroids, *face, dimensions);
x[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 0);
y[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 1);
z[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 2);
std::array< std::vector<double>, dim > X;
{
const std::vector<double>::size_type
nf = std::distance(c2f[cell].begin(),
c2f[cell].end ());
for (int d = 0; d < dim; ++d) {
X[d].reserve(nf);
}
}
cube[0] = *max_element(x.begin(), x.end()) - *min_element(x.begin(), x.end());
cube[1] = *max_element(y.begin(), y.end()) - *min_element(y.begin(), y.end());
cube[2] = *max_element(z.begin(), z.end()) - *min_element(z.begin(), z.end());
typedef typename C2F::row_type::const_iterator FI;
for (FI f = c2f[cell].begin(), e = c2f[cell].end(); f != e; ++f) {
using Opm::UgGridHelpers::increment;
using Opm::UgGridHelpers::getCoordinate;
const FC& fc = increment(begin_face_centroids, *f, dim);
for (int d = 0; d < dim; ++d) {
X[d].push_back(getCoordinate(fc, d));
}
}
std::array<double, dim> cube;
for (int d = 0; d < dim; ++d) {
typedef std::vector<double>::iterator VI;
typedef std::pair<VI,VI> PVI;
const PVI m = std::minmax_element(X[d].begin(), X[d].end());
cube[d] = *m.second - *m.first;
}
return cube;
}
} // end namespace
} // end namespace WellsManagerDetail
namespace Opm
{
template<class C2F, class CC, class FC>
template<class C2F, class CC, class FC, class NTG>
void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t timeStep,
const C2F& c2f,
const int* cart_dims,
@ -87,9 +111,16 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
std::vector<WellData>& well_data,
std::map<std::string, int>& well_names_to_index,
const PhaseUsage& phaseUsage,
std::map<int,int> cartesian_to_compressed,
const double* permeability)
const std::map<int,int>& cartesian_to_compressed,
const double* permeability,
const NTG& ntg)
{
if (dimensions != 3) {
OPM_THROW(std::domain_error,
"WellsManager::createWellsFromSpecs() only "
"supported in three space dimensions");
}
std::vector<std::vector<PerfData> > wellperf_data;
wellperf_data.resize(wells.size());
@ -141,10 +172,16 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
std::array<double, 3> cubical = WellsManagerDetail::getCubeDim(c2f, begin_face_centroids,
dimensions, cell);
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(c2f, begin_face_centroids, cell);
const double* cell_perm = &permeability[dimensions*dimensions*cell];
pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, completion->getSkinFactor());
pd.well_index =
WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
completion->getSkinFactor(),
completion->getDirection(),
ntg[cell]);
}
wellperf_data[well_index].push_back(pd);
}
@ -157,7 +194,7 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
const int num_wells = well_data.size();
int num_perfs = 0;
assert(dimensions == 3);
assert (dimensions == 3);
for (int w = 0; w < num_wells; ++w) {
num_perfs += wellperf_data[w].size();
if (well_data[w].reference_bhp_depth == -1e100) {
@ -165,13 +202,19 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
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);
using UgGridHelpers::increment;
using UgGridHelpers::getCoordinate;
const CC& cc =
increment(begin_cell_centroids,
wellperf_data[w][perf].cell,
dimensions);
const double depth = getCoordinate(cc, 2);
min_depth = std::min(min_depth, depth);
}
well_data[w].reference_bhp_depth = min_depth;
}
}
@ -185,66 +228,86 @@ void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t
// Add wells.
for (int w = 0; w < num_wells; ++w) {
const int w_num_perf = wellperf_data[w].size();
std::vector<int> perf_cells(w_num_perf);
const int w_num_perf = wellperf_data[w].size();
std::vector<int> perf_cells (w_num_perf);
std::vector<double> perf_prodind(w_num_perf);
for (int perf = 0; perf < w_num_perf; ++perf) {
perf_cells[perf] = wellperf_data[w][perf].cell;
perf_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_);
const 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.");
OPM_THROW(std::runtime_error,
"Failed adding well "
<< well_names[w]
<< " to Wells data structure.");
}
}
}
template<class CC, class C2F, class FC>
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& cell_to_faces,
FC begin_face_centroids,
const double* permeability)
template <class CC, class C2F, class FC>
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& cell_to_faces,
FC begin_face_centroids,
const double* permeability)
: w_(0)
{
init(eclipseState, timeStep, number_of_cells, global_cell, cart_dims, dimensions,
begin_cell_centroids, cell_to_faces, begin_face_centroids, permeability);
init(eclipseState, timeStep, number_of_cells, global_cell,
cart_dims, dimensions, begin_cell_centroids,
cell_to_faces, begin_face_centroids, permeability);
}
/// Construct wells from deck.
template<class CC, class C2F, class FC>
void WellsManager::init(const Opm::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& cell_to_faces,
FC begin_face_centroids,
const double* permeability)
template <class CC, class C2F, class FC>
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& 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.");
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");
OPM_MESSAGE("No wells specified in Schedule section, "
"initializing no wells");
return;
}
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(global_cell,
number_of_cells, cartesian_to_compressed);
setupCompressedToCartesian(global_cell, number_of_cells,
cartesian_to_compressed);
// Obtain phase usage data.
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
@ -258,35 +321,46 @@ void WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
// For easy lookup:
std::map<std::string, int> well_names_to_index;
ScheduleConstPtr schedule = eclipseState->getSchedule();
std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
ScheduleConstPtr schedule = eclipseState->getSchedule();
std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
well_names.reserve(wells.size());
well_data.reserve(wells.size());
typedef GridPropertyAccess::ArrayPolicy::ExtractFromDeck<double> DoubleArray;
typedef GridPropertyAccess::Compressed<DoubleArray, GridPropertyAccess::Tag::NTG> NTGArray;
DoubleArray ntg_glob(eclipseState, "NTG", 1.0);
NTGArray ntg(ntg_glob, global_cell);
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);
well_names, well_data, well_names_to_index,
pu, cartesian_to_compressed, permeability, ntg);
setupWellControls(wells, timeStep, well_names, pu);
{
GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
GroupTreeNodeConstPtr fieldNode =
schedule->getGroupTree(timeStep)->getNode("FIELD");
GroupConstPtr fieldGroup =
schedule->getGroup(fieldNode->name());
well_collection_.addField(fieldGroup, timeStep, pu);
addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
addChildGroups(fieldNode, schedule, timeStep, pu);
}
for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
well_collection_.addWell((*wellIter), timeStep, pu);
for (auto w = wells.begin(), e = wells.end(); w != e; ++w) {
well_collection_.addWell(*w, timeStep, pu);
}
well_collection_.setWellsPointer(w_);
well_collection_.applyGroupControls();
setupGuideRates(wells, timeStep, well_data, well_names_to_index);
// Debug output.