OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #433 from GitPaean/optimizing_interleaved_implementation

Browse Source 
optimizing the formInterleavedSystem().
This commit is contained in:
Atgeirr Flø Rasmussen 2015-08-17 16:06:02 +02:00
commit dd87f45dfb
2 changed files with 21 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -105,14 +105,26 @@ namespace Opm
eqs[phase] = eqs[phase] * matbalscale[phase]; eqs[phase] = eqs[phase] * matbalscale[phase];
} }
// Form modified system. // calculating the size for b
Eigen::SparseMatrix<double, Eigen::RowMajor> A; int size_b = 0;
V b; for (int elem = 0; elem < np; ++elem) {
formEllipticSystem(np, eqs, A, b); const int loc_size = eqs[elem].size();
size_b += loc_size;
}
V b(size_b);
int pos = 0;
for (int elem = 0; elem < np; ++elem) {
const int loc_size = eqs[elem].size();
b.segment(pos, loc_size) = eqs[elem].value();
pos += loc_size;
}
assert(pos == size);
// Create ISTL matrix with interleaved rows and columns (block structured). // Create ISTL matrix with interleaved rows and columns (block structured).
Mat istlA; Mat istlA;
formInterleavedSystem(eqs, A, istlA); formInterleavedSystem(eqs, istlA);
// Solve reduced system. // Solve reduced system.
SolutionVector dx(SolutionVector::Zero(b.size())); SolutionVector dx(SolutionVector::Zero(b.size()));
@ -185,7 +197,6 @@ namespace Opm
void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs, void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs,
const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
Mat& istlA) const Mat& istlA) const
{ {
const int np = eqs.size(); const int np = eqs.size();
@ -194,9 +205,10 @@ namespace Opm
// corresponding to the jacobians with respect to pressure. // corresponding to the jacobians with respect to pressure.
// Use addition to get to the union structure. // Use addition to get to the union structure.
Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0]; Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0];
for (int phase = 0; phase < np; ++phase) { for (int phase = 1; phase < np; ++phase) {
structure += eqs[phase].derivative()[0]; structure += eqs[phase].derivative()[0];
} }
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure; Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;
// Create ISTL matrix with interleaved rows and columns (block structured). // Create ISTL matrix with interleaved rows and columns (block structured).
@ -211,17 +223,15 @@ namespace Opm
row.insert(ja[i]); row.insert(ja[i]);
} }
} }
const int size = s.rows(); const int size = s.rows();
Span span[3] = { Span(size, 1, 0),
Span(size, 1, size),
Span(size, 1, 2*size) };
for (int row = 0; row < size; ++row) { for (int row = 0; row < size; ++row) {
for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) { for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
const int col = ja[col_ix]; const int col = ja[col_ix];
MatrixBlockType block; MatrixBlockType block;
for (int p1 = 0; p1 < np; ++p1) { for (int p1 = 0; p1 < np; ++p1) {
for (int p2 = 0; p2 < np; ++p2) { for (int p2 = 0; p2 < np; ++p2) {
block[p1][p2] = A.coeff(span[p1][row], span[p2][col]); block[p1][p2] = eqs[p1].derivative()[p2].coeff(row, col);
} }
} }
istlA[row][col] = block; istlA[row][col] = block;

1
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@ -156,7 +156,6 @@ namespace Opm
} }
void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs, void formInterleavedSystem(const std::vector<LinearisedBlackoilResidual::ADB>& eqs,
const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
Mat& istlA) const; Mat& istlA) const;
mutable int iterations_; mutable int iterations_;