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