My Project
OwningTwoLevelPreconditioner.hpp
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
23#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
24#include <opm/simulators/linalg/PressureSolverPolicy.hpp>
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, "Could not write weights to file " << filename << ".");
111 }
112 Dune::writeMatrixMarket(weights_, outfile);
113 }
114 }
115
116 OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
117 const std::function<VectorType()> weightsCalculator,
118 std::size_t pressureIndex, const Communication& comm)
119 : linear_operator_(linearoperator)
120 , finesmoother_(PrecFactory::create(linearoperator,
121 prm.get_child_optional("finesmoother") ?
122 prm.get_child("finesmoother"): Opm::PropertyTree(),
123 std::function<VectorType()>(),
124 comm, pressureIndex))
125 , comm_(&comm)
126 , weightsCalculator_(weightsCalculator)
127 , weights_(weightsCalculator())
128 , levelTransferPolicy_(*comm_, weights_, prm, pressureIndex)
129 , coarseSolverPolicy_(prm.get_child_optional("coarsesolver") ? prm.get_child("coarsesolver") : Opm::PropertyTree())
130 , twolevel_method_(linearoperator,
131 finesmoother_,
132 levelTransferPolicy_,
133 coarseSolverPolicy_,
134 prm.get<int>("pre_smooth", 0),
135 prm.get<int>("post_smooth", 1))
136 , prm_(prm)
137 {
138 if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
139 auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
140 std::ofstream outfile(filename);
141 if (!outfile) {
142 OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
143 }
144 Dune::writeMatrixMarket(weights_, outfile);
145 }
146 }
147
148 virtual void pre(VectorType& x, VectorType& b) override
149 {
150 twolevel_method_.pre(x, b);
151 }
152
153 virtual void apply(VectorType& v, const VectorType& d) override
154 {
155 twolevel_method_.apply(v, d);
156 }
157
158 virtual void post(VectorType& x) override
159 {
160 twolevel_method_.post(x);
161 }
162
163 virtual void update() override
164 {
165 weights_ = weightsCalculator_();
166 updateImpl(comm_);
167 }
168
169 virtual Dune::SolverCategory::Category category() const override
170 {
171 return linear_operator_.category();
172 }
173
174private:
175 using CoarseOperator = typename LevelTransferPolicy::CoarseOperator;
178 LevelTransferPolicy>;
179 using TwoLevelMethod
181
182 // Handling parallel vs serial instantiation of preconditioner factory.
183 template <class Comm>
184 void updateImpl(const Comm*)
185 {
186 // Parallel case.
187 auto child = prm_.get_child_optional("finesmoother");
188 finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree(), *comm_);
189 twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
190 }
191
192 void updateImpl(const Dune::Amg::SequentialInformation*)
193 {
194 // Serial case.
195 auto child = prm_.get_child_optional("finesmoother");
196 finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree());
197 twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
198 }
199
200 const OperatorType& linear_operator_;
201 std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
202 const Communication* comm_;
203 std::function<VectorType()> weightsCalculator_;
204 VectorType weights_;
205 LevelTransferPolicy levelTransferPolicy_;
206 CoarseSolverPolicy coarseSolverPolicy_;
207 TwoLevelMethod twolevel_method_;
209 Communication dummy_comm_;
210};
211
212} // namespace Dune
213
214
215
216
217#endif // OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:45
A version of the two-level preconditioner that is:
Definition: OwningTwoLevelPreconditioner.hpp:79
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
This is an object factory for creating preconditioners.
Definition: PreconditionerFactory.hpp:69
static PrecPtr create(const Operator &op, const PropertyTree &prm, const std::function< Vector()> &weightsCalculator={}, std::size_t pressureIndex=std::numeric_limits< std::size_t >::max())
Create a new serial preconditioner and return a pointer to it.
Definition: PreconditionerFactory_impl.hpp:500
Definition: PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Algebraic twolevel methods.