1#ifndef DUNE_FEM_ISTLPRECONDITIONERWRAPPER_HH
2#define DUNE_FEM_ISTLPRECONDITIONERWRAPPER_HH
10#include <dune/common/version.hh>
12#include <dune/istl/operators.hh>
13#include <dune/istl/preconditioners.hh>
15#include <dune/istl/paamg/amg.hh>
16#include <dune/istl/paamg/pinfo.hh>
28 template <
class MatrixImp>
31 template <
class MatrixImp>
34 template <
class MatrixImp>
68 template<
class MatrixObject,
class X,
class Y >
69 class FemDiagonalPreconditioner :
public Preconditioner<X,Y>
72 typedef typename MatrixObject :: ColumnDiscreteFunctionType DiscreteFunctionType ;
74 typedef DiagonalPreconditionerBase< DiscreteFunctionType, MatrixObject, true > PreconditionerType ;
77 PreconditionerType diagonalPrecon_;
80 typedef typename X::field_type field_type;
81 FemDiagonalPreconditioner(
const MatrixObject& mObj )
82 : diagonalPrecon_( mObj )
86 void pre (X& x, Y& b)
override {}
89 void apply (X& v,
const Y& d)
override
91 diagonalPrecon_.applyToISTLBlockVector( d, v );
95 void post (X& x)
override {}
97#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
99 SolverCategory::Category category ()
const override {
return SolverCategory::sequential; }
104 template<
class X,
class Y >
105 class IdentityPreconditionerWrapper :
public Preconditioner<X,Y>
107 template <
class XImp,
class YImp>
110 inline static void copy(XImp& v,
const YImp& d)
115 template <
class XImp>
116 struct Apply<XImp,XImp>
118 inline static void copy(X& v,
const Y& d)
126 typedef X domain_type;
128 typedef Y range_type;
130 typedef typename X::field_type field_type;
132#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
135 category=SolverCategory::sequential };
139 IdentityPreconditionerWrapper(){}
142 void pre (X& x, Y& b)
override {}
145 void apply (X& v,
const Y& d)
override
147 Apply< X, Y> :: copy( v, d );
151 void post (X& x)
override {}
153#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
154 SolverCategory::Category category ()
const override {
return SolverCategory::sequential; }
162 template<
class MatrixImp>
163 class PreconditionerWrapper
164 :
public Preconditioner<typename MatrixImp :: RowBlockVectorType,
165 typename MatrixImp :: ColBlockVectorType>
167 typedef MatrixImp MatrixType;
168 typedef typename MatrixType :: RowBlockVectorType X;
169 typedef typename MatrixType :: ColBlockVectorType Y;
172 typedef typename MatrixType :: BaseType ISTLMatrixType ;
174 typedef typename MatrixType :: CollectiveCommunictionType
175 CollectiveCommunictionType;
182 typedef MatrixAdapter< ISTLMatrixType, X, Y > OperatorType;
184 mutable std::shared_ptr< OperatorType > op_;
187 typedef Preconditioner<X,Y> PreconditionerInterfaceType;
188 mutable std::shared_ptr<PreconditionerInterfaceType> preconder_;
194 template <
class XImp,
class YImp>
197 inline static void apply(PreconditionerInterfaceType& preconder,
198 XImp& v,
const YImp& d)
202 inline static void copy(XImp& v,
const YImp& d)
207 template <
class XImp>
208 struct Apply<XImp,XImp>
210 inline static void apply(PreconditionerInterfaceType& preconder,
211 XImp& v,
const XImp& d)
213 preconder.apply( v ,d );
216 inline static void copy(X& v,
const Y& d)
223 typedef X domain_type;
225 typedef Y range_type;
227 typedef typename X::field_type field_type;
229#if ! DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
232 category=SolverCategory::sequential };
236 PreconditionerWrapper (
const PreconditionerWrapper& org,
bool verbose)
238 , preconder_(org.preconder_)
240 , verbose_( verbose )
245 PreconditionerWrapper(
bool verbose)
249 , verbose_( verbose )
253 template <
class PreconditionerType>
254 PreconditionerWrapper(MatrixType & matrix,
258 const PreconditionerType* p)
260 , preconder_( new PreconditionerType( matrix, iter, relax ) )
262 , verbose_( verbose )
269 template <
class PreconditionerType>
270 PreconditionerWrapper(MatrixType & matrix,
bool verbose,
271 PreconditionerType* p)
275 , verbose_( verbose )
280 template <
class PreconditionerType>
281 PreconditionerWrapper(MatrixType & matrix,
285 const PreconditionerType* p ,
286 const CollectiveCommunictionType& comm)
290 : op_( new OperatorType( matrix ) )
292 , preconder_( createAMGPreconditioner(comm, iter, relax, p) )
294 , verbose_( verbose )
299 void pre (X& x, Y& b)
override
306 preconder_->pre(x,b);
312 void apply (X& v,
const Y& d)
override
317 Apply<X,Y>::apply( *preconder_ , v, d );
321 Apply<X,Y>::copy( v, d );
326 void post (X& x)
override
337#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 6)
338 SolverCategory::Category category ()
const override
340 return (preconder_ ? preconder_->category() : SolverCategory::sequential);
345 template <
class Smoother>
346 PreconditionerInterfaceType*
347 createAMGPreconditioner(
const CollectiveCommunictionType& comm,
348 int iter, field_type relax,
const Smoother* )
350 typedef typename Dune::FieldTraits< field_type>::real_type real_type;
351 typedef typename std::conditional< std::is_convertible<field_type,real_type>::value,
352 Dune::Amg::FirstDiagonal, Dune::Amg::RowSum >::type Norm;
353 typedef Dune::Amg::CoarsenCriterion<
354 Dune::Amg::UnSymmetricCriterion<ISTLMatrixType, Norm > > Criterion;
356 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
358 SmootherArgs smootherArgs;
360 smootherArgs.iterations = iter;
361 smootherArgs.relaxationFactor = relax ;
363 int coarsenTarget=1200;
364 Criterion criterion(15,coarsenTarget);
365 criterion.setDefaultValuesIsotropic(2);
366 criterion.setAlpha(.67);
367 criterion.setBeta(1.0e-8);
368 criterion.setMaxLevel(10);
369 if( verbose_ && Parameter :: verbose() )
370 criterion.setDebugLevel( 1 );
372 criterion.setDebugLevel( 0 );
386 typedef Dune::Amg::AMG<OperatorType, X, Smoother> AMG;
387 return new AMG(*op_, criterion, smootherArgs);
397 template <
class MatrixObject >
398 class ISTLMatrixAdapterFactory ;
400 template <
class DomainSpace,
class RangeSpace,
class DomainBlock,
class RangeBlock,
401 template <
class,
class,
class,
class>
class MatrixObject >
402 class ISTLMatrixAdapterFactory< MatrixObject< DomainSpace, RangeSpace, DomainBlock, RangeBlock > >
404 typedef DomainSpace DomainSpaceType ;
405 typedef RangeSpace RangeSpaceType;
408 typedef MatrixObject< DomainSpaceType, RangeSpaceType, DomainBlock, RangeBlock > MatrixObjectType;
409 typedef typename MatrixObjectType :: MatrixType MatrixType;
410 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterInterfaceType;
413 static std::unique_ptr< MatrixAdapterInterfaceType >
414 matrixAdapter(
const MatrixObjectType& matrixObj,
415 const ISTLSolverParameter& param)
417 std::unique_ptr< MatrixAdapterInterfaceType > ptr;
418 if( matrixObj.domainSpace().continuous() )
420 typedef LagrangeParallelMatrixAdapter< MatrixType > MatrixAdapterImplementation;
421 ptr.reset( matrixAdapterObject( matrixObj, (MatrixAdapterImplementation *)
nullptr, param ) );
425 typedef DGParallelMatrixAdapter< MatrixType > MatrixAdapterImplementation;
426 ptr.reset( matrixAdapterObject( matrixObj, (MatrixAdapterImplementation *)
nullptr, param ) );
432 template <
class MatrixAdapterType>
433 static MatrixAdapterType*
434 matrixAdapterObject(
const MatrixObjectType& matrixObj,
435 const MatrixAdapterType*,
436 const ISTLSolverParameter& param )
438 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
439 return new MatrixAdapterType( matrixObj.exportMatrix(),
440 matrixObj.domainSpace(), matrixObj.rangeSpace(), PreConType(param.verbose()) );
445#ifndef DISABLE_ISTL_PRECONDITIONING
447 template <
class Space,
class DomainBlock,
class RangeBlock,
448 template <
class,
class,
class,
class>
class MatrixObject >
449 class ISTLMatrixAdapterFactory< MatrixObject< Space, Space, DomainBlock, RangeBlock > >
452 typedef Space DomainSpaceType ;
453 typedef Space RangeSpaceType;
455 typedef MatrixObject< DomainSpaceType, RangeSpaceType, DomainBlock, RangeBlock > MatrixObjectType;
456 typedef typename MatrixObjectType :: MatrixType MatrixType;
457 typedef typename MatrixType :: BaseType ISTLMatrixType ;
458 typedef Fem::PreconditionerWrapper< MatrixType > PreconditionAdapterType;
461 template <
class MatrixAdapterType,
class PreconditionerType>
462 static MatrixAdapterType*
463 createMatrixAdapter(
const MatrixAdapterType*,
464 const PreconditionerType* preconditioning,
466 const DomainSpaceType& domainSpace,
467 const RangeSpaceType& rangeSpace,
468 const double relaxFactor,
469 std::size_t numIterations,
472 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
473 PreConType preconAdapter(matrix, numIterations, relaxFactor, verbose, preconditioning );
474 return new MatrixAdapterType(matrix, domainSpace, rangeSpace, preconAdapter );
477 template <
class MatrixAdapterType,
class PreconditionerType>
478 static MatrixAdapterType*
479 createAMGMatrixAdapter(
const MatrixAdapterType*,
480 const PreconditionerType* preconditioning,
482 const DomainSpaceType& domainSpace,
483 const RangeSpaceType& rangeSpace,
484 const double relaxFactor,
485 std::size_t numIterations,
488 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
489 PreConType preconAdapter(matrix, numIterations, relaxFactor, solverVerbose,
490 preconditioning, domainSpace.gridPart().comm() );
491 return new MatrixAdapterType(matrix, domainSpace, rangeSpace, preconAdapter );
494 template <
class MatrixAdapterType>
495 static MatrixAdapterType*
496 matrixAdapterObject(
const MatrixObjectType& matrixObj,
497 const MatrixAdapterType*,
498 const ISTLSolverParameter& param )
500 int preconditioning = param.preconditionMethod(
501 { SolverParameter::none,
502 SolverParameter::ssor,
503 SolverParameter::sor ,
504 SolverParameter::ilu ,
505 SolverParameter::gauss_seidel,
506 SolverParameter::jacobi,
507 SolverParameter::amg_ilu,
508 SolverParameter::amg_jacobi,
509 SolverParameter::ildl
511 const double relaxFactor = param.relaxation();
512 const size_t numIterations = param.preconditionerIteration();
514 const DomainSpaceType& domainSpace = matrixObj.domainSpace();
515 const RangeSpaceType& rangeSpace = matrixObj.rangeSpace();
517 MatrixType& matrix = matrixObj.exportMatrix();
518 const auto procs = domainSpace.gridPart().comm().size();
520 typedef typename MatrixType :: BaseType ISTLMatrixType;
521 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
523 typedef typename MatrixObjectType :: RowBlockVectorType RowBlockVectorType;
524 typedef typename MatrixObjectType :: ColumnBlockVectorType ColumnBlockVectorType;
527 if( preconditioning == SolverParameter::none )
529 return new MatrixAdapterType( matrix, domainSpace, rangeSpace, PreConType(param.verbose()) );
532 else if( preconditioning == SolverParameter::ssor )
535 DUNE_THROW(InvalidStateException,
"ISTL::SeqSSOR not working in parallel computations");
537 typedef SeqSSOR<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
538 return createMatrixAdapter( (MatrixAdapterType *)
nullptr,
539 (PreconditionerType*)
nullptr,
540 matrix, domainSpace, rangeSpace, relaxFactor, numIterations, param.verbose() );
543 else if(preconditioning == SolverParameter::sor )
546 DUNE_THROW(InvalidStateException,
"ISTL::SeqSOR not working in parallel computations");
548 typedef SeqSOR<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
549 return createMatrixAdapter( (MatrixAdapterType *)
nullptr,
550 (PreconditionerType*)
nullptr,
551 matrix, domainSpace, rangeSpace, relaxFactor, numIterations, param.verbose() );
554 else if(preconditioning == SolverParameter::ilu)
557 DUNE_THROW(InvalidStateException,
"ISTL::SeqILU not working in parallel computations");
559 typedef SeqILU<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
560 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
562 PreConType preconAdapter( matrix, param.verbose(),
new PreconditionerType( matrix, numIterations, relaxFactor, param.fastILUStorage()) );
563 return new MatrixAdapterType( matrix, domainSpace, rangeSpace, preconAdapter );
566 else if(preconditioning == SolverParameter::gauss_seidel)
569 DUNE_THROW(InvalidStateException,
"ISTL::SeqGS not working in parallel computations");
571 typedef SeqGS<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
572 return createMatrixAdapter( (MatrixAdapterType *)
nullptr,
573 (PreconditionerType*)
nullptr,
574 matrix, domainSpace, rangeSpace, relaxFactor, numIterations, param.verbose() );
577 else if(preconditioning == SolverParameter::jacobi)
579 if( numIterations == 1 )
581 typedef FemDiagonalPreconditioner< MatrixObjectType, RowBlockVectorType, ColumnBlockVectorType > PreconditionerType;
582 typedef typename MatrixAdapterType :: PreconditionAdapterType PreConType;
583 PreConType preconAdapter( matrix, param.verbose(),
new PreconditionerType( matrixObj ) );
584 return new MatrixAdapterType( matrix, domainSpace, rangeSpace, preconAdapter );
586 else if ( procs == 1 )
588 typedef SeqJac<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
589 return createMatrixAdapter( (MatrixAdapterType *)
nullptr,
590 (PreconditionerType*)
nullptr,
591 matrix, domainSpace, rangeSpace, relaxFactor, numIterations,
596 DUNE_THROW(InvalidStateException,
"ISTL::SeqJac(Jacobi) only working with istl.preconditioning.iterations: 1 in parallel computations");
600 else if(preconditioning == SolverParameter::amg_ilu)
603 typedef SeqILU<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType> PreconditionerType;
604 return createAMGMatrixAdapter( (MatrixAdapterType *)
nullptr,
605 (PreconditionerType*)
nullptr,
606 matrix, domainSpace, rangeSpace, relaxFactor, numIterations,
610 else if(preconditioning == SolverParameter::amg_jacobi)
612 typedef SeqJac<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType,1> PreconditionerType;
613 return createAMGMatrixAdapter( (MatrixAdapterType *)
nullptr,
614 (PreconditionerType*)
nullptr,
615 matrix, domainSpace, rangeSpace, relaxFactor, numIterations,
619 else if(preconditioning == SolverParameter::ildl)
622 DUNE_THROW(InvalidStateException,
"ISTL::SeqILDL not working in parallel computations");
624 PreConType preconAdapter( matrix, param.verbose(),
new SeqILDL<ISTLMatrixType,RowBlockVectorType,ColumnBlockVectorType>( matrix , relaxFactor ) );
625 return new MatrixAdapterType( matrix, domainSpace, rangeSpace, preconAdapter );
629 preConErrorMsg(preconditioning);
631 return new MatrixAdapterType(matrix, domainSpace, rangeSpace, PreConType(param.verbose()) );
634 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterInterfaceType;
636 static void preConErrorMsg(
int preCon)
639 std::cerr <<
"ERROR: Wrong precoditioning number (p = " << preCon;
640 std::cerr <<
") in ISTLMatrixObject! \n";
641 std::cerr <<
"Valid values are: \n";
642 std::cerr <<
"0 == no \n";
643 std::cerr <<
"1 == SSOR \n";
644 std::cerr <<
"2 == SOR \n";
645 std::cerr <<
"3 == ILU-0 \n";
646 std::cerr <<
"4 == ILU-n \n";
647 std::cerr <<
"5 == Gauss-Seidel \n";
648 std::cerr <<
"6 == Jacobi \n";
651 DUNE_THROW(NotImplemented,
"Wrong precoditioning selected");
656 static std::unique_ptr< MatrixAdapterInterfaceType >
657 matrixAdapter(
const MatrixObjectType& matrixObj,
658 const ISTLSolverParameter& param)
660 const ISTLSolverParameter* parameter =
dynamic_cast< const ISTLSolverParameter*
> (¶m);
661 std::unique_ptr< ISTLSolverParameter > paramPtr;
664 paramPtr.reset(
new ISTLSolverParameter( param ) );
665 parameter = paramPtr.operator->();
668 std::unique_ptr< MatrixAdapterInterfaceType > ptr;
669 if( matrixObj.domainSpace().continuous() )
671 typedef LagrangeParallelMatrixAdapter< MatrixType > MatrixAdapterImplementation;
672 ptr.reset( matrixAdapterObject( matrixObj, (MatrixAdapterImplementation *)
nullptr, *parameter ) );
676 typedef DGParallelMatrixAdapter< MatrixType > MatrixAdapterImplementation;
677 ptr.reset( matrixAdapterObject( matrixObj, (MatrixAdapterImplementation *)
nullptr, *parameter ) );
Definition: bindguard.hh:11
static ParameterContainer & container()
Definition: io/parameter.hh:193
Definition: io/parameter.hh:551
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:159
Definition: istlpreconditioner.hh:29
Definition: istlpreconditioner.hh:35
Definition: istlpreconditioner.hh:38
virtual bool fastILUStorage() const
Definition: istlpreconditioner.hh:61
LocalParameter< SolverParameter, ISTLSolverParameter > BaseType
Definition: istlpreconditioner.hh:39
ISTLSolverParameter(const std::string &keyPrefix, const ParameterReader ¶meter=Parameter::container())
Definition: istlpreconditioner.hh:48
virtual double overflowFraction() const
Definition: istlpreconditioner.hh:56
ISTLSolverParameter(const ParameterReader ¶meter=Parameter::container())
Definition: istlpreconditioner.hh:44
ISTLSolverParameter(const SolverParameter &sp)
Definition: istlpreconditioner.hh:52
Definition: solver/parameter.hh:15
const ParameterReader & parameter() const
Definition: solver/parameter.hh:68
const std::string & keyPrefix() const
Definition: solver/parameter.hh:66