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
27namespace 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
261private:
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 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 PressureBhpTransferPolicy * clone() const override
Clone the current object.
Definition: PressureBhpTransferPolicy.hpp:251
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.