My Project
PressureBhpTransferPolicy.hpp
1 /*
2  Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2019 Dr. Blatt - HPC-Simulation-Software & Services.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #pragma once
22 
23 #include <opm/simulators/linalg/matrixblock.hh>
25 #include <dune/istl/paamg/pinfo.hh>
26 
27 namespace Opm
28 {
29  template <class Communication>
30  void extendCommunicatorWithWells(const Communication& comm,
31  std::shared_ptr<Communication>& commRW,
32  const int nw)
33  {
34  // used for extending the coarse communicator pattern
35  using IndexSet = typename Communication::ParallelIndexSet;
36  using LocalIndex = typename IndexSet::LocalIndex;
37  const IndexSet& indset = comm.indexSet();
38  IndexSet& indset_rw = commRW->indexSet();
39  const int max_nw = comm.communicator().max(nw);
40  const int rank = comm.communicator().rank();
41  int glo_max = 0;
42  size_t loc_max = 0;
43  indset_rw.beginResize();
44  for (auto ind = indset.begin(), indend = indset.end(); ind != indend; ++ind) {
45  indset_rw.add(ind->global(), LocalIndex(ind->local(), ind->local().attribute(), true));
46  const int glo = ind->global();
47  const size_t loc = ind->local().local();
48  glo_max = std::max(glo_max, glo);
49  loc_max = std::max(loc_max, loc);
50  }
51  const int global_max = comm.communicator().max(glo_max);
52  // used append the welldofs at the end
53  assert(loc_max + 1 == indset.size());
54  size_t local_ind = loc_max + 1;
55  for (int i = 0; i < nw; ++i) {
56  // need to set unique global number
57  const size_t v = global_max + max_nw * rank + i + 1;
58  // set to true to not have problems with higher levels if growing of domains is used
59  indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner, true));
60  ++local_ind;
61  }
62  indset_rw.endResize();
63  assert(indset_rw.size() == indset.size() + nw);
64  // assume same communication pattern
65  commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
66  commRW->remoteIndices().template rebuild<true>();
67  //commRW->clearInterfaces(); may need this for correct rebuild
68  }
69 
70  namespace Details
71  {
72  using PressureMatrixType = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
73  using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
74  using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
75  template <class Comm>
76  using ParCoarseOperatorType
77  = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Comm>;
78  template <class Comm>
79  using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
80  SeqCoarseOperatorType,
81  ParCoarseOperatorType<Comm>>;
82  } // namespace Details
83 
84  template <class FineOperator, class Communication, bool transpose = false>
85  class PressureBhpTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
86  {
87  public:
88  typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
90  typedef Communication ParallelInformation;
91  typedef typename FineOperator::domain_type FineVectorType;
92 
93  public:
94  PressureBhpTransferPolicy(const Communication& comm,
95  const FineVectorType& weights,
96  const Opm::PropertyTree& prm,
97  const std::size_t pressureIndex)
98  : communication_(&const_cast<Communication&>(comm))
99  , weights_(weights)
100  , prm_(prm)
101  , pressure_var_index_(pressureIndex)
102  {
103  }
104 
105  virtual void createCoarseLevelSystem(const FineOperator& fineOperator) override
106  {
107  using CoarseMatrix = typename CoarseOperator::matrix_type;
108  const auto& fineLevelMatrix = fineOperator.getmat();
109  const auto& nw = fineOperator.getNumberOfExtraEquations();
110  if (prm_.get<bool>("add_wells")) {
111  const size_t average_elements_per_row
112  = static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
113  const double overflow_fraction = 1.2;
114  coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N() + nw,
115  fineLevelMatrix.M() + nw,
116  average_elements_per_row,
117  overflow_fraction,
118  CoarseMatrix::implicit));
119  int rownum = 0;
120  for (const auto& row : fineLevelMatrix) {
121  for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
122  coarseLevelMatrix_->entry(rownum, col.index()) = 0.0;
123  }
124  ++rownum;
125  }
126  } else {
127  coarseLevelMatrix_.reset(
128  new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
129  auto createIter = coarseLevelMatrix_->createbegin();
130  for (const auto& row : fineLevelMatrix) {
131  for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
132  createIter.insert(col.index());
133  }
134  ++createIter;
135  }
136  }
137 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
138  if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
139  coarseLevelCommunication_ = std::make_shared<Communication>();
140  } else {
141  coarseLevelCommunication_ = std::make_shared<Communication>(
142  communication_->communicator(), communication_->category(), false);
143  }
144 #else
145  if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
146  coarseLevelCommunication_ = std::make_shared<Communication>();
147  } else {
148  coarseLevelCommunication_ = std::make_shared<Communication>(
149  communication_->communicator(), communication_->getSolverCategory(), false);
150  }
151 #endif
152  if (prm_.get<bool>("add_wells")) {
153  fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
154  coarseLevelMatrix_->compress(); // all elemenst should be set
155  if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
156  extendCommunicatorWithWells(*communication_, coarseLevelCommunication_, nw);
157  }
158  }
159  calculateCoarseEntries(fineOperator);
160 
161  this->lhs_.resize(this->coarseLevelMatrix_->M());
162  this->rhs_.resize(this->coarseLevelMatrix_->N());
163  using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
164 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
165  OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
166  this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
167 #else
168  OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
169  this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
170 #endif
171  }
172 
173  virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
174  {
175  const auto& fineMatrix = fineOperator.getmat();
176  *coarseLevelMatrix_ = 0;
177  auto rowCoarse = coarseLevelMatrix_->begin();
178  for (auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
179  assert(row.index() == rowCoarse.index());
180  auto entryCoarse = rowCoarse->begin();
181  for (auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
182  assert(entry.index() == entryCoarse.index());
183  double matrix_el = 0;
184  if (transpose) {
185  const auto& bw = weights_[entry.index()];
186  for (size_t i = 0; i < bw.size(); ++i) {
187  matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
188  }
189  } else {
190  const auto& bw = weights_[row.index()];
191  for (size_t i = 0; i < bw.size(); ++i) {
192  matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
193  }
194  }
195  (*entryCoarse) = matrix_el;
196  }
197  }
198  if (prm_.get<bool>("add_wells")) {
199  assert(transpose == false); // not implemented
200  bool use_well_weights = prm_.get<bool>("use_well_weights");
201  fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
202 #ifndef NDEBUG
203  std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
204  assert(rowCoarse == coarseLevelMatrix_->end());
205 #endif
206 
207  }
208  }
209 
210  virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
211  {
212  //NB we iterate over fine assumming welldofs is at the end
213  // Set coarse vector to zero
214  this->rhs_ = 0;
215 
216  auto end = fine.end(), begin = fine.begin();
217 
218  for (auto block = begin; block != end; ++block) {
219  const auto& bw = weights_[block.index()];
220  double rhs_el = 0.0;
221  if (transpose) {
222  rhs_el = (*block)[pressure_var_index_];
223  } else {
224  for (size_t i = 0; i < block->size(); ++i) {
225  rhs_el += (*block)[i] * bw[i];
226  }
227  }
228  this->rhs_[block - begin] = rhs_el;
229  }
230 
231  this->lhs_ = 0;
232  }
233 
234  virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
235  {
236  //NB we iterate over fine assumming welldofs is at the end
237  auto end = fine.end(), begin = fine.begin();
238 
239  for (auto block = begin; block != end; ++block) {
240  if (transpose) {
241  const auto& bw = weights_[block.index()];
242  for (size_t i = 0; i < block->size(); ++i) {
243  (*block)[i] = this->lhs_[block - begin] * bw[i];
244  }
245  } else {
246  (*block)[pressure_var_index_] = this->lhs_[block - begin];
247  }
248  }
249  }
250 
251  virtual PressureBhpTransferPolicy* clone() const override
252  {
253  return new PressureBhpTransferPolicy(*this);
254  }
255 
256  const Communication& getCoarseLevelCommunication() const
257  {
258  return *coarseLevelCommunication_;
259  }
260 
261 private:
262  Communication* communication_;
263  const FineVectorType& weights_;
264  PropertyTree prm_;
265  const int pressure_var_index_;
266  std::shared_ptr<Communication> coarseLevelCommunication_;
267  std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
268 };
269 
270 } // namespace Opm
271 
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
Definition: PressureBhpTransferPolicy.hpp:86
virtual PressureBhpTransferPolicy * clone() const override
Clone the current object.
Definition: PressureBhpTransferPolicy.hpp:251
virtual void moveToFineLevel(typename ParentType::FineDomainType &fine) override
Updates the fine level linear system after the correction of the coarse levels system.
Definition: PressureBhpTransferPolicy.hpp:234
virtual void calculateCoarseEntries(const FineOperator &fineOperator) override
???.
Definition: PressureBhpTransferPolicy.hpp:173
virtual void createCoarseLevelSystem(const FineOperator &fineOperator) override
Algebraically creates the coarse level system.
Definition: PressureBhpTransferPolicy.hpp:105
Definition: PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Algebraic twolevel methods.