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
40namespace 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_ = std::make_shared<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_ = std::make_shared<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_ = std::make_shared<Dune::RestartedGMResSolver<VectorType>>(*linearoperator_for_solver_,
167 *scalarproduct_,
168 *preconditioner_,
169 tol,// desired residual reduction factor
170 restart,
171 maxiter, // maximum number of iterations
172 verbosity);
173 } else if (solver_type == "flexgmres") {
174 int restart = prm.get<int>("restart", 15);
175 linsolver_ = std::make_shared<Dune::RestartedFlexibleGMResSolver<VectorType>>(*linearoperator_for_solver_,
176 *scalarproduct_,
177 *preconditioner_,
178 tol,// desired residual reduction factor
179 restart,
180 maxiter, // maximum number of iterations
181 verbosity);
182#if HAVE_SUITESPARSE_UMFPACK
183 } else if (solver_type == "umfpack") {
184 bool dummy = false;
185 using MatrixType = std::remove_const_t<std::remove_reference_t<decltype(linearoperator_for_solver_->getmat())>>;
186 linsolver_ = std::make_shared<Dune::UMFPack<MatrixType>>(linearoperator_for_solver_->getmat(), verbosity, dummy);
187#endif
188 } else {
189 OPM_THROW(std::invalid_argument,
190 "Properties: Solver " + solver_type + " not known.");
191 }
192 }
193
194
195 // Main initialization routine.
196 // Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
197 template <class Operator>
198 template <class Comm>
199 void
200 FlexibleSolver<Operator>::
201 init(Operator& op,
202 const Comm& comm,
203 const Opm::PropertyTree& prm,
204 const std::function<VectorType()> weightsCalculator,
205 std::size_t pressureIndex)
206 {
207 initOpPrecSp(op, prm, weightsCalculator, comm, pressureIndex);
208 initSolver(prm, comm.communicator().rank() == 0);
209 }
210
211} // namespace Dune
212
213
214// Macros to simplify explicit instantiation of FlexibleSolver for various block sizes.
215
216// Vectors and matrices.
217template <int N>
218using BV = Dune::BlockVector<Dune::FieldVector<double, N>>;
219template <int N>
220using OBM = Dune::BCRSMatrix<Opm::MatrixBlock<double, N, N>>;
221
222// Sequential operators.
223template <int N>
224using SeqOpM = Dune::MatrixAdapter<OBM<N>, BV<N>, BV<N>>;
225template <int N>
226using SeqOpW = Opm::WellModelMatrixAdapter<OBM<N>, BV<N>, BV<N>, false>;
227
228#if HAVE_MPI
229
230// Parallel communicator and operators.
231using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
232template <int N>
233using ParOpM = Dune::OverlappingSchwarzOperator<OBM<N>, BV<N>, BV<N>, Comm>;
234template <int N>
235using ParOpW = Opm::WellModelGhostLastMatrixAdapter<OBM<N>, BV<N>, BV<N>, true>;
236
237// Note: we must instantiate the constructor that is a template.
238// This is only needed in the parallel case, since otherwise the Comm type is
239// not a template argument but always SequentialInformation.
240
241#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
242template class Dune::FlexibleSolver<Operator>; \
243template Dune::FlexibleSolver<Operator>::FlexibleSolver(Operator& op, \
244 const Comm& comm, \
245 const Opm::PropertyTree& prm, \
246 const std::function<typename Operator::domain_type()>& weightsCalculator, \
247 std::size_t pressureIndex);
248#define INSTANTIATE_FLEXIBLESOLVER(N) \
249INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
250INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>); \
251INSTANTIATE_FLEXIBLESOLVER_OP(ParOpM<N>); \
252INSTANTIATE_FLEXIBLESOLVER_OP(ParOpW<N>);
253
254#else // HAVE_MPI
255
256#define INSTANTIATE_FLEXIBLESOLVER_OP(Operator) \
257template class Dune::FlexibleSolver<Operator>;
258
259#define INSTANTIATE_FLEXIBLESOLVER(N) \
260INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpM<N>); \
261INSTANTIATE_FLEXIBLESOLVER_OP(SeqOpW<N>);
262
263#endif // HAVE_MPI
264
265
266
267#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
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:486
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