diff --git a/flowexperimental/comp/wells/CompWellEquations_impl.hpp b/flowexperimental/comp/wells/CompWellEquations_impl.hpp index fe1190867f..484ce60a26 100644 --- a/flowexperimental/comp/wells/CompWellEquations_impl.hpp +++ b/flowexperimental/comp/wells/CompWellEquations_impl.hpp @@ -25,9 +25,6 @@ template CompWellEquations:: CompWellEquations() { - duneB_.setBuildMode(OffDiagMatWell::row_wise); - duneC_.setBuildMode(OffDiagMatWell::row_wise), - duneD_.setBuildMode(DiagMatWell::row_wise); invDuneD_.setBuildMode(DiagMatWell::row_wise); } @@ -37,18 +34,22 @@ void CompWellEquations:: init(const int num_conn, const std::vector& cells) { + duneB_.setBuildMode(OffDiagMatWell::row_wise); + duneC_.setBuildMode(OffDiagMatWell::row_wise), + duneD_.setBuildMode(DiagMatWell::row_wise); + duneD_.setSize(1, 1, 1); duneB_.setSize(1, num_conn, num_conn); duneC_.setSize(1, num_conn, num_conn); - for (auto row = duneD_.createbegin(), - end = duneD_.createend(); row != end; ++row) { + auto endD = duneD_.createend(); + for (auto row = duneD_.createbegin(); row != endD; ++row) { // Add nonzeros for diagonal row.insert(row.index()); } - for (auto row = duneB_.createbegin(), - end = duneB_.createend(); row != end; ++row) { + auto endB = duneB_.createend(); + for (auto row = duneB_.createbegin(); row != endB; ++row) { for (int con = 0 ; con < num_conn; ++con) { row.insert(con); } @@ -56,8 +57,8 @@ init(const int num_conn, const std::vector& cells) // make the C^T matrix // TODO: let us see whether we should change the naming of DuneC_ - for (auto row = duneC_.createbegin(), - end = duneC_.createend(); row != end; ++row) { + auto endC = duneC_.createend(); + for (auto row = duneC_.createbegin(); row != endC; ++row) { for (int con = 0; con < num_conn; ++con) { row.insert(con); } diff --git a/opm/simulators/linalg/DILU.hpp b/opm/simulators/linalg/DILU.hpp index 9bd1abe241..7e4ce5f7a6 100644 --- a/opm/simulators/linalg/DILU.hpp +++ b/opm/simulators/linalg/DILU.hpp @@ -283,15 +283,16 @@ private: OPM_TIMEBLOCK(upper_solve); // upper triangular solve: (D + U_A) v = Dy - auto rendi = A_.beforeBegin(); - for (auto row = A_.beforeEnd(); row != rendi; --row) { - const auto row_i = row.index(); + auto rindex = [](auto it) { return std::prev(it.base()).index(); }; + auto rbegini = std::make_reverse_iterator(A_.begin()); + for (auto row = std::make_reverse_iterator(A_.end()); row != rbegini; ++row) { + const auto row_i = rindex(row); // rhs = 0 Xblock rhs(0.0); - for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij) { + for (auto a_ij = std::make_reverse_iterator(row->end()); rindex(a_ij) > row_i; ++a_ij) { // if A[i][j] != 0 // rhs += A[i][j]*v[j] - const auto col_j = a_ij.index(); + const auto col_j = rindex(a_ij); a_ij->umv(v[col_j], rhs); } // calculate update v = M^-1*d @@ -346,9 +347,10 @@ private: auto row = A_reordered_->begin() + level_start_idx + row_idx_in_level; const auto row_i = reordered_to_natural_[row.index()]; Xblock rhs(0.0); - for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij) { + auto rindex = [](auto it) { return std::prev(it.base()).index(); }; + for (auto a_ij = std::make_reverse_iterator((*row).end()); rindex(a_ij) > row_i; ++a_ij) { // rhs += A[i][j]*v[j] - const auto col_j = a_ij.index(); + const auto col_j = rindex(a_ij); a_ij->umv(v[col_j], rhs); } // calculate update v = M^-1*d diff --git a/opm/simulators/linalg/GraphColoring.hpp b/opm/simulators/linalg/GraphColoring.hpp index cc85fc6897..0e2724b070 100644 --- a/opm/simulators/linalg/GraphColoring.hpp +++ b/opm/simulators/linalg/GraphColoring.hpp @@ -273,14 +273,15 @@ getMatrixRowColoring(const M& matrix, ColoringType coloringType) } } } else if (coloringType == ColoringType::UPPER) { - for (auto i = matrix.beforeEnd(); i != matrix.beforeBegin(); --i) { - for (auto a_ij = ++(i->find(i.index())); a_ij != i->end(); ++a_ij) { - color[i.index()] = std::max({color[i.index()], color[a_ij.index()] + 1}); + auto rindex = [](auto it) { return std::prev(it.base()).index(); }; + for (auto i = std::make_reverse_iterator(matrix.end()); i != std::make_reverse_iterator(matrix.begin()); ++i) { + for (auto a_ij = ++(i->find(rindex(i))); a_ij != i->end(); ++a_ij) { + color[rindex(i)] = std::max({color[rindex(i)], color[a_ij.index()] + 1}); } - if (color[i.index()] >= colorCnt.size()) { + if (color[rindex(i)] >= colorCnt.size()) { colorCnt.push_back(1); } else { - ++colorCnt[color[i.index()]]; + ++colorCnt[color[rindex(i)]]; } } } else if (coloringType == ColoringType::LOWER) { diff --git a/opm/simulators/linalg/MILU.cpp b/opm/simulators/linalg/MILU.cpp index f793a9a645..c1a32556f6 100644 --- a/opm/simulators/linalg/MILU.cpp +++ b/opm/simulators/linalg/MILU.cpp @@ -167,23 +167,21 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU, auto iluRow = ILU.createbegin(); - for(std::size_t i = 0, iend = A.N(); i < iend; ++i) + for (auto&& [_, i] : sparseRange(A)) { - auto& orow = A[inverseOrdering[i]]; - Map rowPattern; - for ( auto col = orow.begin(), cend = orow.end(); col != cend; ++col) + for (auto&& [_, j] : sparseRange(A[inverseOrdering[i]])) { - rowPattern[ordering[col.index()]] = 0; + rowPattern[ordering[j]] = 0; } for(auto ik = rowPattern.begin(); ik->first < i; ++ik) { if ( ik->second < n ) { - auto& rowk = ILU[ik->first]; + auto rowk = ILU.begin() + ik->first; - for ( auto kj = rowk.find(ik->first), endk = rowk.end(); + for ( auto kj = rowk->find(ik->first), endk = rowk->end(); kj != endk; ++kj) { // Assume double and block_type FieldMatrix @@ -220,14 +218,14 @@ void milun_decomposition(const M& A, int n, MILU_VARIANT milu, M& ILU, // copy Entries from A for(auto iter=A.begin(), iend = A.end(); iter != iend; ++iter) { - auto& newRow = ILU[ordering[iter.index()]]; // reset stored generation - std::ranges::fill(newRow, 0); + auto newRowIt = ILU.begin() + ordering[iter.index()]; + std::ranges::fill(*newRowIt, 0); // copy row. for(auto col = iter->begin(), cend = iter->end(); col != cend; ++col) { - newRow[ordering[col.index()]] = *col; + (*newRowIt)[ordering[col.index()]] = *col; } } // call decomposition on pattern diff --git a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp index 4e1f0fa779..87b27898d7 100644 --- a/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp +++ b/opm/simulators/linalg/ParallelOverlappingILU0_impl.hpp @@ -155,7 +155,7 @@ void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv) assert(colcount == numLower); - const auto rendi = A.beforeBegin(); + const auto rbegini = std::make_reverse_iterator(A.begin()); row = 0; colcount = 0; upper.rows_[ 0 ] = colcount ; @@ -164,21 +164,22 @@ void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv) // NOTE: upper and inv store entries in reverse order, reverse here // relative to ILU - for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row ) + auto rindex = [](auto it) { return std::prev(it.base()).index(); }; + for (auto i=std::make_reverse_iterator(A.end()); i!=rbegini; ++i, ++ row ) { - const size_type iIndex = i.index(); + const size_type iIndex = rindex(i); // store in reverse row order // eliminate entries left of diagonal; store L factor - for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j ) + for (auto j=std::make_reverse_iterator(i->end()); rindex(j)>=iIndex; ++j ) { - const size_type jIndex = j.index(); - if( j.index() == iIndex ) + const size_type jIndex = rindex(j); + if( rindex(j) == iIndex ) { inv[ row ] = (*j); break; } - else if ( j.index() >= i.index() ) + else if ( rindex(j) >= rindex(i) ) { upper.push_back( (*j), jIndex ); ++colcount ; @@ -465,9 +466,8 @@ update() // have changed, so we must copy all elements. for (std::size_t row = 0; row < A_->N(); ++row) { const auto& Arow = (*A_)[row]; - auto& ILUrow = (*ILU_)[row]; auto Ait = Arow.begin(); - auto Iit = ILUrow.begin(); + auto Iit = (*ILU_)[row].begin(); for (; Ait != Arow.end(); ++Ait, ++Iit) { *Iit = *Ait; } @@ -484,7 +484,8 @@ update() A_->nonzeroes(), Matrix::row_wise); auto& newA = *ILU_; // Create sparsity pattern - for (auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter) + auto endcreateA = newA.createend(); + for (auto iter = newA.createbegin(); iter != endcreateA; ++iter) { const auto& row = (*A_)[inverseOrdering[iter.index()]]; for (auto col = row.begin(), cend = row.end(); col != cend; ++col) @@ -493,12 +494,12 @@ update() } } // Copy values - for (auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter) + for (auto iter = A_->begin(); iter != A_->end(); ++iter) { - auto& newRow = newA[ordering_[iter.index()]]; - for (auto col = iter->begin(), cend = iter->end(); col != cend; ++col) + auto newRow = newA.begin() + ordering_[iter.index()]; + for (auto&& [A_ij, j] : sparseRange(*newRow)) { - newRow[ordering_[col.index()]] = *col; + (*newRow)[ordering_[j]] = A_ij; } } } diff --git a/opm/simulators/linalg/gpubridge/GpuBridge.cpp b/opm/simulators/linalg/gpubridge/GpuBridge.cpp index 40411cb87a..8d49c22213 100644 --- a/opm/simulators/linalg/gpubridge/GpuBridge.cpp +++ b/opm/simulators/linalg/gpubridge/GpuBridge.cpp @@ -160,7 +160,7 @@ int replaceZeroDiagonal(BridgeMatrix& mat, } else { for (typename BridgeMatrix::iterator r = mat.begin(); r != mat.end(); ++r) { typename BridgeMatrix::size_type offset = diag_indices[r.index()]; - auto& diag_block = r->getptr()[offset]; // diag_block is a reference to MatrixBlock, located on column r of row r + auto& diag_block = *std::next(r->begin(), offset); // diag_block is a reference to MatrixBlock, located on column r of row r for (int rr = 0; rr < dim; ++rr) { auto& val = diag_block[rr][rr]; if (val == 0.0) { // could be replaced by '< 1e-30' or similar diff --git a/opm/simulators/linalg/istlsparsematrixadapter.hh b/opm/simulators/linalg/istlsparsematrixadapter.hh index ac59e12188..cd6fd977d9 100644 --- a/opm/simulators/linalg/istlsparsematrixadapter.hh +++ b/opm/simulators/linalg/istlsparsematrixadapter.hh @@ -138,15 +138,11 @@ public: for (int i = 0; i < diagBlock.rows; ++i) diagBlock[i][i] = diag; - auto& matRow = (*istlMatrix_)[row]; - auto colIt = matRow.begin(); - const auto& colEndIt = matRow.end(); - for (; colIt != colEndIt; ++colIt) { - if (colIt.index() == row) - *colIt = diagBlock; + for (auto&& [a_ij, j] : sparseRange((*istlMatrix_)[row])) + if (j == row) + a_ij = diagBlock; else - *colIt = Scalar(0.0); - } + a_ij = Scalar(0.0); } /*! diff --git a/opm/simulators/linalg/overlappingbcrsmatrix.hh b/opm/simulators/linalg/overlappingbcrsmatrix.hh index d6f32fc336..2da3c0d148 100644 --- a/opm/simulators/linalg/overlappingbcrsmatrix.hh +++ b/opm/simulators/linalg/overlappingbcrsmatrix.hh @@ -328,8 +328,8 @@ private: size_t numDomestic = overlap_->numDomestic(); // allocate the rows - this->setSize(numDomestic, numDomestic); this->setBuildMode(ParentType::random); + this->setSize(numDomestic, numDomestic); // communicate the entries buildIndices_(nativeMatrix); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index b52189bcb5..86da9b0e5f 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1493,15 +1493,12 @@ namespace Opm { { int nw = this->numLocalWellsEnd(); int rdofs = local_num_cells_; - for (int i = 0; i < nw; ++i) { - int wdof = rdofs + i; - jacobian.entry(wdof,wdof) = 1.0;// better scaling ? - } const auto wellconnections = this->getMaxWellConnections(); for (int i = 0; i < nw; ++i) { + int wdof = rdofs + i; + jacobian.entry(wdof,wdof) = 0.0; const auto& perfcells = wellconnections[i]; for (int perfcell : perfcells) { - int wdof = rdofs + i; jacobian.entry(wdof, perfcell) = 0.0; jacobian.entry(perfcell, wdof) = 0.0; } diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index def1efa40f..a57448f15a 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -86,8 +86,8 @@ init(const int numPerfs, duneC_.setSize(well_.numberOfSegments(), numPerfs, numPerfs); // we need to add the off diagonal ones - for (auto row = duneD_.createbegin(), - end = duneD_.createend(); row != end; ++row) { + auto endD = duneD_.createend(); + for (auto row = duneD_.createbegin(); row != endD; ++row) { // the number of the row corrspnds to the segment now const int seg = row.index(); // adding the item related to outlet relation @@ -108,8 +108,8 @@ init(const int numPerfs, } // make the C matrix - for (auto row = duneC_.createbegin(), - end = duneC_.createend(); row != end; ++row) { + auto endC = duneC_.createend(); + for (auto row = duneC_.createbegin(); row != endC; ++row) { // the number of the row corresponds to the segment number now. for (const int& perf : perforations[row.index()]) { const int local_perf_index = parallel_well_info.activePerfToLocalPerf(perf); @@ -120,8 +120,8 @@ init(const int numPerfs, } // make the B^T matrix - for (auto row = duneB_.createbegin(), - end = duneB_.createend(); row != end; ++row) { + auto endB = duneB_.createend(); + for (auto row = duneB_.createbegin(); row != endB; ++row) { // the number of the row corresponds to the segment number now. for (const int& perf : perforations[row.index()]) { const int local_perf_index = parallel_well_info.activePerfToLocalPerf(perf); diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp index 39c0d88b3b..9efa27445f 100644 --- a/opm/simulators/wells/StandardWellEquations.cpp +++ b/opm/simulators/wells/StandardWellEquations.cpp @@ -47,8 +47,6 @@ StandardWellEquations:: StandardWellEquations(const ParallelWellInfo& parallel_well_info) : parallelB_(duneB_, parallel_well_info) { - duneB_.setBuildMode(OffDiagMatWell::row_wise); - duneC_.setBuildMode(OffDiagMatWell::row_wise), invDuneD_.setBuildMode(DiagMatWell::row_wise); } @@ -62,20 +60,24 @@ init(const int numWellEq, //[A C^T [x = [ res // B D] x_well] res_well] // set the size of the matrices + duneB_.setBuildMode(OffDiagMatWell::row_wise); + duneC_.setBuildMode(OffDiagMatWell::row_wise), + duneD_.setBuildMode(DiagMatWell::row_wise), + duneD_.setSize(1, 1, 1); duneB_.setSize(1, numPerfs, numPerfs); duneC_.setSize(1, numPerfs, numPerfs); - for (auto row = duneD_.createbegin(), - end = duneD_.createend(); row != end; ++row) { + auto endD = duneD_.createend(); + for (auto row = duneD_.createbegin(); row != endD; ++row) { // Add nonzeros for diagonal row.insert(row.index()); } // the block size is run-time determined now duneD_[0][0].resize(numWellEq, numWellEq); - for (auto row = duneB_.createbegin(), - end = duneB_.createend(); row != end; ++row) { + auto endB = duneB_.createend(); + for (auto row = duneB_.createbegin(); row != endB; ++row) { for (int perf = 0 ; perf < numPerfs; ++perf) { row.insert(perf); } @@ -87,8 +89,8 @@ init(const int numWellEq, } // make the C^T matrix - for (auto row = duneC_.createbegin(), - end = duneC_.createend(); row != end; ++row) { + auto endC = duneC_.createend(); + for (auto row = duneC_.createbegin(); row != endC; ++row) { for (int perf = 0; perf < numPerfs; ++perf) { row.insert(perf); } diff --git a/opm/simulators/wells/WellHelpers.cpp b/opm/simulators/wells/WellHelpers.cpp index ab2e5673c9..b3b51a4ea5 100644 --- a/opm/simulators/wells/WellHelpers.cpp +++ b/opm/simulators/wells/WellHelpers.cpp @@ -31,6 +31,8 @@ #include +#include + #include #include