OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implemented matrix() method.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-01-05 21:39:33 +01:00
parent f278f87f14
commit 0c15624617
3 changed files with 66 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
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) {

4
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
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_;
};

11
tests/bo_resprop_test.cpp
View File

@ -22,6 +22,7 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/eclipse/EclipseGridInspector.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <iterator>
@ -35,8 +36,14 @@ int main(int argc, char** argv)
// Parser.
std::string ecl_file = param.get<std::string>("filename");
Dune::EclipseGridParser deck(ecl_file);
Opm::BlackoilPropertiesFromDeck props(deck);
Dune::EclipseGridInspector insp(deck);
std::tr1::array<int, 3> gs = insp.gridSize();
int num_cells = gs[0]*gs[1]*gs[2];
std::vector<int> global_cell(num_cells);
for (int i = 0; i < num_cells; ++i) {
global_cell[i] = i;
}
Opm::BlackoilPropertiesFromDeck props(deck, global_cell);
const int n = 1;
double p[n] = { 150e5 };