My Project
Loading...
Searching...
No Matches
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/common/TimingMacros.hpp>
24
25#include <opm/simulators/linalg/matrixblock.hh>
26#include <opm/simulators/linalg/PropertyTree.hpp>
28
29#include <dune/istl/paamg/pinfo.hh>
30
31namespace Opm
32{
33 template <class Communication>
34 void extendCommunicatorWithWells(const Communication& comm,
35 std::shared_ptr<Communication>& commRW,
36 const int nw)
37 {
38 OPM_TIMEBLOCK(extendCommunicatorWithWells);
39 // used for extending the coarse communicator pattern
40 using IndexSet = typename Communication::ParallelIndexSet;
41 using LocalIndex = typename IndexSet::LocalIndex;
42 const IndexSet& indset = comm.indexSet();
43 IndexSet& indset_rw = commRW->indexSet();
44 const int max_nw = comm.communicator().max(nw);
45 const int rank = comm.communicator().rank();
46 int glo_max = 0;
47 size_t loc_max = 0;
48 indset_rw.beginResize();
49 for (auto ind = indset.begin(), indend = indset.end(); ind != indend; ++ind) {
50 indset_rw.add(ind->global(), LocalIndex(ind->local(), ind->local().attribute(), true));
51 const int glo = ind->global();
52 const size_t loc = ind->local().local();
53 glo_max = std::max(glo_max, glo);
54 loc_max = std::max(loc_max, loc);
55 }
56 const int global_max = comm.communicator().max(glo_max);
57 // used append the welldofs at the end
58 assert(loc_max + 1 == indset.size());
59 size_t local_ind = loc_max + 1;
60 for (int i = 0; i < nw; ++i) {
61 // need to set unique global number
62 const size_t v = global_max + max_nw * rank + i + 1;
63 // set to true to not have problems with higher levels if growing of domains is used
64 indset_rw.add(v, LocalIndex(local_ind, Dune::OwnerOverlapCopyAttributeSet::owner, true));
65 ++local_ind;
66 }
67 indset_rw.endResize();
68 assert(indset_rw.size() == indset.size() + nw);
69 // assume same communication pattern
70 commRW->remoteIndices().setNeighbours(comm.remoteIndices().getNeighbours());
71 commRW->remoteIndices().template rebuild<true>();
72 //commRW->clearInterfaces(); may need this for correct rebuild
73 }
74
75 namespace Details
76 {
77 using PressureMatrixType = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
78 using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
79 using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
80 template <class Comm>
81 using ParCoarseOperatorType
82 = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Comm>;
83 template <class Comm>
84 using CoarseOperatorType = std::conditional_t<std::is_same<Comm, Dune::Amg::SequentialInformation>::value,
85 SeqCoarseOperatorType,
86 ParCoarseOperatorType<Comm>>;
87 } // namespace Details
88
89 template <class FineOperator, class Communication, bool transpose = false>
90 class PressureBhpTransferPolicy : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
91 {
92 public:
93 typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
95 typedef Communication ParallelInformation;
96 typedef typename FineOperator::domain_type FineVectorType;
97
98 public:
99 PressureBhpTransferPolicy(const Communication& comm,
100 const FineVectorType& weights,
101 const Opm::PropertyTree& prm,
102 const std::size_t pressureIndex)
103 : communication_(&const_cast<Communication&>(comm))
104 , weights_(weights)
105 , prm_(prm)
106 , pressure_var_index_(pressureIndex)
107 {
108 }
109
111 {
113 using CoarseMatrix = typename CoarseOperator::matrix_type;
114 const auto& fineLevelMatrix = fineOperator.getmat();
115 const auto& nw = fineOperator.getNumberOfExtraEquations();
116 if (prm_.get<bool>("add_wells")) {
117 const size_t average_elements_per_row
118 = static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
119 const double overflow_fraction = 1.2;
120 coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N() + nw,
121 fineLevelMatrix.M() + nw,
124 CoarseMatrix::implicit));
125 int rownum = 0;
126 for (const auto& row : fineLevelMatrix) {
127 for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
128 coarseLevelMatrix_->entry(rownum, col.index()) = 0.0;
129 }
130 ++rownum;
131 }
132 } else {
133 coarseLevelMatrix_.reset(
134 new CoarseMatrix(fineLevelMatrix.N(), fineLevelMatrix.M(), CoarseMatrix::row_wise));
135 auto createIter = coarseLevelMatrix_->createbegin();
136 for (const auto& row : fineLevelMatrix) {
137 for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
138 createIter.insert(col.index());
139 }
140 ++createIter;
141 }
142 }
143 if constexpr (std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
144 coarseLevelCommunication_ = std::make_shared<Communication>();
145 } else {
146 coarseLevelCommunication_ = std::make_shared<Communication>(
147 communication_->communicator(), communication_->category(), false);
148 }
149 if (prm_.get<bool>("add_wells")) {
150 fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
151 coarseLevelMatrix_->compress(); // all elemenst should be set
152 if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
153 extendCommunicatorWithWells(*communication_, coarseLevelCommunication_, nw);
154 }
155 }
157
158 this->lhs_.resize(this->coarseLevelMatrix_->M());
159 this->rhs_.resize(this->coarseLevelMatrix_->N());
160 using OperatorArgs = typename Dune::Amg::ConstructionTraits<CoarseOperator>::Arguments;
161 OperatorArgs oargs(coarseLevelMatrix_, *coarseLevelCommunication_);
162 this->operator_ = Dune::Amg::ConstructionTraits<CoarseOperator>::construct(oargs);
163 }
164
166 {
168 const auto& fineMatrix = fineOperator.getmat();
169 *coarseLevelMatrix_ = 0;
170 auto rowCoarse = coarseLevelMatrix_->begin();
171 for (auto row = fineMatrix.begin(), rowEnd = fineMatrix.end(); row != rowEnd; ++row, ++rowCoarse) {
172 assert(row.index() == rowCoarse.index());
173 auto entryCoarse = rowCoarse->begin();
174 for (auto entry = row->begin(), entryEnd = row->end(); entry != entryEnd; ++entry, ++entryCoarse) {
175 assert(entry.index() == entryCoarse.index());
176 double matrix_el = 0;
177 if (transpose) {
178 const auto& bw = weights_[entry.index()];
179 for (size_t i = 0; i < bw.size(); ++i) {
180 matrix_el += (*entry)[pressure_var_index_][i] * bw[i];
181 }
182 } else {
183 const auto& bw = weights_[row.index()];
184 for (size_t i = 0; i < bw.size(); ++i) {
185 matrix_el += (*entry)[i][pressure_var_index_] * bw[i];
186 }
187 }
188 (*entryCoarse) = matrix_el;
189 }
190 }
191 if (prm_.get<bool>("add_wells")) {
193 assert(transpose == false); // not implemented
194 bool use_well_weights = prm_.get<bool>("use_well_weights");
195 fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
196#ifndef NDEBUG
197 std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
198 assert(rowCoarse == coarseLevelMatrix_->end());
199#endif
200
201 }
202 }
203
204 virtual void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override
205 {
206 OPM_TIMEBLOCK(moveToCoarseLevel);
207 //NB we iterate over fine assumming welldofs is at the end
208 // Set coarse vector to zero
209 this->rhs_ = 0;
210
211 auto end = fine.end(), begin = fine.begin();
212
213 for (auto block = begin; block != end; ++block) {
214 const auto& bw = weights_[block.index()];
215 double rhs_el = 0.0;
216 if (transpose) {
217 rhs_el = (*block)[pressure_var_index_];
218 } else {
219 for (size_t i = 0; i < block->size(); ++i) {
220 rhs_el += (*block)[i] * bw[i];
221 }
222 }
223 this->rhs_[block - begin] = rhs_el;
224 }
225
226 this->lhs_ = 0;
227 }
228
229 virtual void moveToFineLevel(typename ParentType::FineDomainType& fine) override
230 {
232 //NB we iterate over fine assumming welldofs is at the end
233 auto end = fine.end(), begin = fine.begin();
234
235 for (auto block = begin; block != end; ++block) {
236 if (transpose) {
237 const auto& bw = weights_[block.index()];
238 for (size_t i = 0; i < block->size(); ++i) {
239 (*block)[i] = this->lhs_[block - begin] * bw[i];
240 }
241 } else {
242 (*block)[pressure_var_index_] = this->lhs_[block - begin];
243 }
244 }
245 }
246
247 virtual PressureBhpTransferPolicy* clone() const override
248 {
249 return new PressureBhpTransferPolicy(*this);
250 }
251
252 const Communication& getCoarseLevelCommunication() const
253 {
254 return *coarseLevelCommunication_;
255 }
256
257private:
258 Communication* communication_;
259 const FineVectorType& weights_;
260 PropertyTree prm_;
261 const int pressure_var_index_;
262 std::shared_ptr<Communication> coarseLevelCommunication_;
263 std::shared_ptr<typename CoarseOperator::matrix_type> coarseLevelMatrix_;
264};
265
266} // namespace Opm
267
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 AquiferInterface.hpp:35
Definition PressureBhpTransferPolicy.hpp:91
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:229
virtual void calculateCoarseEntries(const FineOperator &fineOperator) override
???.
Definition PressureBhpTransferPolicy.hpp:165
virtual PressureBhpTransferPolicy * clone() const override
Clone the current object.
Definition PressureBhpTransferPolicy.hpp:247
virtual void createCoarseLevelSystem(const FineOperator &fineOperator) override
Algebraically creates the coarse level system.
Definition PressureBhpTransferPolicy.hpp:110
Definition PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Algebraic twolevel methods.