incorperate the review comments/decisions for multi-region PVT

the largest change is that all classes below opm/core/props/pvt take
the PVT region index as an argument, the higher-level ones (i.e.,
BlackoilProps*) take cell indices.
This commit is contained in:
Andreas Lauser
2014-05-13 12:59:11 +02:00
parent 9a7b068d15
commit 749d0e838d
13 changed files with 72 additions and 24 deletions

View File

@@ -157,9 +157,11 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesBasic::density(const int n,
const double* A,
const int* /*cells*/,
double* rho) const
{
const int np = numPhases();
@@ -177,7 +179,7 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesBasic::surfaceDensity() const
const double* BlackoilPropertiesBasic::surfaceDensity(int /*cellIdx*/) const
{
return pvt_.surfaceDensities();
}

View File

@@ -59,6 +59,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return NULL; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@@ -114,14 +119,17 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \param[in] cellIdx The index of the cell for which the surface density ought to be calculated
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.

View File

@@ -96,14 +96,19 @@ namespace Opm
void BlackoilPropertiesFromDeck::viscosity(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* mu,
double* dmudp) const
{
if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented.");
} else {
pvt_.mu(n, p, z, mu);
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
pvt_.mu(n, &pvtTableIdx[0], p, z, mu);
}
}
@@ -120,21 +125,27 @@ namespace Opm
void BlackoilPropertiesFromDeck::matrix(const int n,
const double* p,
const double* z,
const int* /*cells*/,
const int* cells,
double* A,
double* dAdp) const
{
const int np = numPhases();
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
B_.resize(n*np);
R_.resize(n*np);
if (dAdp) {
dB_.resize(n*np);
dR_.resize(n*np);
pvt_.dBdp(n, p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, p, z, &R_[0], &dR_[0]);
pvt_.dBdp(n, &pvtTableIdx[0], p, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
} else {
pvt_.B(n, p, z, &B_[0]);
pvt_.R(n, p, z, &R_[0]);
pvt_.B(n, &pvtTableIdx[0], p, z, &B_[0]);
pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
@@ -207,15 +218,19 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
void BlackoilPropertiesFromDeck::density(const int n,
const double* A,
const int* cells,
double* rho) const
{
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int cellIdx = cells?cells[i]:i;
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
const double* sdens = pvt_.surfaceDensities(pvtRegionIdx);
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
for (int comp = 0; comp < np; ++comp) {
@@ -227,9 +242,10 @@ namespace Opm
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
const double* BlackoilPropertiesFromDeck::surfaceDensity() const
const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const
{
return pvt_.surfaceDensities();
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
return pvt_.surfaceDensities(pvtRegionIdx);
}
/// \param[in] n Number of data points.

View File

@@ -97,6 +97,11 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const
{ return &cellPvtRegionIdx_[0]; }
/// \return Array of N porosity values.
virtual const double* porosity() const;
@@ -152,14 +157,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const;
virtual const double* surfaceDensity(int cellIdx = 0) const;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
@@ -206,6 +213,13 @@ namespace Opm
double* smax) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
template<class CentroidIterator>
void init(Opm::DeckConstPtr deck,
int number_of_cells,
@@ -224,6 +238,7 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock);
RockFromDeck rock_;
std::vector<int> cellPvtRegionIdx_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;
mutable std::vector<double> B_;

View File

@@ -47,6 +47,10 @@ namespace Opm
/// \return N, the number of cells.
virtual int numCells() const = 0;
/// Return an array containing the PVT table index for each
/// grid cell
virtual const int* cellPvtRegionIndex() const = 0;
/// \return Array of N porosity values.
virtual const double* porosity() const = 0;
@@ -102,14 +106,16 @@ namespace Opm
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
/// are assumed to be in Fortran order, and are typically the result
/// of a call to the method matrix().
/// \param[in] cells The index of the grid cell of each data point.
/// \param[out] rho Array of nP density values, array must be valid before calling.
virtual void density(const int n,
const double* A,
const int* cells,
double* rho) const = 0;
/// Densities of stock components at surface conditions.
/// \return Array of P density values.
virtual const double* surfaceDensity() const = 0;
virtual const double* surfaceDensity(int regionIdx = 0) const = 0;
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.

View File

@@ -95,6 +95,8 @@ namespace Opm
double* output_dRdp) const;
private:
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> density_;
std::vector<double> viscosity_;
std::vector<double> formation_volume_factor_;

View File

@@ -33,9 +33,9 @@ namespace Opm
{
}
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck )
void PvtPropertiesIncompFromDeck::init(Opm::DeckConstPtr deck)
{
// If we need multiple regions, this class and the Pvt* classes must change.
// So far, this class only supports a single PVT region. TODO?
int region_number = 0;
PhaseUsage phase_usage = phaseUsageFromDeck(deck);