OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-18 21:43:27 -06:00
Code Issues Packages Projects Releases Wiki Activity

matrix(): Honour chain rule of differentiation.

Browse Source 
The original implementation of change-set 90d8dd8c8040 contained a
crucial misprint leading to incorrect results for all compressible
fluids.
This commit is contained in:
Bård Skaflestad 2012-05-11 14:26:48 +02:00
parent 0677f58026
commit 29cb30ed58
1 changed files with 29 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -147,19 +147,40 @@ namespace Opm
} }
// Derivative of A matrix. // Derivative of A matrix.
// A = R*inv(B) whence
//
// dA/dp = (dR/dp*inv(B) + R*d(inv(B))/dp)
// = (dR/dp*inv(B) - R*inv(B)*(dB/dp)*inv(B))
// = (dR/dp - A*(dB/dp)) * inv(B)
//
// The B matrix is diagonal and that fact is exploited in the
// following implementation.
if (dAdp) { if (dAdp) {
// #pragma omp parallel for // #pragma omp parallel for
std::copy(A, A + n*np*np, dAdp);
for (int i = 0; i < n; ++i) { for (int i = 0; i < n; ++i) {
double* m = dAdp + i*np*np; double* m = dAdp + i*np*np;
std::fill(m, m + np*np, 0.0);
// Diagonal entries. const double* dB = & dB_[i * np];
for (int phase = 0; phase < np; ++phase) { for (int i2 = 0; i2 < np; ++i2) {
m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase]; for (int i1 = 0; i1 < np; ++i1) {
m[i2*np + i1] *= - dB[ i1 ]; // Note sign.
}
} }
// Off-diagonal entries.
if (oil_and_gas) { if (oil_and_gas) {
m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g]; const double* dR = & dR_[i * np];
m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o];
m[o*np + g] += dR[ o ];
m[g*np + o] += dR[ g ];
}
const double* B = & B_[i * np];
for (int i2 = 0; i2 < np; ++i2) {
for (int i1 = 0; i1 < np; ++i1) {
m[i2*np + i1] /= B[ i1 ];
}
} }
} }
} }