fix eigen bug when compute the capPress.

This commit is contained in:
Liu Ming 2014-01-27 16:53:45 +08:00
parent 075e16dc36
commit f01c4dac10

View File

@ -137,42 +137,42 @@ namespace Opm
const ADB& so, const ADB& so,
const Cells& cells) const const Cells& cells) const
{ {
const int numCells = cells.size();
const int n = cells.size(); const int numActivePhases = numPhases();
const int np = numPhases(); const int numBlocks = so.numBlocks();
const int num_blocks = so.numBlocks(); assert(sw.value().size() == numCells);
Block s_all(n, np); assert(so.value().size() == numCells);
assert(sw.size() == n && so.size() == n); Block s_all(numCells, numActivePhases);
s_all.col(0) = sw.value(); s_all.col(0) = sw.value();
s_all.col(1) = so.value(); s_all.col(1) = so.value();
Block pc(n, np); Block pc(numCells, numActivePhases);
Block dpc(n, np * np); Block dpc(numCells, numActivePhases*numActivePhases);
satprops_.capPress(numCells, s_all.data(), cells.data(), pc.data(), dpc.data());
satprops_.capPress(n, s_all.data(), cells.data(), pc.data(), dpc.data()); std::vector<ADB> adbCapPressures;
adbCapPressures.reserve(2);
std::vector<ADB> capPressures;
capPressures.reserve(2);
const ADB* s[2] = { &sw, &so}; const ADB* s[2] = { &sw, &so};
for (int phase1 = 0; phase1 < 3; ++phase1) { for (int phase1 = 0; phase1 < 2; ++phase1) {
const int phase1_pos = phase1; const int phase1_pos = phase1;
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(numBlocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < numBlocks; ++block) {
jacs[block] = ADB::M(n, s[phase1]->derivative()[block].cols()); jacs[block] = ADB::M(numCells, s[phase1]->derivative()[block].cols());
} }
for (int phase2 = 0; phase2 < 3; ++phase2) { for (int phase2 = 0; phase2 < 2; ++phase2) {
const int phase2_pos = phase2; const int phase2_pos = phase2;
// Assemble dpc1/ds2. // Assemble dpc1/ds2.
const int column = phase1_pos + phase2_pos; // Recall: Fortran ordering from props_.relperm() const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column)); ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < numBlocks; ++block) {
jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block]; jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
} }
} }
capPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs)); adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));
} }
return capPressures; return adbCapPressures;
} }
} //namespace Opm } //namespace Opm