dune-fem  2.8-git
amgistl.hh
Go to the documentation of this file.
1 #ifndef DUNE_ASH_SOLVER_ISTL_HH
2 #define DUNE_ASH_SOLVER_ISTL_HH
3 
4 #include <cassert>
5 
6 #include <limits>
7 #include <map>
8 #include <memory>
9 #include <type_traits>
10 #include <tuple>
11 #include <utility>
12 #include <vector>
13 
14 #include <dune/common/exceptions.hh>
15 #include <dune/common/ftraits.hh>
16 #include <dune/common/hybridutilities.hh>
17 #include <dune/common/parallel/remoteindices.hh>
18 
19 #include <dune/geometry/dimension.hh>
20 
21 #include <dune/grid/common/datahandleif.hh>
22 #include <dune/grid/common/rangegenerators.hh>
23 #include <dune/grid/common/partitionset.hh>
24 
25 #if HAVE_DUNE_ISTL
26 #include <dune/istl/bcrsmatrix.hh>
27 #include <dune/istl/operators.hh>
28 #include <dune/istl/paamg/amg.hh>
29 #include <dune/istl/preconditioners.hh>
30 #include <dune/istl/schwarz.hh>
31 #include <dune/istl/solvers.hh>
32 
33 #include <dune/fem/io/parameter.hh>
37 
38 
39 namespace Dune
40 {
41 
42  namespace Amg
43  {
44 
45  template< class M, class X, class Y >
46  struct ConstructionTraits< SeqILDL< M, X, Y > >
47  {
48  typedef DefaultConstructionArgs< SeqILDL< M, X, Y > > Arguments;
49 
50  static SeqILDL< M, X, Y > *construct ( Arguments &args )
51  {
52  return new SeqILDL< M, X, Y >( args.getMatrix(), args.getArgs().relaxationFactor );
53  }
54 
55  static void deconstruct ( SeqILDL< M, X, Y > *p ) { delete p; }
56  };
57 
58  } // namespace Amg
59 
60 
61 
62  namespace Fem
63  {
64 
65  // External Forward Declarations
66  // -----------------------------
67 
68  template< class DiscreteFunctionSpace >
69  class HierarchicalDiscreteFunction;
70 
71  template< class DomainFunction, class RangeFunction >
72  class ISTLLinearOperator;
73 
74  } // namespace Fem
75 
76 
77 
78  namespace Fem
79  {
80 
81  namespace ISTL
82  {
83 
84  // VectorType
85  // ----------
86 
87  template< class DiscreteFunction >
88  using VectorType = std::decay_t< decltype( std::declval< const DiscreteFunction & >().dofVector().array() ) >;
89 
90 
91 
92  // MatrixType
93  // ----------
94 
95  template< class LinearOperator >
96  struct __MatrixType
97  {
98  typedef std::decay_t< decltype( std::declval< const LinearOperator & >().matrix() ) > Type;
99  };
100 
101  template< class DomainFunction, class RangeFunction >
102  struct __MatrixType< Dune::Fem::ISTLLinearOperator< DomainFunction, RangeFunction > >
103  {
104  typedef Dune::BCRSMatrix< typename Dune::Fem::ISTLLinearOperator< DomainFunction, RangeFunction >::LittleBlockType > Type;
105  };
106 
107 
108  template< class LinearOperator >
109  using MatrixType = typename __MatrixType< LinearOperator >::Type;
110 
111 
112 
113  // Symmetry
114  // --------
115 
116  enum Symmetry : bool { unsymmetric = false, symmetric = true };
117 
118 
119 
120  // FemCommunication
121  // ----------------
122 
123  template< class DiscreteFunctionSpace >
124  class FemCommunication
125  {
126  typedef FemCommunication< DiscreteFunctionSpace > ThisType;
127 
128  public:
130 
131  typedef int GlobalLookupIndexSet;
132 
133  explicit FemCommunication ( const DiscreteFunctionSpaceType &dfSpace, Dune::SolverCategory::Category solverCategory = Dune::SolverCategory::sequential )
134  : dfSpace_( dfSpace ), solverCategory_( solverCategory )
135  {}
136 
137  const typename DiscreteFunctionSpace::GridPartType::CollectiveCommunicationType &communicator () const { return dfSpace_.gridPart().comm(); }
138 
139  template< class T >
140  void copyOwnerToAll ( const T &x, T &y ) const
141  {
142  y = x;
144  DiscreteFunctionType z( "", dfSpace_, reinterpret_cast< typename DiscreteFunctionType::DofVectorType & >( y ) );
145  z.communicate();
146  }
147 
148  template< class T >
149  void project ( T &x ) const
150  {
151  typedef typename T::field_type field_type;
152 
153  // clear auxiliary DoFs
154  const auto &auxiliaryDofs = dfSpace_.auxiliaryDofs();
155  for( int i : auxiliaryDofs )
156  x[ i ] = field_type( 0 );
157  }
158 
159  template< class T, class F >
160  void dot ( const T &x, const T &y, F &scp ) const
161  {
162  const auto &auxiliaryDofs = dfSpace_.auxiliaryDofs();
163 
164  const int numAuxiliarys = auxiliaryDofs.size();
165  for( int auxiliary = 0, i = 0; auxiliary < numAuxiliarys; ++auxiliary, ++i )
166  {
167  const int nextAuxiliary = auxiliaryDofs[ auxiliary ];
168  for( ; i < nextAuxiliary; ++i )
169  scp += x[ i ] * y[ i ];
170  }
171 
172  scp = communicator().sum( scp );
173  }
174 
175  template< class T >
176  typename Dune::FieldTraits< typename T::field_type >::real_type norm ( const T &x ) const
177  {
178  using std::sqrt;
179  typename Dune::FieldTraits< typename T::field_type >::real_type norm2( 0 );
180  dot( x, x, norm2 );
181  return sqrt( norm2 );
182  }
183 
184  Dune::SolverCategory::Category getSolverCategory () const { return solverCategory_; }
185 
186  private:
187  const DiscreteFunctionSpaceType &dfSpace_;
188  Dune::SolverCategory::Category solverCategory_;
189  };
190 
191 
192 
193  // BuildRemoteIndicesDataHandle
194  // ----------------------------
195 
196  template< class Mapper, class GlobalLookup >
197  struct BuildRemoteIndicesDataHandle
198  : public Dune::CommDataHandleIF< BuildRemoteIndicesDataHandle< Mapper, GlobalLookup >, int >
199  {
200  typedef typename GlobalLookup::GlobalIndex GlobalIndexType;
201  typedef typename GlobalLookup::LocalIndex::Attribute AttributeType;
202 
203  BuildRemoteIndicesDataHandle ( int rank, const Mapper &mapper, const GlobalLookup &globalLookup )
204  : rank_( rank ), mapper_( mapper ), globalLookup_( globalLookup )
205  {}
206 
207  bool contains ( int dim, int codim ) const { return mapper_.contains( codim ); }
208  bool fixedSize( int dim, int codim ) const { return true; }
209 
210  template< class Buffer, class Entity >
211  void gather ( Buffer &buffer, const Entity &entity ) const
212  {
213  buffer.write( rank_ );
214  int attribute = -1;
215  mapper_.mapEachEntityDof( entity, [ this, &attribute ] ( int, auto index ) {
216  auto *pair = globalLookup_.pair( index );
217  assert( pair && ((attribute == -1) || (attribute == pair->local().attribute())) );
218  attribute = pair->local().attribute();
219  } );
220  buffer.write( static_cast< int >( attribute ) );
221  }
222 
223  template< class Buffer, class Entity >
224  void scatter ( Buffer &buffer, const Entity &entity, std::size_t n )
225  {
226  int rank, attribute;
227  buffer.read( rank );
228  buffer.read( attribute );
229  assert( (attribute != -1) || (mapper_.numEntityDofs( entity ) == 0) );
230  mapper_.mapEachEntityDof( entity, [ this, rank, attribute ] ( int, auto index ) {
231  auto *pair = globalLookup_.pair( index );
232  assert( pair );
233  remotes[ rank ].emplace_back( static_cast< AttributeType >( attribute ), pair );
234  } );
235  }
236 
237  template< class Entity >
238  std::size_t size ( const Entity &entity ) const
239  {
240  return 2;
241  }
242 
243  std::map< int, std::vector< Dune::RemoteIndex< GlobalIndexType, AttributeType > > > remotes;
244 
245  private:
246  int rank_;
247  const Mapper &mapper_;
248  const GlobalLookup &globalLookup_;
249  };
250 
251 
252 
253  // buildCommunication
254  // ------------------
255 
256  template< class DiscreteFunction >
257  void buildCommunication ( const typename DiscreteFunction::DiscreteFunctionSpaceType &dfSpace,
258  Dune::SolverCategory::Category solverCategory,
259  std::shared_ptr< FemCommunication< DiscreteFunction > > &communication )
260  {
261  communication.reset( new FemCommunication< DiscreteFunction >( dfSpace, solverCategory ) );
262  }
263 
264  template< class DiscreteFunctionSpace, class GlobalId, class LocalId >
265  void buildCommunication ( const DiscreteFunctionSpace &dfSpace,
266  Dune::SolverCategory::Category solverCategory,
267  std::shared_ptr< Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId > > &communication )
268  {
269  typedef typename DiscreteFunctionSpace::GridPartType GridPartType;
270  typedef typename DiscreteFunctionSpace::BlockMapperType LocalMapperType;
271 
272  typedef typename Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId >::GlobalLookupIndexSet GlobalLookupType;
273 
274  typedef typename GlobalLookupType::LocalIndex LocalIndexType;
275 
276  communication.reset( new Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId >( solverCategory ) );
277 
278  const GridPartType &gridPart = dfSpace.gridPart();
279  LocalMapperType &localMapper = dfSpace.blockMapper();
280 
281  // create global index mapping
283 
284  // construct local attributes
285  std::vector< typename LocalIndexType::Attribute > attribute( localMapper.size(), Dune::OwnerOverlapCopyAttributeSet::owner );
286  for( const auto &auxiliary : dfSpace.auxiliaryDofs() )
287  attribute[ auxiliary ] = Dune::OwnerOverlapCopyAttributeSet::copy;
288 
289  // build parallel index set
290  communication->indexSet().beginResize();
291  for( LocalId i = 0, size = localMapper.size(); i < size; ++i )
292  communication->indexSet().add( globalMapper.mapping()[ i ], LocalIndexType( i, attribute[ i ] ) );
293  communication->indexSet().endResize();
294 
295  // build remote indices
296  communication->buildGlobalLookup();
297  BuildRemoteIndicesDataHandle< LocalMapperType, GlobalLookupType > buildRemoteIndicesDataHandle( gridPart.comm().rank(), localMapper, communication->globalLookup() );
298  gridPart.communicate( buildRemoteIndicesDataHandle, Dune::All_All_Interface, Dune::ForwardCommunication );
299  communication->freeGlobalLookup();
300 
301  communication->remoteIndices().setIndexSets( communication->indexSet(), communication->indexSet(), communication->communicator() );
302  if( !buildRemoteIndicesDataHandle.remotes.empty() )
303  {
304  for( auto &remote : buildRemoteIndicesDataHandle.remotes )
305  {
306  std::sort( remote.second.begin(), remote.second.end(), [] ( const auto &a, const auto &b ) { return (a.localIndexPair().global() < b.localIndexPair().global()); } );
307  auto modifier = communication->remoteIndices().template getModifier< false, true >( remote.first );
308  for( const auto &remoteIndex : remote.second )
309  modifier.insert( remoteIndex );
310  }
311  }
312  else
313  communication->remoteIndices().template getModifier< false, true >( 0 );
314  }
315 
316 
317 
318  // NamedType
319  // ---------
320 
321  template< class T >
322  struct NamedType
323  {
324  typedef T Type;
325 
326  NamedType ( std::string name ) : name( name ) {}
327 
328  std::string name;
329  };
330 
331 
332 
333  // getEnum
334  // -------
335 
336  template< class... T, class F >
337  void getEnum ( const Dune::Fem::ParameterReader &parameter, const std::string &key, const std::tuple< NamedType< T >... > &types, const std::string &defaultValue, F &&f )
338  {
339  const std::string &value = parameter.getValue( key, defaultValue );
340  bool success = false;
341  Dune::Hybrid::forEach( std::index_sequence_for< T... >(), [ &types, &f, &value, &success ] ( auto &&i ) -> void {
342  const auto &type = std::get< std::decay_t< decltype( i ) >::value >( types );
343  if( value != type.name )
344  return;
345 
346  assert( !success );
347  success = true;
348 
349  f( type );
350  } );
351  if( !success )
352  DUNE_THROW( Dune::Fem::ParameterInvalid, "Parameter '" << key << "' invalid." );
353  }
354 
355  template< class... T, class F >
356  void getEnum ( const std::string &key, const std::tuple< NamedType< T >... > &types, const std::string &defaultValue, F &&f )
357  {
358  getEnum( Dune::Fem::Parameter::container(), key, types, defaultValue, std::forward< F >( f ) );
359  }
360 
361 
362 
363  // makePreconditioner
364  // ------------------
365 
366  template< class AssembledOperator, class Communication >
367  inline std::shared_ptr< Dune::Preconditioner< typename AssembledOperator::domain_type, typename AssembledOperator::range_type > >
368  makePreconditioner ( const std::shared_ptr< AssembledOperator > op, const Communication &comm, Symmetry symmetry )
369  {
370  typedef typename AssembledOperator::matrix_type matrix_type;
371  typedef typename AssembledOperator::domain_type domain_type;
372  typedef typename AssembledOperator::range_type range_type;
373  typedef typename Dune::FieldTraits< typename AssembledOperator::field_type >::real_type real_type;
374 
375  Dune::Amg::DefaultSmootherArgs< real_type > smootherArgs;
376  smootherArgs.relaxationFactor = Dune::Fem::Parameter::getValue( "istl.preconditioner.relax", real_type( 1 ) );
377 
378  const auto smootherTypes = std::make_tuple( NamedType< Dune::SeqJac < matrix_type, domain_type, range_type > >( "jacobi" ),
379  NamedType< Dune::SeqSSOR< matrix_type, domain_type, range_type > >( "ssor" ),
380  NamedType< Dune::SeqILU < matrix_type, domain_type, range_type > >( "ilu" ),
381  NamedType< Dune::SeqILDL< matrix_type, domain_type, range_type > >( "ildl" ) );
382 
383  std::shared_ptr< Dune::Preconditioner< domain_type, range_type > > preconditioner;
384  const std::string preconditionerTypes[] = { "richardson", "smoother", "amg" };
385  switch( Dune::Fem::Parameter::getEnum( "istl.preconditioner.type", preconditionerTypes, 1 ) )
386  {
387  case 0:
388  {
389  auto sp = std::make_shared< Dune::Richardson< domain_type, range_type > >( smootherArgs.relaxationFactor );
390  auto *bp = new Dune::BlockPreconditioner< domain_type, range_type, Communication, std::decay_t< decltype( *sp ) > >( *sp, comm );
391  preconditioner.reset( bp, [ op, sp ] ( decltype( bp ) p ) { delete p; } );
392  }
393  break;
394 
395  case 1:
396  getEnum( "istl.preconditioner.smoother", smootherTypes, "jacobi", [ op, &comm, &smootherArgs, &preconditioner ] ( auto namedType ) {
397  typedef Dune::BlockPreconditioner< domain_type, range_type, Communication, typename decltype( namedType )::Type > SmootherType;
398  typedef Dune::Amg::ConstructionTraits< SmootherType > ConstructionTraits;
399 
400  typename ConstructionTraits::Arguments args;
401  args.setArgs( smootherArgs );
402  args.setComm( comm );
403  args.setMatrix( op->getmat() );
404 
405  preconditioner.reset( ConstructionTraits::construct( args ), [ op ] ( SmootherType *p ) { ConstructionTraits::deconstruct( p ); } );
406  } );
407  break;
408 
409  case 2:
410  getEnum( "istl.preconditioner.smoother", smootherTypes, "jacobi", [ op, &comm, symmetry, &smootherArgs, &preconditioner ] ( auto namedType ) {
411  typedef Dune::BlockPreconditioner< domain_type, range_type, Communication, typename decltype( namedType )::Type > SmootherType;
412 
413  Dune::Amg::Parameters amgParams;
414 
415  // coarsening parameters
416  amgParams.setMaxLevel( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.maxlevel", 100 ) );
417  amgParams.setCoarsenTarget( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.coarsentarget", 1000 ) );
418  amgParams.setMinCoarsenRate( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.mincoarsenrate", 1.2 ) );
419  const std::string accumulationModeNames[] = { "none", "once", "successive" };
420  amgParams.setAccumulate( static_cast< Dune::Amg::AccumulationMode >( Dune::Fem::Parameter::getEnum( "istl.preconditioner.amg.accumulate", accumulationModeNames, 2 ) ) );
421  amgParams.setProlongationDampingFactor( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.prolongation.dampingfactor", 1.6 ) );
422 
423  // parameters
424  amgParams.setDebugLevel( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.debuglevel", 0 ) );
425  amgParams.setNoPreSmoothSteps( Dune::Fem::Parameter::getValue< std::size_t >( "istl.preconditioner.amg.presmoothsteps", 2 ) );
426  amgParams.setNoPostSmoothSteps( Dune::Fem::Parameter::getValue< std::size_t >( "istl.preconditioner.amg.postsmoothsteps", 2 ) );
427  const std::string cycleNames[] = { "v-cycle", "w-cycle" };
428  amgParams.setGamma( 1 + Dune::Fem::Parameter::getEnum( "istl.preconditioner.amg.cycle", cycleNames, 0 ) );
429  amgParams.setAdditive( Dune::Fem::Parameter::getValue( "istl.preconditioner.amg.additive", false ) );
430 
431  typedef Dune::Amg::AMG< AssembledOperator, domain_type, SmootherType, Communication > AMG;
432  if( symmetry == symmetric )
433  {
434  Dune::Amg::CoarsenCriterion< Dune::Amg::SymmetricCriterion< matrix_type, Dune::Amg::RowSum > > criterion( amgParams );
435  preconditioner.reset( new AMG( *op, criterion, smootherArgs, comm ), [ op ] ( AMG *p ) { delete p; } );
436  }
437  else
438  {
439  Dune::Amg::CoarsenCriterion< Dune::Amg::UnSymmetricCriterion< matrix_type, Dune::Amg::RowSum > > criterion( amgParams );
440  preconditioner.reset( new AMG( *op, criterion, smootherArgs, comm ), [ op ] ( AMG *p ) { delete p; } );
441  }
442  } );
443  break;
444 
445  default:
446  DUNE_THROW( Dune::InvalidStateException, "Invalid ISTL preconditioner type selected." );
447  }
448  return preconditioner;
449  }
450 
451 
452 
453  // InverseOperator
454  // ---------------
455 
456  template< class LinearOperator, Symmetry symmetry = unsymmetric >
457  class InverseOperator final
458  : public Dune::Fem::Operator< typename LinearOperator::RangeFunctionType, typename LinearOperator::DomainFunctionType >
459  {
460  static_assert( std::is_same< typename LinearOperator::DomainFunctionType, typename LinearOperator::RangeFunctionType >::value, "Domain function and range function must have the same type." );
461 
462  public:
463  typedef LinearOperator LinearOperatorType;
464  typedef LinearOperatorType OperatorType;
465 
466  typedef typename LinearOperatorType::DomainFunctionType DiscreteFunctionType;
467 
468  typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
469 
470  typedef typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type RealType;
471 
472  typedef SolverParameter SolverParameterType;
473 
474  private:
475  typedef VectorType< DiscreteFunctionType > vector_type;
476  typedef MatrixType< LinearOperatorType > matrix_type;
477 
478  public:
479  // typedef FemCommunication< DiscreteFunctionSpaceType > CommunicationType;
480  typedef Dune::OwnerOverlapCopyCommunication< std::size_t, typename DiscreteFunctionSpaceType::BlockMapperType::GlobalKeyType > CommunicationType;
481 
482  typedef Dune::OverlappingSchwarzOperator< matrix_type, vector_type, vector_type, CommunicationType > AssembledLinearOperatorType;
483  typedef Dune::Preconditioner< vector_type, vector_type > PreconditionerType;
484  typedef Dune::ScalarProduct< vector_type > ScalarProductType;
485 
486  typedef std::function< std::shared_ptr< PreconditionerType > ( std::shared_ptr< AssembledLinearOperatorType >, const CommunicationType &, Symmetry ) > PreconditionerFactory;
487 
488  InverseOperator ( PreconditionerFactory preconditionerFactory, RealType redEps, RealType absLimit, int maxIterations )
489  : preconditionerFactory_( std::move( preconditionerFactory ) ),
490  redEps_( redEps ), absLimit_( absLimit ), maxIterations_( maxIterations )
491  {}
492 
493  InverseOperator ( RealType redEps, RealType absLimit, int maxIterations )
494  : InverseOperator( makePreconditioner< AssembledLinearOperatorType, CommunicationType >, redEps, absLimit, maxIterations )
495  {}
496 
497  InverseOperator ( RealType redEps, RealType absLimit, int maxIterations, bool verbose, const Dune::Fem::SolverParameter& parameter )
498  : InverseOperator( redEps, absLimit, maxIterations )
499  {}
500 
501  InverseOperator ( PreconditionerFactory preconditionerFactory, RealType redEps, RealType absLimit )
502  : InverseOperator( std::move( preconditionerFactory ), redEps, absLimit, std::numeric_limits< int >::max() )
503  {}
504 
505  InverseOperator ( RealType redEps, RealType absLimit )
506  : InverseOperator( redEps, absLimit, std::numeric_limits< int >::max() )
507  {}
508 
509  InverseOperator ( LinearOperator &op, RealType redEps, RealType absLimit, int maxIterations )
510  : InverseOperator( redEps, absLimit, maxIterations )
511  {
512  bind( op );
513  }
514 
515  InverseOperator ( LinearOperator &op, RealType redEps, RealType absLimit )
516  : InverseOperator( op, redEps, absLimit, std::numeric_limits< int >::max() )
517  {}
518 
519  InverseOperator ( const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
520  : InverseOperator( parameter.tolerance(), parameter.tolerance(), parameter.maxIterations() )
521  {}
522 
523  void operator() ( const DiscreteFunctionType &u, DiscreteFunctionType &w ) const override
524  {
525  Dune::InverseOperatorResult result;
526  vector_type b( u.dofVector().array() );
527  if( absLimit_ < std::numeric_limits< RealType >::max() )
528  {
529  vector_type tmp( b );
530  linearOperator_->apply( w.dofVector().array(), tmp );
531  tmp -= b;
532  RealType res = scalarProduct_->norm( tmp );
533  RealType red = (res >0)? absLimit_ / res : 1e-3;
534  solver_->apply( w.dofVector().array(), b, red, result );
535  }
536  else
537  solver_->apply( w.dofVector().array(), b, result );
538  iterations_ = result.iterations;
539  }
540 
541  void bind ( LinearOperator &op )
542  {
543  buildCommunication( op.domainSpace(), Dune::SolverCategory::overlapping, communication_ );
544 
545  linearOperator_.reset( new AssembledLinearOperatorType( op.matrix(), *communication_ ) );
546  scalarProduct_.reset( new Dune::OverlappingSchwarzScalarProduct< vector_type, CommunicationType >( *communication_ ) );
547 
548  // create preconditioner
549  preconditioner_ = preconditionerFactory_( linearOperator_, *communication_, symmetry );
550 
551  // choose absolute error or reduction error
552  const std::string reductionType[] = { "absolute", "relative" };
553  int errorType = Dune::Fem::Parameter::getEnum( "istl.solver.errormeasure", reductionType, 0 );
554  if( errorType != 0 )
556 
557  // create linear solver
558  const std::string verbosityTable[] = { "off", "on", "full" };
559  int verbosity = Dune::Fem::Parameter::getEnum( "istl.solver.verbosity", verbosityTable, 0 );
560  if( op.domainSpace().gridPart().comm().rank() != 0 )
561  verbosity = 0;
562 
563  const std::string solverTypes[] = { "cg", "gcg", "minres", "bicgstab", "gmres" };
564  switch( Dune::Fem::Parameter::getEnum( "istl.solver.type", solverTypes, (symmetry == symmetric ? 0 : 3) ) )
565  {
566  case 0:
567  solver_.reset( new Dune::CGSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
568  break;
569 
570  case 1:
571  {
572  const int restart = Dune::Fem::Parameter::getValue( "istl.solver.restart", 20 );
573  solver_.reset( new Dune::GeneralizedPCGSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity, restart ) );
574  }
575  break;
576 
577  case 2:
578  solver_.reset( new Dune::MINRESSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
579  break;
580 
581  case 3:
582  solver_.reset( new Dune::BiCGSTABSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
583  break;
584 
585  case 4:
586  {
587  const int restart = Dune::Fem::Parameter::getValue( "istl.solver.restart", 20 );
588  solver_.reset( new Dune::RestartedGMResSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, restart, maxIterations_, verbosity ) );
589  }
590  break;
591  }
592  }
593 
594  void unbind ()
595  {
596  solver_.reset();
597  preconditioner_.reset();
598  scalarProduct_.reset();
599  linearOperator_.reset();
600  communication_.reset();
601  }
602 
603  int iterations () const { return iterations_; }
604 
605  void setMaxIterations ( int maxIterations ) { maxIterations_ = maxIterations; }
606 
607  private:
608  PreconditionerFactory preconditionerFactory_;
609  RealType redEps_, absLimit_;
610  int maxIterations_;
611 
612  std::shared_ptr< CommunicationType > communication_;
613  std::shared_ptr< AssembledLinearOperatorType > linearOperator_;
614  std::shared_ptr< ScalarProductType > scalarProduct_;
615  std::shared_ptr< PreconditionerType > preconditioner_;
616  std::shared_ptr< Dune::InverseOperator< vector_type, vector_type > > solver_;
617 
618  mutable int iterations_;
619  };
620 
621  } // namespace ISTL
622 
623  } // namespace Fem
624 
625 } // namespace Dune
626 
627 #endif // #if HAVE_ISTL
628 
629 #endif // #ifndef DUNE_ASH_SOLVER_ISTL_HH
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
static double sqrt(const Double &v)
Definition: double.hh:886
static double max(const Double &v, const double p)
Definition: double.hh:398
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
static constexpr T max(T a)
Definition: utility.hh:77
void buildCommunication(const DiscreteFunctionSpace &dfSpace, Dune::SolverCategory::Category solverCategory, std::shared_ptr< FemCommunication< DiscreteFunctionSpace > > &communication)
Definition: fem.hh:143
Definition: hierarchical/function.hh:106
Definition: io/parameter/exceptions.hh:26
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:159
static T getValue(const std::string &key)
get a mandatory parameter from the container
Definition: io/parameter.hh:333
static int getEnum(const std::string &key, const std::string(&values)[n])
Definition: io/parameter.hh:389
static ParameterContainer & container()
Definition: io/parameter.hh:193
abstract operator
Definition: operator.hh:34
FemCommunication(const DiscreteFunctionSpaceType &dfSpace, Dune::SolverCategory::Category solverCategory=Dune::SolverCategory::sequential)
Definition: fem.hh:80
DiscreteFunctionSpace DiscreteFunctionSpaceType
Definition: fem.hh:78
void dot(const T &x, const T &y, F &scp) const
Definition: fem.hh:106
void project(T &x) const
Definition: fem.hh:96
const DiscreteFunctionSpace::GridPartType::CollectiveCommunicationType & communicator() const
Definition: fem.hh:84
Dune::SolverCategory::Category getSolverCategory() const
Definition: fem.hh:130
Dune::FieldTraits< typename T::field_type >::real_type norm(const T &x) const
Definition: fem.hh:122
void copyOwnerToAll(const T &x, T &y) const
Definition: fem.hh:87
GlobalLookup::LocalIndex::Attribute AttributeType
Definition: owneroverlapcopy.hh:46
void gather(Buffer &buffer, const Entity &entity) const
Definition: owneroverlapcopy.hh:56
BuildRemoteIndicesDataHandle(int rank, const Mapper &mapper, const GlobalLookup &globalLookup)
Definition: owneroverlapcopy.hh:48
std::size_t size(const Entity &entity) const
Definition: owneroverlapcopy.hh:83
bool contains(int dim, int codim) const
Definition: owneroverlapcopy.hh:52
std::map< int, std::vector< Dune::RemoteIndex< GlobalIndexType, AttributeType > > > remotes
Definition: owneroverlapcopy.hh:88
GlobalLookup::GlobalIndex GlobalIndexType
Definition: owneroverlapcopy.hh:45
void scatter(Buffer &buffer, const Entity &entity, std::size_t n)
Definition: owneroverlapcopy.hh:69
bool fixedSize(int dim, int codim) const
Definition: owneroverlapcopy.hh:53
Definition: solver/parameter.hh:15
Definition: parallel.hh:87
discrete function space