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
30namespace Opm
31{
32
33namespace 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
49template <class FineOperator, class Communication, bool transpose = false>
51 : public Dune::Amg::LevelTransferPolicyCpr<FineOperator, Details::CoarseOperatorType<Communication>>
52{
53public:
54 typedef typename Details::CoarseOperatorType<Communication> CoarseOperator;
56 typedef Communication ParallelInformation;
57 typedef typename FineOperator::domain_type FineVectorType;
58
59public:
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 }
180private:
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 PressureTransferPolicy * clone() const override
Clone the current object.
Definition: PressureTransferPolicy.hpp:166
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 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.