23#include <opm/common/TimingMacros.hpp>
29#include <dune/istl/paamg/pinfo.hh>
36 template <
class Communication>
38 std::shared_ptr<Communication>& commRW,
43 using IndexSet =
typename Communication::ParallelIndexSet;
44 using LocalIndex =
typename IndexSet::LocalIndex;
45 const IndexSet& indset = comm.indexSet();
46 IndexSet& indset_rw = commRW->indexSet();
47 const int max_nw = comm.communicator().max(nw);
48 const int rank = comm.communicator().rank();
50 std::size_t loc_max = 0;
51 indset_rw.beginResize();
52 for (
auto ind = indset.begin(), indend = indset.end(); ind != indend; ++ind) {
53 indset_rw.add(ind->global(), LocalIndex(ind->local(), ind->local().attribute(),
true));
54 const int glo = ind->global();
55 const std::size_t loc = ind->local().local();
56 glo_max = std::max(glo_max, glo);
57 loc_max = std::max(loc_max, loc);
59 const int global_max = comm.communicator().max(glo_max);
61 assert(loc_max + 1 == indset.size());
62 std::size_t local_ind = loc_max + 1;
63 for (
int i = 0; i < nw; ++i) {
65 const std::size_t v = global_max + max_nw * rank + i + 1;
67 indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner,
true));
70 indset_rw.endResize();
71 assert(indset_rw.size() == indset.size() + nw);
73 commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
74 commRW->remoteIndices().template rebuild<true>();
81 template<
class Scalar>
using PressureVectorType = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
85 template<
class Scalar,
class Comm>
91 template<
class Scalar,
class Comm>
92 using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
97 template<
class FineOperator,
class Communication,
class Scalar,
bool transpose = false>
110 const std::size_t pressureIndex)
114 , pressure_var_index_(pressureIndex)
121 using CoarseMatrix =
typename CoarseOperator::matrix_type;
122 const auto& fineLevelMatrix = fineOperator.getmat();
123 const auto& nw = fineOperator.getNumberOfExtraEquations();
124 if (prm_.
get<
bool>(
"add_wells")) {
125 const std::size_t average_elements_per_row
126 =
static_cast<std::size_t
>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
127 const double overflow_fraction = 1.2;
128 coarseLevelMatrix_.reset(
new CoarseMatrix(fineLevelMatrix.N() + nw,
129 fineLevelMatrix.M() + nw,
130 average_elements_per_row,
132 CoarseMatrix::implicit));
134 for (
const auto& row : fineLevelMatrix) {
135 for (
auto col = row.begin(), cend = row.end(); col != cend; ++col) {
136 coarseLevelMatrix_->entry(rownum, col.index()) = 0.0;
141 coarseLevelMatrix_.reset(
142 new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
143 auto createIter = coarseLevelMatrix_->createbegin();
144 for (
const auto& row : fineLevelMatrix) {
145 for (
auto col = row.begin(), cend = row.end(); col != cend; ++col) {
146 createIter.insert(col.index());
151 if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
152 coarseLevelCommunication_ = std::make_shared<Communication>();
154 coarseLevelCommunication_ = std::make_shared<Communication>(
155 communication_->communicator(), communication_->category(),
false);
157 if (prm_.
get<
bool>(
"add_wells")) {
158 fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
159 coarseLevelMatrix_->compress();
160 if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
166 this->
lhs_.resize(this->coarseLevelMatrix_->M());
167 this->
rhs_.resize(this->coarseLevelMatrix_->N());
168 using OperatorArgs =
typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
169 OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
170 this->
operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
176 const auto& fineMatrix = fineOperator.getmat();
177 *coarseLevelMatrix_ = 0;
178 auto rowCoarse = coarseLevelMatrix_->begin();
179 for (
auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
180 assert(row.index() == rowCoarse.index());
181 auto entryCoarse = rowCoarse->begin();
182 for (
auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
183 assert(entry.index() == entryCoarse.index());
184 Scalar matrix_el = 0;
186 const auto& bw = weights_[entry.index()];
187 for (std::size_t i = 0; i < bw.size(); ++i) {
188 matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
191 const auto& bw = weights_[row.index()];
192 for (std::size_t i = 0; i < bw.size(); ++i) {
193 matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
196 (*entryCoarse) = matrix_el;
199 if (prm_.
get<
bool>(
"add_wells")) {
200 OPM_TIMEBLOCK(cprwAddWellEquation);
201 assert(transpose ==
false);
202 bool use_well_weights = prm_.
get<
bool>(
"use_well_weights");
203 fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
205 std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
206 assert(rowCoarse == coarseLevelMatrix_->end());
219 auto end = fine.end(), begin = fine.begin();
221 for (
auto block = begin; block != end; ++block) {
222 const auto& bw = weights_[block.index()];
225 rhs_el = (*block)[pressure_var_index_];
227 for (std::size_t i = 0; i < block->size(); ++i) {
228 rhs_el += (*block)[i] * bw[i];
231 this->
rhs_[block - begin] = rhs_el;
241 auto end = fine.end(), begin = fine.begin();
243 for (
auto block = begin; block != end; ++block) {
245 const auto& bw = weights_[block.index()];
246 for (std::size_t i = 0; i < block->size(); ++i) {
247 (*block)[i] = this->
lhs_[block - begin] * bw[i];
250 (*block)[pressure_var_index_] = this->
lhs_[block - begin];
262 return *coarseLevelCommunication_;
269 const int pressure_var_index_;
270 std::shared_ptr<Communication> coarseLevelCommunication_;
271 std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:285
Abstract base class for transfer between levels and creation of the coarse level system.
Definition: twolevelmethodcpr.hh:44
FineOperatorType::range_type FineRangeType
The type of the range of the fine level operator.
Definition: twolevelmethodcpr.hh:54
CoarseDomainType lhs_
The coarse level lhs.
Definition: twolevelmethodcpr.hh:141
CoarseRangeType rhs_
The coarse level rhs.
Definition: twolevelmethodcpr.hh:139
FineOperatorType::domain_type FineDomainType
The type of the domain of the fine level operator.
Definition: twolevelmethodcpr.hh:58
std::shared_ptr< CoarseOperatorType > operator_
the coarse level linear operator.
Definition: twolevelmethodcpr.hh:143
Dune linear operator that assumes ghost rows are ordered after interior rows. Avoids some computation...
Definition: WellOperators.hpp:340
Definition: PressureBhpTransferPolicy.hpp:99
typename Details::CoarseOperatorType< Scalar, Communication > CoarseOperator
Definition: PressureBhpTransferPolicy.hpp:101
typename FineOperator::domain_type FineVectorType
Definition: PressureBhpTransferPolicy.hpp:104
Communication ParallelInformation
Definition: PressureBhpTransferPolicy.hpp:103
PressureBhpTransferPolicy * clone() const override
Clone the current object.
Definition: PressureBhpTransferPolicy.hpp:255
void moveToFineLevel(typename ParentType::FineDomainType &fine) override
Updates the fine level linear system after the correction of the coarse levels system.
Definition: PressureBhpTransferPolicy.hpp:237
void createCoarseLevelSystem(const FineOperator &fineOperator) override
Algebraically creates the coarse level system.
Definition: PressureBhpTransferPolicy.hpp:118
PressureBhpTransferPolicy(const Communication &comm, const FineVectorType &weights, const PropertyTree &prm, const std::size_t pressureIndex)
Definition: PressureBhpTransferPolicy.hpp:107
const Communication & getCoarseLevelCommunication() const
Definition: PressureBhpTransferPolicy.hpp:260
void calculateCoarseEntries(const FineOperator &fineOperator) override
???.
Definition: PressureBhpTransferPolicy.hpp:173
void moveToCoarseLevel(const typename ParentType::FineRangeType &fine) override
Definition: PressureBhpTransferPolicy.hpp:212
Definition: PropertyTree.hpp:37
T get(const std::string &key) const
Dune::BCRSMatrix< MatrixBlock< Scalar, 1, 1 > > PressureMatrixType
Definition: PressureBhpTransferPolicy.hpp:80
std::conditional_t< std::is_same< Comm, Dune::Amg::SequentialInformation >::value, SeqCoarseOperatorType< Scalar >, ParCoarseOperatorType< Scalar, Comm > > CoarseOperatorType
Definition: PressureBhpTransferPolicy.hpp:94
Dune::BlockVector< Dune::FieldVector< Scalar, 1 > > PressureVectorType
Definition: PressureBhpTransferPolicy.hpp:81
Dune::MatrixAdapter< PressureMatrixType< Scalar >, PressureVectorType< Scalar >, PressureVectorType< Scalar > > SeqCoarseOperatorType
Definition: PressureBhpTransferPolicy.hpp:84
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:37
void extendCommunicatorWithWells(const Communication &comm, std::shared_ptr< Communication > &commRW, const int nw)
Definition: PressureBhpTransferPolicy.hpp:37
Algebraic twolevel methods.