Use fastSparseProduct() in place of Eigen's product.

This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-05 12:59:04 +01:00
parent c43b9f4a22
commit 5b6765f9d8
3 changed files with 38 additions and 17 deletions

View File

@ -392,7 +392,7 @@ namespace Opm
const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmudp_diag * pw.derivative()[block];
fastSparseProduct(dmudp_diag, pw.derivative()[block], jacs[block]);
}
return ADB::function(mu, jacs);
}
@ -427,7 +427,10 @@ namespace Opm
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmudp_diag * po.derivative()[block] + dmudr_diag * rs.derivative()[block];
fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
jacs[block] += temp;
}
return ADB::function(mu, jacs);
}
@ -458,7 +461,7 @@ namespace Opm
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmudp_diag * pg.derivative()[block];
fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]);
}
return ADB::function(mu, jacs);
}
@ -493,7 +496,10 @@ namespace Opm
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmudp_diag * pg.derivative()[block] + dmudr_diag * rv.derivative()[block];
fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
jacs[block] += temp;
}
return ADB::function(mu, jacs);
}
@ -653,7 +659,7 @@ namespace Opm
const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dbdp_diag * pw.derivative()[block];
fastSparseProduct(dbdp_diag, pw.derivative()[block], jacs[block]);
}
return ADB::function(b, jacs);
}
@ -689,7 +695,10 @@ namespace Opm
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dbdp_diag * po.derivative()[block] + dbdr_diag * rs.derivative()[block];
fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
jacs[block] += temp;
}
return ADB::function(b, jacs);
}
@ -721,7 +730,7 @@ namespace Opm
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dbdp_diag * pg.derivative()[block];
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
}
return ADB::function(b, jacs);
}
@ -753,11 +762,14 @@ namespace Opm
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
ADB::M dmudr_diag = spdiag(dbdr);
ADB::M dbdr_diag = spdiag(dbdr);
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dbdp_diag * pg.derivative()[block] + dmudr_diag * rv.derivative()[block];;
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
jacs[block] += temp;
}
return ADB::function(b, jacs);
}
@ -817,7 +829,7 @@ namespace Opm
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = drbubdp_diag * po.derivative()[block];
fastSparseProduct(drbubdp_diag, po.derivative()[block], jacs[block]);
}
return ADB::function(rbub, jacs);
}
@ -889,7 +901,7 @@ namespace Opm
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = drvdp_diag * po.derivative()[block];
fastSparseProduct(drvdp_diag, po.derivative()[block], jacs[block]);
}
return ADB::function(rv, jacs);
}
@ -1004,7 +1016,9 @@ namespace Opm
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm()
ADB::M dkr1_ds2_diag = spdiag(dkr.col(column));
for (int block = 0; block < num_blocks; ++block) {
jacs[block] += dkr1_ds2_diag * s[phase2]->derivative()[block];
ADB::M temp;
fastSparseProduct(dkr1_ds2_diag, s[phase2]->derivative()[block], temp);
jacs[block] += temp;
}
}
relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs));
@ -1062,7 +1076,9 @@ namespace Opm
const int column = phase1_pos + numActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm()
ADB::M dpc1_ds2_diag = spdiag(dpc.col(column));
for (int block = 0; block < numBlocks; ++block) {
jacs[block] += dpc1_ds2_diag * s[phase2]->derivative()[block];
ADB::M temp;
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp);
jacs[block] += temp;
}
}
adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs));

View File

@ -2059,7 +2059,7 @@ namespace {
const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dpm_diag * p.derivative()[block];
fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(pm, jacs);
} else {
@ -2087,7 +2087,7 @@ namespace {
const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dtm_diag * p.derivative()[block];
fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(tm, jacs);
} else {

View File

@ -279,7 +279,9 @@ namespace Opm
continue;
}
// solve Du = C
const M u = Di * Jn[var]; // solver.solve(Jn[var]);
// const M u = Di * Jn[var]; // solver.solve(Jn[var]);
M u;
fastSparseProduct(Di, Jn[var], u); // solver.solve(Jn[var]);
for (int eq = 0; eq < num_eq; ++eq) {
if (eq == n) {
continue;
@ -292,7 +294,9 @@ namespace Opm
jacs[eq].push_back(Je[var]);
M& J = jacs[eq].back();
// Subtract Bu (B*inv(D)*C)
J -= B * u;
M Bu;
fastSparseProduct(B, u, Bu);
J -= Bu;
}
}
@ -397,6 +401,7 @@ namespace Opm
void formEllipticSystem(const int num_phases,
const std::vector<ADB>& eqs_in,
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
// M& A,
V& b)
{
if (num_phases != 3) {