Implemented matrix() method.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-01-05 21:39:33 +01:00
parent d5e7b4740c
commit 1a0e068f44
2 changed files with 57 additions and 4 deletions

View File

@ -104,13 +104,62 @@ 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
{
THROW("BlackoilPropertiesFromDeck::matrix() not yet implemented.");
}
const int np = numPhases();
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]);
} else {
pvt_.B(n, p, z, &B_[0]);
pvt_.R(n, p, z, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
pvt_.phaseUsed()[BlackoilPhases::Vapour];
const int o = phase_pos[BlackoilPhases::Liquid];
const int g = phase_pos[BlackoilPhases::Vapour];
// Compute A matrix
#pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = A + i*np*np;
std::fill(m, m + np*np, 0.0);
// Diagonal entries.
for (int phase = 0; phase < np; ++phase) {
m[phase + phase*np] = 1.0/B_[i*np + phase];
}
// Off-diagonal entries.
if (oil_and_gas) {
m[o + g*np] = R_[i*np + g]/B_[i*np + g];
m[g + o*np] = R_[i*np + o]/B_[i*np + o];
}
}
// Derivative of A matrix.
if (dAdp) {
#pragma omp parallel for
for (int i = 0; i < n; ++i) {
double* m = dAdp + i*np*np;
std::fill(m, m + np*np, 0.0);
// Diagonal entries.
for (int phase = 0; phase < np; ++phase) {
m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase];
}
// Off-diagonal entries.
if (oil_and_gas) {
m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g];
m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o];
}
}
}
}
/// \param[in] n Number of data points.
/// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the
@ -122,7 +171,7 @@ namespace Opm
const double* A,
double* rho) const
{
int np = numPhases();
const int np = numPhases();
const double* sdens = pvt_.surfaceDensities();
#pragma omp parallel for
for (int i = 0; i < n; ++i) {

View File

@ -145,6 +145,10 @@ namespace Opm
RockFromDeck rock_;
BlackoilPvtProperties pvt_;
SaturationPropsFromDeck satprops_;
mutable std::vector<double> B_;
mutable std::vector<double> dB_;
mutable std::vector<double> R_;
mutable std::vector<double> dR_;
};