From 3680c21f329024336ecfbd42282e7be89eb7d90a Mon Sep 17 00:00:00 2001 From: jakobtorben Date: Mon, 7 Oct 2024 17:44:23 +0200 Subject: [PATCH] Map perforated cell index to global index for CPR weights and add well contrib to matrix --- .../wells/MultisegmentWellEquations.cpp | 20 +++++++++++----- .../wells/MultisegmentWellEquations.hpp | 3 +++ .../wells/StandardWellEquations.cpp | 24 ++++++++++++------- .../wells/StandardWellEquations.hpp | 3 +++ 4 files changed, 36 insertions(+), 14 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index a117e7f20..411b06a9c 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -122,6 +122,9 @@ init(const int num_cells, } resWell_.resize(well_.numberOfSegments()); + + // Store the global index of well perforated cells + cells_ = cells; } template @@ -297,11 +300,12 @@ extract(SparseMatrixAdapter& jacobian) const for (std::size_t rowC = 0; rowC < duneC_.N(); ++rowC) { for (auto colC = duneC_[rowC].begin(), endC = duneC_[rowC].end(); colC != endC; ++colC) { - const auto row_index = colC.index(); + // map the well perforated cell index to global cell index + const auto row_index = cells_[colC.index()]; for (std::size_t rowB = 0; rowB < duneB_.N(); ++rowB) { for (auto colB = duneB_[rowB].begin(), endB = duneB_[rowB].end(); colB != endB; ++colB) { - const auto col_index = colB.index(); + const auto col_index = cells_[colB.index()]; OffDiagMatrixBlockWellType tmp1; detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); typename SparseMatrixAdapter::MatrixBlock tmp2; @@ -327,12 +331,14 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, // Add the pressure contribution to the cpr system for the well // Add for coupling from well to reservoir - const int welldof_ind = duneC_.M() + well.indexOfWell(); + const int number_cells = weights.size(); + const int welldof_ind = number_cells + well.indexOfWell(); if (!well.isPressureControlled(well_state)) { for (std::size_t rowC = 0; rowC < duneC_.N(); ++rowC) { for (auto colC = duneC_[rowC].begin(), endC = duneC_[rowC].end(); colC != endC; ++colC) { - const auto row_index = colC.index(); + // map the well perforated cell index to global cell index + const auto row_index = cells_[colC.index()]; const auto& bw = weights[row_index]; double matel = 0.0; @@ -352,7 +358,8 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, for (std::size_t rowB = 0; rowB < duneB_.N(); ++rowB) { for (auto colB = duneB_[rowB].begin(), endB = duneB_[rowB].end(); colB != endB; ++colB) { - const auto col_index = colB.index(); + // map the well perforated cell index to global cell index + const auto col_index = cells_[colB.index()]; const auto& bw = weights[col_index]; well_weight += bw; num_perfs += 1; @@ -368,7 +375,8 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, const auto& bw = well_weight; for (auto colB = duneB_[rowB].begin(), endB = duneB_[rowB].end(); colB != endB; ++colB) { - const auto col_index = colB.index(); + // map the well perforated cell index to global cell index + const auto col_index = cells_[colB.index()]; double matel = 0.0; for (std::size_t i = 0; i< bw.size(); ++i) { matel += bw[i] *(*colB)[i][pressureVarIndex]; diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index f204c0699..ae2944f9c 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -145,6 +145,9 @@ public: BVectorWell resWell_; const MultisegmentWellGeneric& well_; //!< Reference to well + + // Store the global index of well perforated cells + std::vector cells_; }; } diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp index b45440ebd..b2e8d0c74 100644 --- a/opm/simulators/wells/StandardWellEquations.cpp +++ b/opm/simulators/wells/StandardWellEquations.cpp @@ -111,6 +111,9 @@ init(const int num_cells, for (unsigned i = 0; i < duneD_.N(); ++i) { invDrw_[i].resize(numWellEq); } + + // Store the global index of well perforated cells + cells_ = cells; } template @@ -261,14 +264,17 @@ extract(SparseMatrixAdapter& jacobian) const for (auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) { - const auto row_index = colC.index(); + // map the well perforated cell index to global cell index + const auto row_index = this->cells_[colC.index()]; for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) { + // map the well perforated cell index to global cell index + const auto col_index = this->cells_[colB.index()]; detail::multMatrix(invDuneD_[0][0], (*colB), tmp); detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat); - jacobian.addToBlock(row_index, colB.index(), tmpMat); + jacobian.addToBlock(row_index, col_index, tmpMat); } } } @@ -307,22 +313,23 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, int nperf = 0; auto cell_weights = weights[0];// not need for not(use_well_weights) cell_weights = 0.0; - assert(duneC_.M() == weights.size()); - const int welldof_ind = duneC_.M() + well.indexOfWell(); + const int number_cells = weights.size(); + const int welldof_ind = number_cells + well.indexOfWell(); // do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also) if (!well.isPressureControlled(well_state) || use_well_weights) { // make coupling for reservoir to well for (auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC) { - const auto row_ind = colC.index(); - const auto& bw = weights[row_ind]; + // map the well perforated cell index to global cell index + const auto row_index = cells_[colC.index()]; + const auto& bw = weights[row_index]; Scalar matel = 0; assert((*colC).M() == bw.size()); for (std::size_t i = 0; i < bw.size(); ++i) { matel += (*colC)[bhp_var_index][i] * bw[i]; } - jacobian[row_ind][welldof_ind] = matel; + jacobian[row_index][welldof_ind] = matel; cell_weights += bw; nperf += 1; } @@ -381,7 +388,8 @@ extractCPRPressureMatrix(PressureMatrix& jacobian, if (!well.isPressureControlled(well_state) || use_well_weights) { for (auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB) { - const auto col_index = colB.index(); + // map the well perforated cell index to global cell index + const auto col_index = cells_[colB.index()]; const auto& bw = bweights[0]; Scalar matel = 0; for (std::size_t i = 0; i < bw.size(); ++i) { diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp index 736d7ded1..76d21ef05 100644 --- a/opm/simulators/wells/StandardWellEquations.hpp +++ b/opm/simulators/wells/StandardWellEquations.hpp @@ -150,6 +150,9 @@ private: // several vector used in the matrix calculation mutable BVectorWell Bx_; mutable BVectorWell invDrw_; + + // Store the global index of well perforated cells + std::vector cells_; }; }