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