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 
37 namespace 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.
42 template <class Operator, class Comm = Dune::Amg::SequentialInformation>
43 class PreconditionerFactory;
44 }
45 
46 namespace Dune
47 {
48 
49 
50 // Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
51 template <class Operator>
52 class FlexibleSolver;
53 
54 template <typename T, typename A, int i>
55 std::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 
74 template <class OperatorType,
75  class VectorType,
76  class LevelTransferPolicy,
77  class Communication = Dune::Amg::SequentialInformation>
78 class OwningTwoLevelPreconditioner : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
79 {
80 public:
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 
174 private:
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_;
208  Opm::PropertyTree prm_;
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.