Improve well matrix products

This commit is contained in:
Markus Blatt
2017-05-17 15:54:07 +02:00
parent 22dfa6f135
commit 01b90bce05

View File

@@ -2105,18 +2105,23 @@ namespace Detail
void void
StandardWell<TypeTag>::addWellContributions(Mat& mat) const StandardWell<TypeTag>::addWellContributions(Mat& mat) const
{ {
auto colC = duneC_[0].begin(); // We need to change matrx A as follows
for(int perf1 = 0; perf1 < number_of_perforations_; ++perf1, ++colC) { // A -= C^T D^-1 B
const auto row_index = well_cells_[perf1]; // D is diagonal
// B and C have 1 row, nc colums and nonzero
// at (i,j) only if well i has perforation at cell j.
for(auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) {
const auto row_index = colC.index();
auto& row = mat[row_index]; auto& row = mat[row_index];
auto colB = duneB_[0].begin(); auto colB = duneB_[0].begin();
assert(colC.index() == row_index); auto col = row.begin();
for(int perf2 = 0 ; perf2 < number_of_perforations_; ++perf2, ++colB) { for(auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) {
const auto col_index = well_cells_[perf2]; const auto col_index = colB.index();
auto col = row.find(col_index); // Move col to index pj
assert(col != row.end()); while(col != row.end() && col.index() < col_index) ++col;
assert(colB.index() == col_index); assert(col != row.end() && col.index() == col_index);
Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp; Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp;
typename Mat::block_type tmp1; typename Mat::block_type tmp1;