Map perforated cell index to global index for CPR weights and add well contrib to matrix

This commit is contained in:
jakobtorben
2024-10-07 17:44:23 +02:00
parent 58a8d5ee41
commit 3680c21f32
4 changed files with 36 additions and 14 deletions

View File

@@ -122,6 +122,9 @@ init(const int num_cells,
}
resWell_.resize(well_.numberOfSegments());
// Store the global index of well perforated cells
cells_ = cells;
}
template<class Scalar, int numWellEq, int numEq>
@@ -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];

View File

@@ -145,6 +145,9 @@ public:
BVectorWell resWell_;
const MultisegmentWellGeneric<Scalar>& well_; //!< Reference to well
// Store the global index of well perforated cells
std::vector<int> cells_;
};
}

View File

@@ -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<class Scalar, int numEq>
@@ -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) {

View File

@@ -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<int> cells_;
};
}