1#ifndef DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
2#define DUNE_FEM_SOLVER_ISTLINVERSEOPERATORS_HH
14#include <dune/common/version.hh>
19#include <dune/istl/scalarproducts.hh>
20#include <dune/istl/solvers.hh>
21#include <dune/istl/preconditioner.hh>
30 template<
class Preconditioner >
31 class ISTLPreconditionAdapter
32 :
public Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType >
34 typedef ISTLPreconditionAdapter< Preconditioner > ThisType;
35 typedef Dune::Preconditioner< typename Preconditioner::RangeFunctionType::DofStorageType, typename Preconditioner::DomainFunctionType::DofStorageType > BaseType;
37 typedef typename Preconditioner::DomainFunctionType DomainFunctionType;
38 typedef typename Preconditioner::RangeFunctionType RangeFunctionType;
41#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
42 enum {category=SolverCategory::sequential};
45 typedef typename BaseType::domain_type domain_type;
46 typedef typename BaseType::range_type range_type;
47 typedef typename BaseType::field_type field_type;
49 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainFunctionSpaceType;
50 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeFunctionSpaceType;
52 ISTLPreconditionAdapter (
const Preconditioner *precon,
const DomainFunctionSpaceType &domainSpace,
const RangeFunctionSpaceType &rangeSpace )
54 domainSpace_( domainSpace ),
55 rangeSpace_( rangeSpace )
59 virtual void pre ( domain_type &x, range_type &y )
override {}
60 virtual void post ( domain_type &x )
override {}
62 virtual void apply ( domain_type &x,
const range_type &y )
override
73 RangeFunctionType px(
"ISTLPreconditionAdapter::apply::x", rangeSpace_, x );
74 DomainFunctionType py(
"ISTLPreconditionAdapter::apply::y", domainSpace_, y );
80#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
81 SolverCategory::Category category ()
const override {
return SolverCategory::sequential; }
85 const Preconditioner *precon_;
86 const DomainFunctionSpaceType &domainSpace_;
87 const RangeFunctionSpaceType &rangeSpace_;
91 template<
class BlockVector >
92 struct ISTLSolverReduction
94 ISTLSolverReduction ( std::shared_ptr< ISTLSolverParameter > parameter )
95 : parameter_( parameter )
98 double operator() (
const Dune::LinearOperator< BlockVector, BlockVector > &op,
99 Dune::ScalarProduct< BlockVector > &scp,
100 const BlockVector &rhs,
const BlockVector &x )
const
103 if( parameter_->errorMeasure() == 0)
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);
111 return parameter_->tolerance();
115 std::shared_ptr<ISTLSolverParameter> parameter_;
118 template<
int method,
120 class Reduction = ISTLSolverReduction< X > >
121 struct ISTLSolverAdapter
123 typedef Reduction ReductionType;
125 typedef X domain_type;
126 typedef X range_type;
128 ISTLSolverAdapter (
const ReductionType &reduction, std::shared_ptr<ISTLSolverParameter> parameter )
129 : reduction_( reduction ),
130 method_( method < 0 ?
131 parameter->solverMethod({ SolverParameter::gmres,
133 SolverParameter::bicgstab,
134 SolverParameter::minres,
135 SolverParameter::gradient,
136 SolverParameter::loop,
137 SolverParameter::superlu
140 parameter_( parameter )
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
149 int maxIterations =
std::min( std::numeric_limits< int >::max(), parameter_->maxIterations() );
150 if( method_ == SolverParameter::cg )
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 );
157 else if( method_ == SolverParameter::bicgstab )
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 );
164 else if( method_ == SolverParameter::gmres )
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 );
171 else if( method_ == SolverParameter::minres )
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 );
177 else if( method_ == SolverParameter::gradient )
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 );
184 else if( method_ == SolverParameter::loop )
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 );
191 else if( method_ == SolverParameter::superlu )
193 callSuperLU( op, rhs, x, result );
197 DUNE_THROW(NotImplemented,
"ISTLSolverAdapter::operator(): wrong method solver identifier" << method_ );
201 void setMaxLinearIterations(
unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
202 void setMaxIterations(
unsigned int maxIterations ) { parameter_->setMaxIterations(maxIterations); }
204 std::shared_ptr<ISTLSolverParameter> parameter ()
const {
return parameter_; }
207 template<
class ImprovedMatrix >
208 void callSuperLU ( ISTLParallelMatrixAdapterInterface< ImprovedMatrix >& op,
209 range_type &rhs, domain_type &x,
210 Dune::InverseOperatorResult &result )
const
214 typedef typename ImprovedMatrix :: BaseType Matrix;
215 const ImprovedMatrix& matrix = op.getmat();
216 SuperLU< Matrix > solver( matrix, verbosity );
217 solver.apply( x, rhs, result );
219 DUNE_THROW(NotImplemented,
"ISTLSolverAdapter::callSuperLU: SuperLU solver selected but SuperLU not available!");
224 void callSuperLU ( ISTLLinearOperatorAdapter< Op >& op, range_type &rhs, domain_type &x,
225 Dune::InverseOperatorResult &result )
const
227 DUNE_THROW(NotImplemented,
"ISTLSolverAdapter::callSuperLU: SuperLU only works for AssembledLinearOperators!");
230 ReductionType reduction_;
233 std::shared_ptr<ISTLSolverParameter> parameter_;
241 template<
class DiscreteFunction,
int method = -1,
243 class ISTLInverseOperator;
245 template<
class DiscreteFunction,
int method,
class Preconditioner >
246 struct ISTLInverseOperatorTraits
248 typedef DiscreteFunction DiscreteFunctionType;
250 typedef Preconditioner PreconditionerType;
252 typedef Dune::Fem::ISTLLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType;
253 typedef ISTLBlockVectorDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType ;
255 typedef ISTLInverseOperator< DiscreteFunction, method, Preconditioner > InverseOperatorType;
256 typedef ISTLSolverParameter SolverParameterType;
259 template<
class DiscreteFunction,
int method,
class Preconditioner >
260 class ISTLInverseOperator :
public InverseOperatorInterface<
261 ISTLInverseOperatorTraits< DiscreteFunction,
262 method, Preconditioner > >
264 typedef ISTLInverseOperatorTraits< DiscreteFunction, method, Preconditioner > Traits;
265 typedef InverseOperatorInterface< Traits > BaseType;
267 friend class InverseOperatorInterface< Traits >;
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;
275 using BaseType :: bind;
276 using BaseType :: unbind;
277 using BaseType :: setMaxIterations;
278 using BaseType :: setMaxLinearIterations;
281 typedef typename DomainFunctionType :: DiscreteFunctionSpaceType
282 DiscreteFunctionSpaceType;
284 typedef ISTLLinearOperatorAdapter< OperatorType > ISTLOperatorType;
285 typedef ISTLPreconditionAdapter< PreconditionerType > ISTLPreconditionerAdapterType;
287 typedef typename RangeFunctionType :: ScalarProductType ParallelScalarProductType;
288 typedef typename RangeFunctionType::DofStorageType BlockVectorType;
290 typedef ISTLSolverAdapter< method, BlockVectorType > SolverAdapterType;
291 typedef typename SolverAdapterType::ReductionType ReductionType;
295 ISTLInverseOperator (
const ISTLSolverParameter & parameter = ISTLSolverParameter() )
296 : BaseType( parameter ), solverAdapter_( ReductionType( parameter_ ), parameter_ )
299 ISTLInverseOperator (
const OperatorType &op,
300 const ISTLSolverParameter & parameter = ISTLSolverParameter() )
301 : ISTLInverseOperator ( parameter )
306 ISTLInverseOperator (
const OperatorType &op, PreconditionerType &preconditioner,
307 const ISTLSolverParameter & parameter = ISTLSolverParameter() )
308 : ISTLInverseOperator( parameter )
310 bind( op, preconditioner );
316 template <
class DomainFunction>
317 int apply(
const DomainFunction& u, RangeFunctionType& w )
const
319 auto& scp = w.scalarProduct();
321 const DiscreteFunctionSpaceType& space = w.space();
322 ISTLPreconditionerAdapterType istlPreconditioner( preconditioner_, space, space );
324 if( assembledOperator_ )
326 auto& matrix = assembledOperator_->matrixAdapter( *(solverAdapter_.parameter()) );
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 );
337 ISTLOperatorType istlOperator( *operator_, space, space );
338 return solve( istlOperator, scp, istlPreconditioner, u, w );
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
352 rhs_.reset(
new DomainFunctionType(
"ISTLInvOp::rhs", w.space() ) );
353 rightHandSideCopied_ =
false;
356 if( ! rightHandSideCopied_ )
360 rightHandSideCopied_ =
true;
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);
369 using BaseType :: operator_;
370 using BaseType :: assembledOperator_;
371 using BaseType :: preconditioner_;
373 using BaseType :: rhs_;
374 using BaseType :: x_;
376 using BaseType :: rightHandSideCopied_;
377 using BaseType :: parameter_;
379 mutable SolverAdapterType solverAdapter_;
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