dune-fem 2.8.0
Loading...
Searching...
No Matches
istlinverseoperators.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
2#define DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
3
9
12
13#if HAVE_DUNE_ISTL
14#include <dune/common/version.hh>
15
18
19#include <dune/istl/scalarproducts.hh>
20#include <dune/istl/solvers.hh>
21#include <dune/istl/preconditioner.hh>
22
23namespace Dune
24{
25
26 namespace Fem
27 {
28
29 // wrapper for Fem Preconditioners (Operators acting as preconditioners) into ISTL preconditioners
30 template< class Preconditioner >
31 class ISTLPreconditionAdapter
32 : public Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType >
33 {
34 typedef ISTLPreconditionAdapter< Preconditioner > ThisType;
35 typedef Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType > BaseType;
36
37 typedef typename Preconditioner::DomainFunctionType DomainFunctionType;
38 typedef typename Preconditioner::RangeFunctionType RangeFunctionType;
39
40 public:
41#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
42 enum {category=SolverCategory::sequential};
43#endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
44
45 typedef typename BaseType::domain_type domain_type;
46 typedef typename BaseType::range_type range_type;
47 typedef typename BaseType::field_type field_type;
48
49 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
50 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
51
52 ISTLPreconditionAdapter ( const Preconditioner *precon, const DomainFunctionSpaceType &domainSpace, const RangeFunctionSpaceType &rangeSpace )
53 : precon_( precon ),
54 domainSpace_( domainSpace ),
55 rangeSpace_( rangeSpace )
56 {}
57
58 // pre and post do nothing here
59 virtual void pre ( domain_type &x, range_type &y ) override {}
60 virtual void post ( domain_type &x ) override {}
61
62 virtual void apply ( domain_type &x, const range_type &y ) override
63 {
64 // no precon
65 if( !precon_ )
66 {
67 x = y;
68 }
69 else
70 {
71 // note: ISTL switches the arguments !!!
72 // it is assumed that we have a left preconditioner
73 RangeFunctionType px( "ISTLPreconditionAdapter::apply::x", rangeSpace_, x );
74 DomainFunctionType py( "ISTLPreconditionAdapter::apply::y", domainSpace_, y );
75
76 (*precon_)( px, py );
77 }
78 }
79
80#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
81 SolverCategory::Category category () const override { return SolverCategory::sequential; }
82#endif // #if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
83
84 protected:
85 const Preconditioner *precon_;
86 const DomainFunctionSpaceType &domainSpace_;
87 const RangeFunctionSpaceType &rangeSpace_;
88 };
89
90
91 template< class BlockVector >
92 struct ISTLSolverReduction
93 {
94 ISTLSolverReduction ( std::shared_ptr< ISTLSolverParameter > parameter )
95 : parameter_( parameter )
96 {}
97
98 double operator() ( const Dune::LinearOperator< BlockVector, BlockVector > &op,
99 Dune::ScalarProduct< BlockVector > &scp,
100 const BlockVector &rhs, const BlockVector &x ) const
101 {
102
103 if( parameter_->errorMeasure() == 0)
104 {
105 BlockVector residuum( rhs );
106 op.applyscaleadd( -1., x, residuum );
107 const double res = scp.norm( residuum );
108 return (res > 0 ? parameter_->tolerance() / res : 1e-3);
109 }
110 else
111 return parameter_->tolerance();
112 }
113
114 private:
115 std::shared_ptr<ISTLSolverParameter> parameter_;
116 };
117
118 template< int method,
119 class X,
120 class Reduction = ISTLSolverReduction< X > >
121 struct ISTLSolverAdapter
122 {
123 typedef Reduction ReductionType;
124
125 typedef X domain_type;
126 typedef X range_type;
127
128 ISTLSolverAdapter ( const ReductionType &reduction, std::shared_ptr<ISTLSolverParameter> parameter )
129 : reduction_( reduction ),
130 method_( method < 0 ?
131 parameter->solverMethod({ SolverParameter::gmres,
132 SolverParameter::cg,
133 SolverParameter::bicgstab,
134 SolverParameter::minres,
135 SolverParameter::gradient,
136 SolverParameter::loop,
137 SolverParameter::superlu
138 })
139 : method ),
140 parameter_( parameter )
141 {}
142
143 template<class Op, class ScP, class PC >
144 void operator () ( Op& op, ScP &scp, PC &pc,
145 range_type &rhs, domain_type &x,
146 Dune::InverseOperatorResult &result ) const
147 {
148 const int verbosity = (Dune::Fem::Parameter::verbose() && parameter_->verbose()) ? 2 : 0;
149 int maxIterations = std::min( std::numeric_limits< int >::max(), parameter_->maxIterations() );
150 if( method_ == SolverParameter::cg )
151 {
152 typedef Dune::CGSolver< X > SolverType;
153 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
154 solver.apply( x, rhs, result );
155 return ;
156 }
157 else if( method_ == SolverParameter::bicgstab )
158 {
159 typedef Dune::BiCGSTABSolver< X > SolverType;
160 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
161 solver.apply( x, rhs, result );
162 return ;
163 }
164 else if( method_ == SolverParameter::gmres )
165 {
166 typedef Dune::RestartedGMResSolver< X > SolverType;
167 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), parameter_->gmresRestart(), maxIterations, verbosity );
168 solver.apply( x, rhs, result );
169 return ;
170 }
171 else if( method_ == SolverParameter::minres )
172 {
173 typedef Dune::MINRESSolver< X > SolverType;
174 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
175 solver.apply( x, rhs, result );
176 }
177 else if( method_ == SolverParameter::gradient )
178 {
179 typedef Dune::GradientSolver< X > SolverType;
180 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
181 solver.apply( x, rhs, result );
182 return ;
183 }
184 else if( method_ == SolverParameter::loop )
185 {
186 typedef Dune::LoopSolver< X > SolverType;
187 SolverType solver( op, scp, pc, reduction_( op, scp, rhs, x ), maxIterations, verbosity );
188 solver.apply( x, rhs, result );
189 return ;
190 }
191 else if( method_ == SolverParameter::superlu )
192 {
193 callSuperLU( op, rhs, x, result );
194 }
195 else
196 {
197 DUNE_THROW(NotImplemented,"ISTLSolverAdapter::operator(): wrong method solver identifier" << method_ );
198 }
199 }
200
201 void setMaxLinearIterations( unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
202 void setMaxIterations( unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
203
204 std::shared_ptr<ISTLSolverParameter> parameter () const { return parameter_; }
205
206 protected:
207 template< class ImprovedMatrix >
208 void callSuperLU ( ISTLParallelMatrixAdapterInterface< ImprovedMatrix >& op,
209 range_type &rhs, domain_type &x,
210 Dune::InverseOperatorResult &result ) const
211 {
212#if HAVE_SUPERLU
213 const int verbosity = (Dune::Fem::Parameter::verbose() && parameter_->verbose()) ? 2 : 0;
214 typedef typename ImprovedMatrix :: BaseType Matrix;
215 const ImprovedMatrix& matrix = op.getmat();
216 SuperLU< Matrix > solver( matrix, verbosity );
217 solver.apply( x, rhs, result );
218#else
219 DUNE_THROW(NotImplemented,"ISTLSolverAdapter::callSuperLU: SuperLU solver selected but SuperLU not available!");
220#endif
221 }
222
223 template< class Op >
224 void callSuperLU ( ISTLLinearOperatorAdapter< Op >& op, range_type &rhs, domain_type &x,
225 Dune::InverseOperatorResult &result ) const
226 {
227 DUNE_THROW(NotImplemented,"ISTLSolverAdapter::callSuperLU: SuperLU only works for AssembledLinearOperators!");
228 }
229
230 ReductionType reduction_;
231 const int method_;
232
233 std::shared_ptr<ISTLSolverParameter> parameter_;
234 };
235
236
237
238 // ISTLInverseOperator
239 // -------------------
240
241 template< class DiscreteFunction, int method = -1,
242 class Preconditioner = const Operator< DiscreteFunction, DiscreteFunction > >
243 class ISTLInverseOperator;
244
245 template< class DiscreteFunction, int method, class Preconditioner >
246 struct ISTLInverseOperatorTraits
247 {
248 typedef DiscreteFunction DiscreteFunctionType;
250 typedef Preconditioner PreconditionerType;
251
252 typedef Dune::Fem::ISTLLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType;
253 typedef ISTLBlockVectorDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType ;
254
255 typedef ISTLInverseOperator< DiscreteFunction, method, Preconditioner > InverseOperatorType;
256 typedef ISTLSolverParameter SolverParameterType;
257 };
258
259 template< class DiscreteFunction, int method, class Preconditioner >
260 class ISTLInverseOperator : public InverseOperatorInterface<
261 ISTLInverseOperatorTraits< DiscreteFunction,
262 method, Preconditioner > >
263 {
264 typedef ISTLInverseOperatorTraits< DiscreteFunction, method, Preconditioner > Traits;
265 typedef InverseOperatorInterface< Traits > BaseType;
266
267 friend class InverseOperatorInterface< Traits >;
268 public:
269 typedef typename BaseType::DomainFunctionType DomainFunctionType;
270 typedef typename BaseType::RangeFunctionType RangeFunctionType;
271 typedef typename BaseType::OperatorType OperatorType;
272 typedef typename BaseType::PreconditionerType PreconditionerType;
273 typedef typename BaseType::AssembledOperatorType AssembledOperatorType;
274
275 using BaseType :: bind;
276 using BaseType :: unbind;
277 using BaseType :: setMaxIterations;
278 using BaseType :: setMaxLinearIterations;
279
280 protected:
281 typedef typename DomainFunctionType :: DiscreteFunctionSpaceType
282 DiscreteFunctionSpaceType;
283
284 typedef ISTLLinearOperatorAdapter< OperatorType > ISTLOperatorType;
285 typedef ISTLPreconditionAdapter< PreconditionerType > ISTLPreconditionerAdapterType;
286
287 typedef typename RangeFunctionType :: ScalarProductType ParallelScalarProductType;
288 typedef typename RangeFunctionType::DofStorageType BlockVectorType;
289
290 typedef ISTLSolverAdapter< method, BlockVectorType > SolverAdapterType;
291 typedef typename SolverAdapterType::ReductionType ReductionType;
292 public:
293
294 //non-deprecated constructors
295 ISTLInverseOperator ( const ISTLSolverParameter & parameter = ISTLSolverParameter() )
296 : BaseType( parameter ), solverAdapter_( ReductionType( parameter_ ), parameter_ )
297 {}
298
299 ISTLInverseOperator ( const OperatorType &op,
300 const ISTLSolverParameter & parameter = ISTLSolverParameter() )
301 : ISTLInverseOperator ( parameter )
302 {
303 bind( op );
304 }
305
306 ISTLInverseOperator ( const OperatorType &op, PreconditionerType &preconditioner,
307 const ISTLSolverParameter & parameter = ISTLSolverParameter() )
308 : ISTLInverseOperator( parameter )
309 {
310 bind( op, preconditioner );
311 }
312
313
314 protected:
315 // apply for arbitrary domain function type and matching range function type
316 template <class DomainFunction>
317 int apply( const DomainFunction& u, RangeFunctionType& w ) const
318 {
319 auto& scp = w.scalarProduct();
320 // u may not be a discrete function, therefore use w.space()
321 const DiscreteFunctionSpaceType& space = w.space();
322 ISTLPreconditionerAdapterType istlPreconditioner( preconditioner_, space, space );
323
324 if( assembledOperator_ )
325 {
326 auto& matrix = assembledOperator_->matrixAdapter( *(solverAdapter_.parameter()) );
327 // if preconditioner_ was set use that one, otherwise the one from the matrix object
328 typedef Dune::Preconditioner< BlockVectorType, BlockVectorType > PreconditionerType;
329 PreconditionerType& matrixPre = matrix.preconditionAdapter();
330 PreconditionerType& istlPre = istlPreconditioner;
331 PreconditionerType& precon = ( preconditioner_ ) ? istlPre : matrixPre;
332 return solve( matrix, scp, precon, u, w );
333 }
334 else
335 {
336 assert( operator_ );
337 ISTLOperatorType istlOperator( *operator_, space, space );
338 return solve( istlOperator, scp, istlPreconditioner, u, w );
339 }
340 }
341
343 template< class OperatorAdapter, class ISTLPreconditioner, class DomainFunction >
344 int solve ( OperatorAdapter &istlOperator, ParallelScalarProductType &scp,
345 ISTLPreconditioner &preconditioner,
346 const DomainFunction& u,
347 RangeFunctionType& w ) const
348 {
349 if( ! rhs_ )
350 {
351 // u may not be a discrete function, therefore use w.space()
352 rhs_.reset( new DomainFunctionType( "ISTLInvOp::rhs", w.space() ) );
353 rightHandSideCopied_ = false;
354 }
355
356 if( ! rightHandSideCopied_ )
357 {
358 // copy right hand side since ISTL solvers seem to modify it
359 rhs_->assign( u );
360 rightHandSideCopied_ = true;
361 }
362
363 Dune::InverseOperatorResult result;
364 solverAdapter_.setMaxIterations( parameter_->maxIterations() );
365 solverAdapter_( istlOperator, scp, preconditioner, rhs_->blockVector(), w.blockVector(), result );
366 return (result.converged) ? result.iterations : -(result.iterations);
367 }
368
369 using BaseType :: operator_;
370 using BaseType :: assembledOperator_;
371 using BaseType :: preconditioner_;
372
373 using BaseType :: rhs_;
374 using BaseType :: x_;
375
376 using BaseType :: rightHandSideCopied_;
377 using BaseType :: parameter_;
378
379 mutable SolverAdapterType solverAdapter_;
380 };
381
382 } // namespace Fem
383
384} // namespace Dune
385
386#endif // #if HAVE_ISTL
387
388#endif // #ifndef DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
double min(const Dune::Fem::Double &v, const double p)
Definition: double.hh:953
Definition: bindguard.hh:11
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445