My Project
PressureTransferPolicy.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 #ifndef OPM_PRESSURE_TRANSFER_POLICY_HEADER_INCLUDED
22 #define OPM_PRESSURE_TRANSFER_POLICY_HEADER_INCLUDED
23 
24 
26 #include <opm/simulators/linalg/PropertyTree.hpp>
27 #include <opm/simulators/linalg/matrixblock.hh>
28 
29 
30 namespace Opm
31 {
32 
33 namespace Details
34 {
35  using PressureMatrixType = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
36  using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
37  using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
38  template <class Comm>
39  using ParCoarseOperatorType
40  = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Comm>;
41  template <class Comm>
42  using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
43  SeqCoarseOperatorType,
44  ParCoarseOperatorType<Comm>>;
45 } // namespace Details
46 
47 
48 
49 template <class FineOperator, class Communication, bool transpose = false>
51  : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
52 {
53 public:
54  typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
56  typedef Communication ParallelInformation;
57  typedef typename FineOperator::domain_type FineVectorType;
58 
59 public:
60  PressureTransferPolicy(const Communication& comm,
61  const FineVectorType& weights,
62  const Opm::PropertyTree& /*prm*/,
63  int pressure_var_index)
64  : communication_(&const_cast<Communication&>(comm))
65  , weights_(weights)
66  , pressure_var_index_(pressure_var_index)
67  {
68  }
69 
70  virtual void createCoarseLevelSystem(const FineOperator& fineOperator) override
71  {
72  using CoarseMatrix = typename CoarseOperator::matrix_type;
73  const auto& fineLevelMatrix = fineOperator.getmat();
74  coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
75  auto createIter = coarseLevelMatrix_->createbegin();
76 
77  for (const auto& row : fineLevelMatrix) {
78  for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
79  createIter.insert(col.index());
80  }
81  ++createIter;
82  }
83 
84  calculateCoarseEntries(fineOperator);
85  coarseLevelCommunication_.reset(communication_, [](Communication*) {});
86 
87  this->lhs_.resize(this->coarseLevelMatrix_->M());
88  this->rhs_.resize(this->coarseLevelMatrix_->N());
89  using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
90 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
91  OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
92  this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
93 #else
94  OperatorArgs oargs(*coarseLevelMatrix_, *coarseLevelCommunication_);
95  this->operator_.reset(Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs));
96 #endif
97  }
98 
99  virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
100  {
101  const auto& fineMatrix = fineOperator.getmat();
102  *coarseLevelMatrix_ = 0;
103  auto rowCoarse = coarseLevelMatrix_->begin();
104  for (auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
105  assert(row.index() == rowCoarse.index());
106  auto entryCoarse = rowCoarse->begin();
107  for (auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
108  assert(entry.index() == entryCoarse.index());
109  double matrix_el = 0;
110  if (transpose) {
111  const auto& bw = weights_[entry.index()];
112  for (size_t i = 0; i < bw.size(); ++i) {
113  matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
114  }
115  } else {
116  const auto& bw = weights_[row.index()];
117  for (size_t i = 0; i < bw.size(); ++i) {
118  matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
119  }
120  }
121  (*entryCoarse) = matrix_el;
122  }
123  }
124  assert(rowCoarse == coarseLevelMatrix_->end());
125  }
126 
127  virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
128  {
129  // Set coarse vector to zero
130  this->rhs_ = 0;
131 
132  auto end = fine.end(), begin = fine.begin();
133 
134  for (auto block = begin; block != end; ++block) {
135  const auto& bw = weights_[block.index()];
136  double rhs_el = 0.0;
137  if (transpose) {
138  rhs_el = (*block)[pressure_var_index_];
139  } else {
140  for (size_t i = 0; i < block->size(); ++i) {
141  rhs_el += (*block)[i] * bw[i];
142  }
143  }
144  this->rhs_[block - begin] = rhs_el;
145  }
146 
147  this->lhs_ = 0;
148  }
149 
150  virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
151  {
152  auto end = fine.end(), begin = fine.begin();
153 
154  for (auto block = begin; block != end; ++block) {
155  if (transpose) {
156  const auto& bw = weights_[block.index()];
157  for (size_t i = 0; i < block->size(); ++i) {
158  (*block)[i] = this->lhs_[block - begin] * bw[i];
159  }
160  } else {
161  (*block)[pressure_var_index_] = this->lhs_[block - begin];
162  }
163  }
164  }
165 
166  virtual PressureTransferPolicy* clone() const override
167  {
168  return new PressureTransferPolicy(*this);
169  }
170 
171  const Communication& getCoarseLevelCommunication() const
172  {
173  return *coarseLevelCommunication_;
174  }
175 
176  std::size_t getPressureIndex() const
177  {
178  return pressure_var_index_;
179  }
180 private:
181  Communication* communication_;
182  const FineVectorType& weights_;
183  const std::size_t pressure_var_index_;
184  std::shared_ptr<Communication> coarseLevelCommunication_;
185  std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
186 };
187 
188 } // namespace Opm
189 
190 #endif // OPM_PRESSURE_TRANSFER_POLICY_HEADER_INCLUDED
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: PressureTransferPolicy.hpp:52
virtual void createCoarseLevelSystem(const FineOperator &fineOperator) override
Algebraically creates the coarse level system.
Definition: PressureTransferPolicy.hpp:70
virtual void calculateCoarseEntries(const FineOperator &fineOperator) override
???.
Definition: PressureTransferPolicy.hpp:99
virtual PressureTransferPolicy * clone() const override
Clone the current object.
Definition: PressureTransferPolicy.hpp:166
virtual void moveToFineLevel(typename ParentType::FineDomainType &fine) override
Updates the fine level linear system after the correction of the coarse levels system.
Definition: PressureTransferPolicy.hpp:150
Definition: PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Algebraic twolevel methods.