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>
25 #include <opm/simulators/linalg/PressureTransferPolicy.hpp>
26 #include <opm/simulators/linalg/getQuasiImpesWeights.hpp>
28 
29 #include <opm/common/ErrorMacros.hpp>
30 
31 #include <dune/common/fmatrix.hh>
32 #include <dune/istl/bcrsmatrix.hh>
33 #include <dune/istl/paamg/amg.hh>
34 
35 #include <fstream>
36 #include <type_traits>
37 
38 
39 namespace Opm
40 {
41 // Circular dependency between PreconditionerFactory [which can make an OwningTwoLevelPreconditioner]
42 // and OwningTwoLevelPreconditioner [which uses PreconditionerFactory to choose the fine-level smoother]
43 // must be broken, accomplished by forward-declaration here.
44 template <class Operator, class Comm = Dune::Amg::SequentialInformation>
45 class PreconditionerFactory;
46 }
47 
48 namespace Dune
49 {
50 
51 
52 // Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
53 template <class MatrixTypeT, class VectorTypeT>
54 class FlexibleSolver;
55 
56 template <typename T, typename A, int i>
57 std::ostream& operator<<(std::ostream& out,
58  const BlockVector< FieldVector< T, i >, A >& vector)
59 {
60  Dune::writeMatrixMarket(vector, out);
61  return out;
62 }
63 
64  template <typename T, typename A, int i>
65  std::istream& operator>>(std::istream& input,
66  BlockVector< FieldVector< T, i >, A >& vector)
67 {
68  Dune::readMatrixMarket(vector, input);
69  return input;
70 }
71 
76 template <class OperatorType,
77  class VectorType,
78  bool transpose = false,
79  class Communication = Dune::Amg::SequentialInformation>
80 class OwningTwoLevelPreconditioner : public Dune::PreconditionerWithUpdate<VectorType, VectorType>
81 {
82 public:
83  using MatrixType = typename OperatorType::matrix_type;
85 
86  OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
87  const std::function<VectorType()> weightsCalculator,
88  std::size_t pressureIndex)
89  : linear_operator_(linearoperator)
90  , finesmoother_(PrecFactory::create(linearoperator,
91  prm.get_child_optional("finesmoother") ?
92  prm.get_child("finesmoother") : Opm::PropertyTree(),
93  std::function<VectorType()>(), pressureIndex))
94  , comm_(nullptr)
95  , weightsCalculator_(weightsCalculator)
96  , weights_(weightsCalculator())
97  , levelTransferPolicy_(dummy_comm_, weights_, pressureIndex)
98  , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree())
99  , twolevel_method_(linearoperator,
100  finesmoother_,
101  levelTransferPolicy_,
102  coarseSolverPolicy_,
103  prm.get<int>("pre_smooth", transpose? 1 : 0),
104  prm.get<int>("post_smooth", transpose? 0 : 1))
105  , prm_(prm)
106  {
107  if (prm.get<int>("verbosity", 0) > 10) {
108  std::string filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
109  std::ofstream outfile(filename);
110  if (!outfile) {
111  OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
112  }
113  Dune::writeMatrixMarket(weights_, outfile);
114  }
115  }
116 
117  OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const Opm::PropertyTree& prm,
118  const std::function<VectorType()> weightsCalculator,
119  std::size_t pressureIndex, const Communication& comm)
120  : linear_operator_(linearoperator)
121  , finesmoother_(PrecFactory::create(linearoperator,
122  prm.get_child_optional("finesmoother") ?
123  prm.get_child("finesmoother"): Opm::PropertyTree(),
124  std::function<VectorType()>(),
125  comm, pressureIndex))
126  , comm_(&comm)
127  , weightsCalculator_(weightsCalculator)
128  , weights_(weightsCalculator())
129  , levelTransferPolicy_(*comm_, weights_, pressureIndex)
130  , coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : Opm::PropertyTree())
131  , twolevel_method_(linearoperator,
132  finesmoother_,
133  levelTransferPolicy_,
134  coarseSolverPolicy_,
135  prm.get<int>("pre_smooth", transpose? 1 : 0),
136  prm.get<int>("post_smooth", transpose? 0 : 1))
137  , prm_(prm)
138  {
139  if (prm.get<int>("verbosity", 0) > 10 && comm.communicator().rank() == 0) {
140  auto filename = prm.get<std::string>("weights_filename", "impes_weights.txt");
141  std::ofstream outfile(filename);
142  if (!outfile) {
143  OPM_THROW(std::ofstream::failure, "Could not write weights to file " << filename << ".");
144  }
145  Dune::writeMatrixMarket(weights_, outfile);
146  }
147  }
148 
149  virtual void pre(VectorType& x, VectorType& b) override
150  {
151  twolevel_method_.pre(x, b);
152  }
153 
154  virtual void apply(VectorType& v, const VectorType& d) override
155  {
156  twolevel_method_.apply(v, d);
157  }
158 
159  virtual void post(VectorType& x) override
160  {
161  twolevel_method_.post(x);
162  }
163 
164  virtual void update() override
165  {
166  weights_ = weightsCalculator_();
167  updateImpl(comm_);
168  }
169 
170  virtual Dune::SolverCategory::Category category() const override
171  {
172  return linear_operator_.category();
173  }
174 
175 private:
176  using PressureMatrixType = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
177  using PressureVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
178  using SeqCoarseOperatorType = Dune::MatrixAdapter<PressureMatrixType, PressureVectorType, PressureVectorType>;
179  using ParCoarseOperatorType
180  = Dune::OverlappingSchwarzOperator<PressureMatrixType, PressureVectorType, PressureVectorType, Communication>;
181  using CoarseOperatorType = std::conditional_t<std::is_same<Communication, Dune::Amg::SequentialInformation>::value,
182  SeqCoarseOperatorType,
183  ParCoarseOperatorType>;
185  using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy<CoarseOperatorType,
188  using TwoLevelMethod
190 
191  // Handling parallel vs serial instantiation of preconditioner factory.
192  template <class Comm>
193  void updateImpl(const Comm*)
194  {
195  // Parallel case.
196  auto child = prm_.get_child_optional("finesmoother");
197  finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree(), *comm_);
198  twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
199  }
200 
201  void updateImpl(const Dune::Amg::SequentialInformation*)
202  {
203  // Serial case.
204  auto child = prm_.get_child_optional("finesmoother");
205  finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : Opm::PropertyTree());
206  twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
207  }
208 
209  const OperatorType& linear_operator_;
210  std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
211  const Communication* comm_;
212  std::function<VectorType()> weightsCalculator_;
213  VectorType weights_;
214  LevelTransferPolicy levelTransferPolicy_;
215  CoarseSolverPolicy coarseSolverPolicy_;
216  TwoLevelMethod twolevel_method_;
217  Opm::PropertyTree prm_;
218  Communication dummy_comm_;
219 };
220 
221 } // namespace Dune
222 
223 
224 
225 
226 #endif // OPM_OWNINGTWOLEVELPRECONDITIONER_HEADER_INCLUDED
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
A version of the two-level preconditioner that is:
Definition: OwningTwoLevelPreconditioner.hpp:81
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:51
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.hpp:71
Definition: PressureTransferPolicy.hpp:33
Definition: PropertyTree.hpp:37
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Algebraic twolevel methods.