OwningTwoLevelPreconditioner.hpp
Go to the documentation of this file.
1/*
2 Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20#ifndef OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
21#define OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
22
26
27#include <opm/common/ErrorMacros.hpp>
28
29#include <dune/common/fmatrix.hh>
30#include <dune/istl/bcrsmatrix.hh>
31#include <dune/istl/paamg/amg.hh>
32
33#include <fstream>
34#include <type_traits>
35
36
37namespace Opm
38{
39// Circular dependency between PreconditionerFactory [which can make an OwningTwoLevelPreconditioner]
40// and OwningTwoLevelPreconditioner [which uses PreconditionerFactory to choose the fine-level smoother]
41// must be broken, accomplished by forward-declaration here.
42template <class Operator, class Comm = Dune::Amg::SequentialInformation>
43class PreconditionerFactory;
44}
45
46namespace Dune
47{
48
49
50// Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
51template <class Operator>
52class FlexibleSolver;
53
54template <typename T, typename A, int i>
55std::ostream& operator<<(std::ostream& out,
56 const BlockVector< FieldVector< T, i >, A >& vector)
57{
58 Dune::writeMatrixMarket(vector, out);
59 return out;
60}
61
62 template <typename T, typename A, int i>
63 std::istream& operator>>(std::istream& input,
64 BlockVector< FieldVector< T, i >, A >& vector)
65{
66 Dune::readMatrixMarket(vector, input);
67 return input;
68}
69
74template <class OperatorType,
75 class VectorType,
76 class LevelTransferPolicy,
77 class Communication = Dune::Amg::SequentialInformation>
79{
80public:
81 using MatrixType = typename OperatorType::matrix_type;
83 using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
84
85 OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
86 const std::function<VectorType()> weightsCalculator,
87 std::size_t pressureIndex)
88 : linear_operator_(linearoperator)
89 , finesmoother_(PrecFactory::create(linearoperator,
90 prm.get_child_optional("finesmoother") ?
91 prm.get_child("finesmoother") : Opm::PropertyTree(),
92 std::function<VectorType()>(), pressureIndex))
93 , comm_(nullptr)
94 , weightsCalculator_(weightsCalculator)
95 , weights_(weightsCalculator())
96 , levelTransferPolicy_(dummy_comm_, weights_, prm, pressureIndex)
97 , coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
98 , twolevel_method_(linearoperator,
99 finesmoother_,
100 levelTransferPolicy_,
101 coarseSolverPolicy_,
102 prm.get<int>("pre_smooth", 0),
103 prm.get<int>("post_smooth", 1))
104 , prm_(prm)
105 {
106 if (prm.get<int>("verbosity", 0) > 10) {
107 std::string filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
108 std::ofstream outfile(filename);
109 if (!outfile) {
110 OPM_THROW(std::ofstream::failure,
111 "Could not write weights to file " + filename + ".");
112 }
113 // Unqualified on purpose to enable ADL for GPU/CPU types
114 writeMatrixMarket(weights_, outfile);
115 }
116 }
117
118 OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
119 const std::function<VectorType()> weightsCalculator,
120 std::size_t pressureIndex, const Communication& comm)
121 : linear_operator_(linearoperator)
122 , finesmoother_(PrecFactory::create(linearoperator,
123 prm.get_child_optional("finesmoother") ?
124 prm.get_child("finesmoother"): Opm::PropertyTree(),
125 std::function<VectorType()>(),
126 comm, pressureIndex))
127 , comm_(&comm)
128 , weightsCalculator_(weightsCalculator)
129 , weights_(weightsCalculator())
130 , levelTransferPolicy_(*comm_, weights_, prm, pressureIndex)
131 , coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
132 , twolevel_method_(linearoperator,
133 finesmoother_,
134 levelTransferPolicy_,
135 coarseSolverPolicy_,
136 prm.get<int>("pre_smooth", 0),
137 prm.get<int>("post_smooth", 1))
138 , prm_(prm)
139 {
140 if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
141 auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
142 std::ofstream outfile(filename);
143 if (!outfile) {
144 OPM_THROW(std::ofstream::failure,
145 "Could not write weights to file " + filename + ".");
146 }
147 // Unqualified on purpose to enable ADL for GPU/CPU types
148 writeMatrixMarket(weights_, outfile);
149 }
150 }
151
152 virtual void pre(VectorType& x, VectorType& b) override
153 {
154 twolevel_method_.pre(x, b);
155 }
156
157 virtual void apply(VectorType& v, const VectorType& d) override
158 {
159 twolevel_method_.apply(v, d);
160 }
161
162 virtual void post(VectorType& x) override
163 {
164 twolevel_method_.post(x);
165 }
166
167 virtual void update() override
168 {
169 weights_ = weightsCalculator_();
170 updateImpl(comm_);
171 }
172
173 virtual Dune::SolverCategory::Category category() const override
174 {
175 return linear_operator_.category();
176 }
177
178 virtual bool hasPerfectUpdate() const override {
179 return twolevel_method_.hasPerfectUpdate();
180 }
181
182private:
183 using CoarseOperator = typename LevelTransferPolicy::CoarseOperator;
184 using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperator,
186 LevelTransferPolicy>;
187 using TwoLevelMethod
189
190 // Handling parallel vs serial instantiation of preconditioner factory.
191 template <class Comm>
192 void updateImpl(const Comm*)
193 {
194 // Parallel case.
195 twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
196 }
197
198 void updateImpl(const Dune::Amg::SequentialInformation*)
199 {
200 // Serial case.
201 twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
202 }
203
204 const OperatorType& linear_operator_;
205 std::shared_ptr<PreconditionerWithUpdate<VectorType, VectorType>> finesmoother_;
206 const Communication* comm_;
207 std::function<VectorType()> weightsCalculator_;
208 VectorType weights_;
209 LevelTransferPolicy levelTransferPolicy_;
210 CoarseSolverPolicy coarseSolverPolicy_;
211 TwoLevelMethod twolevel_method_;
213 Communication dummy_comm_;
214};
215
216} // namespace Dune
217
218
219
220
221#endif // OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
Dune::OwnerOverlapCopyCommunication< int, int > Comm
Definition: FlexibleSolver_impl.hpp:304
Definition: PressureSolverPolicy.hpp:35
Definition: twolevelmethodcpr.hh:387
void updatePreconditioner(FineOperatorType &, std::shared_ptr< SmootherType > smoother, CoarseLevelSolverPolicy &coarsePolicy)
Definition: twolevelmethodcpr.hh:463
void apply(FineDomainType &v, const FineRangeType &d)
Definition: twolevelmethodcpr.hh:508
void pre(FineDomainType &x, FineRangeType &b)
Definition: twolevelmethodcpr.hh:487
void post(FineDomainType &x)
Definition: twolevelmethodcpr.hh:498
bool hasPerfectUpdate() const
Definition: twolevelmethodcpr.hh:502
Definition: FlexibleSolver.hpp:45
Definition: OwningTwoLevelPreconditioner.hpp:79
virtual void apply(VectorType &v, const VectorType &d) override
Definition: OwningTwoLevelPreconditioner.hpp:157
virtual bool hasPerfectUpdate() const override
Definition: OwningTwoLevelPreconditioner.hpp:178
typename OperatorType::matrix_type MatrixType
Definition: OwningTwoLevelPreconditioner.hpp:81
OwningTwoLevelPreconditioner(const OperatorType &linearoperator, const Opm::PropertyTree &prm, const std::function< VectorType()> weightsCalculator, std::size_t pressureIndex)
Definition: OwningTwoLevelPreconditioner.hpp:85
virtual void post(VectorType &x) override
Definition: OwningTwoLevelPreconditioner.hpp:162
virtual void pre(VectorType &x, VectorType &b) override
Definition: OwningTwoLevelPreconditioner.hpp:152
OwningTwoLevelPreconditioner(const OperatorType &linearoperator, const Opm::PropertyTree &prm, const std::function< VectorType()> weightsCalculator, std::size_t pressureIndex, const Communication &comm)
Definition: OwningTwoLevelPreconditioner.hpp:118
virtual void update() override
Definition: OwningTwoLevelPreconditioner.hpp:167
Dune::AssembledLinearOperator< MatrixType, VectorType, VectorType > AbstractOperatorType
Definition: OwningTwoLevelPreconditioner.hpp:83
virtual Dune::SolverCategory::Category category() const override
Definition: OwningTwoLevelPreconditioner.hpp:173
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
Definition: PreconditionerFactory.hpp:64
Hierarchical collection of key/value pairs.
Definition: PropertyTree.hpp:39
T get(const std::string &key) const
Definition: fvbaseprimaryvariables.hh:141
std::istream & operator>>(std::istream &input, BlockVector< FieldVector< T, i >, A > &vector)
Definition: OwningTwoLevelPreconditioner.hpp:63
std::ostream & operator<<(std::ostream &out, const BlockVector< FieldVector< T, i >, A > &vector)
Definition: OwningTwoLevelPreconditioner.hpp:55
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
void writeMatrixMarket(const GpuVector< T > &vectorOnDevice, std::ostream &ostr)
Definition: GpuVector.hpp:538
Definition: blackoilbioeffectsmodules.hh:43
Algebraic twolevel methods.