My Project
FlexibleSolver_impl.hpp
1 /*
2  Copyright 2019, 2020 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2020 Equinor.
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_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
22 #define OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
23 
24 #include <opm/simulators/linalg/matrixblock.hh>
25 #include <opm/simulators/linalg/ilufirstelement.hh>
26 #include <opm/simulators/linalg/FlexibleSolver.hpp>
27 #include <opm/simulators/linalg/PreconditionerFactory.hpp>
28 
29 #include <dune/common/fmatrix.hh>
30 #include <dune/istl/bcrsmatrix.hh>
31 #include <dune/istl/solvers.hh>
32 #include <dune/istl/umfpack.hh>
33 #include <dune/istl/owneroverlapcopy.hh>
34 #include <dune/istl/paamg/pinfo.hh>
35 
36 namespace Dune
37 {
39  template <class MatrixType, class VectorType>
42  const Opm::PropertyTree& prm,
43  const std::function<VectorType()>& weightsCalculator,
44  std::size_t pressureIndex)
45  {
46  init(op, Dune::Amg::SequentialInformation(), prm, weightsCalculator,
47  pressureIndex);
48  }
49 
51  template <class MatrixType, class VectorType>
52  template <class Comm>
55  const Comm& comm,
56  const Opm::PropertyTree& prm,
57  const std::function<VectorType()>& weightsCalculator,
58  std::size_t pressureIndex)
59  {
60  init(op, comm, prm, weightsCalculator, pressureIndex);
61  }
62 
63  template <class MatrixType, class VectorType>
64  void
66  apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res)
67  {
68  linsolver_->apply(x, rhs, res);
69  }
70 
71  template <class MatrixType, class VectorType>
72  void
73  FlexibleSolver<MatrixType, VectorType>::
74  apply(VectorType& x, VectorType& rhs, double reduction, Dune::InverseOperatorResult& res)
75  {
76  linsolver_->apply(x, rhs, reduction, res);
77  }
78 
80  template <class MatrixType, class VectorType>
81  auto
84  {
85  return *preconditioner_;
86  }
87 
88  template <class MatrixType, class VectorType>
89  Dune::SolverCategory::Category
91  category() const
92  {
93  return linearoperator_for_solver_->category();
94  }
95 
96  // Machinery for making sequential or parallel operators/preconditioners/scalar products.
97  template <class MatrixType, class VectorType>
98  template <class Comm>
99  void
100  FlexibleSolver<MatrixType, VectorType>::
101  initOpPrecSp(AbstractOperatorType& op,
102  const Opm::PropertyTree& prm,
103  const std::function<VectorType()> weightsCalculator,
104  const Comm& comm,
105  std::size_t pressureIndex)
106  {
107  // Parallel case.
108  using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
109  linearoperator_for_solver_ = &op;
110  auto op_prec = std::make_shared<ParOperatorType>(op.getmat(), comm);
111  auto child = prm.get_child_optional("preconditioner");
113  child ? *child : Opm::PropertyTree(),
114  weightsCalculator,
115  comm,
116  pressureIndex);
117  scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, op.category());
118  linearoperator_for_precond_ = op_prec;
119  }
120 
121  template <class MatrixType, class VectorType>
122  void
123  FlexibleSolver<MatrixType, VectorType>::
124  initOpPrecSp(AbstractOperatorType& op,
125  const Opm::PropertyTree& prm,
126  const std::function<VectorType()> weightsCalculator,
127  const Dune::Amg::SequentialInformation&,
128  std::size_t pressureIndex)
129  {
130  // Sequential case.
131  using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
132  linearoperator_for_solver_ = &op;
133  auto op_prec = std::make_shared<SeqOperatorType>(op.getmat());
134  auto child = prm.get_child_optional("preconditioner");
135  preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*op_prec,
136  child ? *child : Opm::PropertyTree(),
137  weightsCalculator,
138  pressureIndex);
139  scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
140  linearoperator_for_precond_ = op_prec;
141  }
142 
143  template <class MatrixType, class VectorType>
144  void
145  FlexibleSolver<MatrixType, VectorType>::
146  initSolver(const Opm::PropertyTree& prm, const bool is_iorank)
147  {
148  const double tol = prm.get<double>("tol", 1e-2);
149  const int maxiter = prm.get<int>("maxiter", 200);
150  const int verbosity = is_iorank ? prm.get<int>("verbosity", 0) : 0;
151  const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
152  if (solver_type == "bicgstab") {
153  linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_for_solver_,
154  *scalarproduct_,
155  *preconditioner_,
156  tol, // desired residual reduction factor
157  maxiter, // maximum number of iterations
158  verbosity));
159  } else if (solver_type == "loopsolver") {
160  linsolver_.reset(new Dune::LoopSolver<VectorType>(*linearoperator_for_solver_,
161  *scalarproduct_,
162  *preconditioner_,
163  tol, // desired residual reduction factor
164  maxiter, // maximum number of iterations
165  verbosity));
166  } else if (solver_type == "gmres") {
167  int restart = prm.get<int>("restart", 15);
168  linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_for_solver_,
169  *scalarproduct_,
170  *preconditioner_,
171  tol,
172  restart, // desired residual reduction factor
173  maxiter, // maximum number of iterations
174  verbosity));
175 #if HAVE_SUITESPARSE_UMFPACK
176  } else if (solver_type == "umfpack") {
177  bool dummy = false;
178  linsolver_.reset(new Dune::UMFPack<MatrixType>(linearoperator_for_solver_->getmat(), verbosity, dummy));
179 #endif
180  } else {
181  OPM_THROW(std::invalid_argument, "Properties: Solver " << solver_type << " not known.");
182  }
183  }
184 
185 
186  // Main initialization routine.
187  // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
188  template <class MatrixType, class VectorType>
189  template <class Comm>
190  void
191  FlexibleSolver<MatrixType, VectorType>::
192  init(AbstractOperatorType& op,
193  const Comm& comm,
194  const Opm::PropertyTree& prm,
195  const std::function<VectorType()> weightsCalculator,
196  std::size_t pressureIndex)
197  {
198  initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
199  initSolver(prm, comm.communicator().rank() == 0);
200  }
201 
202 } // namespace Dune
203 
204 
205 // Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
206 
207 template <int N>
208 using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
209 template <int N>
210 using BM = Dune::BCRSMatrix<Dune::FieldMatrix<double, N, N>>;
211 template <int N>
212 using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
213 
214 #if HAVE_MPI
215 
216 using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
217 
218 // Note: we must instantiate the constructor that is a template.
219 // This is only needed in the parallel case, since otherwise the Comm type is
220 // not a template argument but always SequentialInformation.
221 
222 #define INSTANTIATE_FLEXIBLESOLVER(N) \
223 template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
224 template class Dune::FlexibleSolver<OBM<N>, BV<N>>; \
225 template Dune::FlexibleSolver<BM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
226  const Comm& comm, \
227  const Opm::PropertyTree& prm, \
228  const std::function<BV<N>()>& weightsCalculator, \
229 std::size_t); \
230 template Dune::FlexibleSolver<OBM<N>, BV<N>>::FlexibleSolver(AbstractOperatorType& op, \
231  const Comm& comm, \
232  const Opm::PropertyTree& prm, \
233  const std::function<BV<N>()>& weightsCalculator, \
234 std::size_t);
235 
236 #else // HAVE_MPI
237 
238 #define INSTANTIATE_FLEXIBLESOLVER(N) \
239 template class Dune::FlexibleSolver<BM<N>, BV<N>>; \
240 template class Dune::FlexibleSolver<OBM<N>, BV<N>>;
241 
242 #endif // HAVE_MPI
243 
244 
245 #endif // OPM_FLEXIBLE_SOLVER_IMPL_HEADER_INCLUDED
A solver class that encapsulates all needed objects for a linear solver (operator,...
Definition: FlexibleSolver.hpp:40
Dune::AssembledLinearOperator< MatrixType, VectorType, VectorType > AbstractOperatorType
Base class type of the operator passed to the solver.
Definition: FlexibleSolver.hpp:46
FlexibleSolver(AbstractOperatorType &op, const Opm::PropertyTree &prm, const std::function< VectorType()> &weightsCalculator, std::size_t pressureIndex)
Create a sequential solver.
Definition: FlexibleSolver_impl.hpp:41
AbstractPrecondType & preconditioner()
Access the contained preconditioner.
Definition: FlexibleSolver_impl.hpp:83
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: PropertyTree.hpp:37