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

Make well equations use only perforated cells

Browse Source
This commit is contained in:
jakobtorben 2024-05-08 12:54:57 +02:00
parent 8ae31797b1
commit 60ece76adf
5 changed files with 76 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -229,6 +229,8 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
{
GlobalEqVector res_local;
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& well: well_container_) {
// Modifiy the Jacobian with explicit Schur complement
@ -238,7 +240,19 @@ namespace Opm {
}
// Apply as Schur complement the well residual to reservoir residuals:
// r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res);
// Well equations B and C uses only the perforated cells, so need to apply on local residual
auto cells = well->cells();
res_local.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
res_local[i] = res[cells[i]];
}
well->apply(res_local);
for (size_t i = 0; i < cells.size(); ++i) {
res[cells[i]] = res_local[i];
}
}
OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ",
simulator_.gridView().comm());
@ -250,6 +264,8 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res)
{
GlobalEqVector res_local;
// Note: no point in trying to do a parallel gathering
// try/catch here, as this function is not called in
// parallel but for each individual domain of each rank.
@ -262,7 +278,19 @@ namespace Opm {
}
// Apply as Schur complement the well residual to reservoir residuals:
// r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res);
// Well equations B and C uses only the perforated cells, so need to apply on local residual
auto cells = well->cells();
res_local.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
res_local[i] = res[cells[i]];
}
well->apply(res_local);
for (size_t i = 0; i < cells.size(); ++i) {
res[cells[i]] = res_local[i];
}
}
}
}
@ -1712,8 +1740,20 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
apply(BVector& r) const
{
BVector r_local;
for (auto& well : well_container_) {
well->apply(r);
auto cells = well->cells();
r_local.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
r_local[i] = r[cells[i]];
}
well->apply(r_local);
for (size_t i = 0; i < cells.size(); ++i) {
r[cells[i]] = r_local[i];
}
}
}
@ -1724,8 +1764,26 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
BVector x_local;
BVector Ax_local;
for (auto& well : well_container_) {
well->apply(x, Ax);
// Well equations B and C uses only the perforated cells, so need to apply on local vectors
auto cells = well->cells();
x_local.resize(cells.size());
Ax_local.resize(cells.size());
for (size_t i = 0; i < cells.size(); ++i) {
x_local[i] = x[cells[i]];
Ax_local[i] = Ax[cells[i]];
}
well->apply(x_local, Ax_local);
for (size_t i = 0; i < cells.size(); ++i) {
// only need to update Ax
Ax[cells[i]] = Ax_local[i];
}
}
}

10
opm/simulators/wells/MultisegmentWellEquations.cpp
View File

@ -78,8 +78,8 @@ init(const int num_cells,
}
duneD_.setSize(well_.numberOfSegments(), well_.numberOfSegments(), nnz_d);
}
duneB_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
duneC_.setSize(well_.numberOfSegments(), num_cells, numPerfs);
duneB_.setSize(well_.numberOfSegments(), numPerfs, numPerfs);
duneC_.setSize(well_.numberOfSegments(), numPerfs, numPerfs);
// we need to add the off diagonal ones
for (auto row = duneD_.createbegin(),
@ -108,8 +108,7 @@ init(const int num_cells,
end = duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
row.insert(perf);
}
}
@ -118,8 +117,7 @@ init(const int num_cells,
end = duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
row.insert(perf);
}
}

2
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -1904,7 +1904,7 @@ namespace Opm
this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
MultisegmentWellAssemble(*this).
assemblePerforationEq(seg, cell_idx, comp_idx, cq_s_effective, this->linSys_);
assemblePerforationEq(seg, perf, comp_idx, cq_s_effective, this->linSys_);
}
}

16
opm/simulators/wells/StandardWellEquations.cpp
View File

@ -62,8 +62,8 @@ init(const int num_cells,
// B D] x_well] res_well]
// set the size of the matrices
duneD_.setSize(1, 1, 1);
duneB_.setSize(1, num_cells, numPerfs);
duneC_.setSize(1, num_cells, numPerfs);
duneB_.setSize(1, numPerfs, numPerfs);
duneC_.setSize(1, numPerfs, numPerfs);
for (auto row = duneD_.createbegin(),
end = duneD_.createend(); row != end; ++row) {
@ -76,29 +76,25 @@ init(const int num_cells,
for (auto row = duneB_.createbegin(),
end = duneB_.createend(); row != end; ++row) {
for (int perf = 0 ; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
row.insert(perf);
}
}
for (int perf = 0 ; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
// the block size is run-time determined now
duneB_[0][cell_idx].resize(numWellEq, numEq);
duneB_[0][perf].resize(numWellEq, numEq);
}
// make the C^T matrix
for (auto row = duneC_.createbegin(),
end = duneC_.createend(); row != end; ++row) {
for (int perf = 0; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
row.insert(perf);
}
}
for (int perf = 0; perf < numPerfs; ++perf) {
const int cell_idx = cells[perf];
duneC_[0][cell_idx].resize(numWellEq, numEq);
duneC_[0][perf].resize(numWellEq, numEq);
}
resWell_.resize(1);

7
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -404,7 +404,6 @@ namespace Opm
water_flux_s, deferred_logger);
}
}
const int cell_idx = this->well_cells_[perf];
for (int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
// the cq_s entering mass balance equations need to consider the efficiency factors.
const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
@ -414,7 +413,7 @@ namespace Opm
StandardWellAssemble<FluidSystem,Indices>(*this).
assemblePerforationEq(cq_s_effective,
componentIdx,
cell_idx,
perf,
this->primary_variables_.numWellEq(),
this->linSys_);
@ -430,7 +429,7 @@ namespace Opm
if constexpr (has_zFraction) {
StandardWellAssemble<FluidSystem,Indices>(*this).
assembleZFracEq(cq_s_zfrac_effective,
cell_idx,
perf,
this->primary_variables_.numWellEq(),
this->linSys_);
}
@ -2132,7 +2131,7 @@ namespace Opm
eq_wat_vel,
pskin_index,
wat_vel_index,
cell_idx,
perf,
this->primary_variables_.numWellEq(),
this->linSys_);
}