1#ifndef DUNE_ASH_SOLVER_ISTL_HH
2#define DUNE_ASH_SOLVER_ISTL_HH
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>
19#include <dune/geometry/dimension.hh>
21#include <dune/grid/common/datahandleif.hh>
22#include <dune/grid/common/rangegenerators.hh>
23#include <dune/grid/common/partitionset.hh>
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>
45 template<
class M,
class X,
class Y >
46 struct ConstructionTraits< SeqILDL< M, X, Y > >
48 typedef DefaultConstructionArgs< SeqILDL< M, X, Y > > Arguments;
50 static SeqILDL< M, X, Y > *construct ( Arguments &args )
52 return new SeqILDL< M, X, Y >( args.getMatrix(), args.getArgs().relaxationFactor );
55 static void deconstruct ( SeqILDL< M, X, Y > *p ) {
delete p; }
68 template<
class DiscreteFunctionSpace >
69 class HierarchicalDiscreteFunction;
71 template<
class DomainFunction,
class RangeFunction >
72 class ISTLLinearOperator;
87 template<
class DiscreteFunction >
88 using VectorType = std::decay_t< decltype( std::declval< const DiscreteFunction & >().dofVector().array() ) >;
95 template<
class LinearOperator >
98 typedef std::decay_t< decltype( std::declval< const LinearOperator & >().matrix() ) > Type;
101 template<
class DomainFunction,
class RangeFunction >
102 struct __MatrixType<
Dune::Fem::ISTLLinearOperator< DomainFunction, RangeFunction > >
104 typedef Dune::BCRSMatrix< typename Dune::Fem::ISTLLinearOperator< DomainFunction, RangeFunction >::LittleBlockType > Type;
108 template<
class LinearOperator >
109 using MatrixType =
typename __MatrixType< LinearOperator >::Type;
116 enum Symmetry :
bool { unsymmetric =
false, symmetric =
true };
123 template<
class DiscreteFunctionSpace >
124 class FemCommunication
126 typedef FemCommunication< DiscreteFunctionSpace > ThisType;
131 typedef int GlobalLookupIndexSet;
134 : dfSpace_( dfSpace ), solverCategory_( solverCategory )
137 const typename DiscreteFunctionSpace::GridPartType::CollectiveCommunicationType &
communicator ()
const {
return dfSpace_.gridPart().comm(); }
144 DiscreteFunctionType z(
"", dfSpace_,
reinterpret_cast< typename DiscreteFunctionType::DofVectorType &
>( y ) );
151 typedef typename T::field_type field_type;
154 const auto &auxiliaryDofs = dfSpace_.auxiliaryDofs();
155 for(
int i : auxiliaryDofs )
156 x[ i ] = field_type( 0 );
159 template<
class T,
class F >
160 void dot (
const T &x,
const T &y, F &scp )
const
162 const auto &auxiliaryDofs = dfSpace_.auxiliaryDofs();
164 const int numAuxiliarys = auxiliaryDofs.size();
165 for(
int auxiliary = 0, i = 0; auxiliary < numAuxiliarys; ++auxiliary, ++i )
167 const int nextAuxiliary = auxiliaryDofs[ auxiliary ];
168 for( ; i < nextAuxiliary; ++i )
169 scp += x[ i ] * y[ i ];
176 typename Dune::FieldTraits< typename T::field_type >::real_type
norm (
const T &x )
const
179 typename Dune::FieldTraits< typename T::field_type >::real_type norm2( 0 );
181 return sqrt( norm2 );
184 Dune::SolverCategory::Category
getSolverCategory ()
const {
return solverCategory_; }
188 Dune::SolverCategory::Category solverCategory_;
196 template<
class Mapper,
class GlobalLookup >
197 struct BuildRemoteIndicesDataHandle
198 :
public Dune::CommDataHandleIF< BuildRemoteIndicesDataHandle< Mapper, GlobalLookup >, int >
201 typedef typename GlobalLookup::LocalIndex::Attribute
AttributeType;
204 : rank_( rank ), mapper_( mapper ), globalLookup_( globalLookup )
207 bool contains (
int dim,
int codim )
const {
return mapper_.contains( codim ); }
208 bool fixedSize(
int dim,
int codim )
const {
return true; }
210 template<
class Buffer,
class Entity >
211 void gather ( Buffer &buffer,
const Entity &entity )
const
213 buffer.write( rank_ );
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();
220 buffer.write(
static_cast< int >( attribute ) );
223 template<
class Buffer,
class Entity >
224 void scatter ( Buffer &buffer,
const Entity &entity, std::size_t n )
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 );
237 template<
class Entity >
238 std::size_t
size (
const Entity &entity )
const
243 std::map< int, std::vector< Dune::RemoteIndex< GlobalIndexType, AttributeType > > >
remotes;
247 const Mapper &mapper_;
248 const GlobalLookup &globalLookup_;
256 template<
class DiscreteFunction >
257 void buildCommunication (
const typename DiscreteFunction::DiscreteFunctionSpaceType &dfSpace,
258 Dune::SolverCategory::Category solverCategory,
259 std::shared_ptr< FemCommunication< DiscreteFunction > > &communication )
261 communication.reset(
new FemCommunication< DiscreteFunction >( dfSpace, solverCategory ) );
264 template<
class DiscreteFunctionSpace,
class GlobalId,
class LocalId >
266 Dune::SolverCategory::Category solverCategory,
267 std::shared_ptr< Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId > > &communication )
269 typedef typename DiscreteFunctionSpace::GridPartType GridPartType;
270 typedef typename DiscreteFunctionSpace::BlockMapperType LocalMapperType;
272 typedef typename Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId >::GlobalLookupIndexSet GlobalLookupType;
274 typedef typename GlobalLookupType::LocalIndex LocalIndexType;
276 communication.reset(
new Dune::OwnerOverlapCopyCommunication< GlobalId, LocalId >( solverCategory ) );
278 const GridPartType &gridPart = dfSpace.gridPart();
279 LocalMapperType &localMapper = dfSpace.blockMapper();
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;
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();
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();
301 communication->remoteIndices().setIndexSets( communication->indexSet(), communication->indexSet(), communication->communicator() );
302 if( !buildRemoteIndicesDataHandle.remotes.empty() )
304 for(
auto &remote : buildRemoteIndicesDataHandle.remotes )
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 );
313 communication->remoteIndices().template getModifier< false, true >( 0 );
326 NamedType ( std::string name ) : name( name ) {}
336 template<
class... T,
class F >
337 void getEnum (
const Dune::Fem::ParameterReader ¶meter,
const std::string &key,
const std::tuple< NamedType< T >... > &types,
const std::string &defaultValue, F &&f )
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 )
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 )
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 )
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;
375 Dune::Amg::DefaultSmootherArgs< real_type > smootherArgs;
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" ) );
383 std::shared_ptr< Dune::Preconditioner< domain_type, range_type > > preconditioner;
384 const std::string preconditionerTypes[] = {
"richardson",
"smoother",
"amg" };
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; } );
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;
400 typename ConstructionTraits::Arguments args;
401 args.setArgs( smootherArgs );
402 args.setComm( comm );
403 args.setMatrix( op->getmat() );
405 preconditioner.reset( ConstructionTraits::construct( args ), [ op ] ( SmootherType *p ) { ConstructionTraits::deconstruct( p ); } );
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;
413 Dune::Amg::Parameters amgParams;
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 ) ) );
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" };
431 typedef Dune::Amg::AMG< AssembledOperator, domain_type, SmootherType, Communication > AMG;
432 if( symmetry == symmetric )
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; } );
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; } );
446 DUNE_THROW( Dune::InvalidStateException,
"Invalid ISTL preconditioner type selected." );
448 return preconditioner;
456 template<
class LinearOperator, Symmetry symmetry = unsymmetric >
457 class InverseOperator final
458 :
public Dune::Fem::Operator< typename LinearOperator::RangeFunctionType, typename LinearOperator::DomainFunctionType >
460 static_assert( std::is_same< typename LinearOperator::DomainFunctionType, typename LinearOperator::RangeFunctionType >::value,
"Domain function and range function must have the same type." );
463 typedef LinearOperator LinearOperatorType;
464 typedef LinearOperatorType OperatorType;
466 typedef typename LinearOperatorType::DomainFunctionType DiscreteFunctionType;
468 typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
470 typedef typename Dune::FieldTraits< typename DiscreteFunctionType::RangeFieldType >::real_type RealType;
472 typedef SolverParameter SolverParameterType;
475 typedef VectorType< DiscreteFunctionType > vector_type;
476 typedef MatrixType< LinearOperatorType > matrix_type;
480 typedef Dune::OwnerOverlapCopyCommunication< std::size_t, typename DiscreteFunctionSpaceType::BlockMapperType::GlobalKeyType > CommunicationType;
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;
486 typedef std::function< std::shared_ptr< PreconditionerType > ( std::shared_ptr< AssembledLinearOperatorType >,
const CommunicationType &, Symmetry ) > PreconditionerFactory;
488 InverseOperator ( PreconditionerFactory preconditionerFactory, RealType redEps, RealType absLimit,
int maxIterations )
489 : preconditionerFactory_(
std::move( preconditionerFactory ) ),
490 redEps_( redEps ), absLimit_( absLimit ), maxIterations_( maxIterations )
493 InverseOperator ( RealType redEps, RealType absLimit,
int maxIterations )
494 : InverseOperator( makePreconditioner< AssembledLinearOperatorType, CommunicationType >, redEps, absLimit, maxIterations )
497 InverseOperator ( RealType redEps, RealType absLimit,
int maxIterations,
bool verbose,
const Dune::Fem::SolverParameter& parameter )
498 : InverseOperator( redEps, absLimit, maxIterations )
501 InverseOperator ( PreconditionerFactory preconditionerFactory, RealType redEps, RealType absLimit )
502 : InverseOperator(
std::move( preconditionerFactory ), redEps, absLimit,
std::numeric_limits< int >::
max() )
505 InverseOperator ( RealType redEps, RealType absLimit )
506 : InverseOperator( redEps, absLimit,
std::numeric_limits< int >::
max() )
509 InverseOperator ( LinearOperator &op, RealType redEps, RealType absLimit,
int maxIterations )
510 : InverseOperator( redEps, absLimit, maxIterations )
515 InverseOperator ( LinearOperator &op, RealType redEps, RealType absLimit )
516 : InverseOperator( op, redEps, absLimit,
std::numeric_limits< int >::
max() )
519 InverseOperator (
const SolverParameter& parameter = SolverParameter( Parameter::container() ) )
520 : InverseOperator( parameter.tolerance(), parameter.tolerance(), parameter.maxIterations() )
523 void operator() (
const DiscreteFunctionType &u, DiscreteFunctionType &w )
const override
525 Dune::InverseOperatorResult result;
526 vector_type b( u.dofVector().array() );
527 if( absLimit_ < std::numeric_limits< RealType >::max() )
529 vector_type tmp( b );
530 linearOperator_->apply( w.dofVector().array(), tmp );
532 RealType res = scalarProduct_->norm( tmp );
533 RealType red = (res >0)? absLimit_ / res : 1e-3;
534 solver_->apply( w.dofVector().array(), b, red, result );
537 solver_->apply( w.dofVector().array(), b, result );
538 iterations_ = result.iterations;
541 void bind ( LinearOperator &op )
543 buildCommunication( op.domainSpace(), Dune::SolverCategory::overlapping, communication_ );
545 linearOperator_.reset(
new AssembledLinearOperatorType( op.matrix(), *communication_ ) );
546 scalarProduct_.reset(
new Dune::OverlappingSchwarzScalarProduct< vector_type, CommunicationType >( *communication_ ) );
549 preconditioner_ = preconditionerFactory_( linearOperator_, *communication_, symmetry );
552 const std::string reductionType[] = {
"absolute",
"relative" };
555 absLimit_ = std::numeric_limits< RealType >::max();
558 const std::string verbosityTable[] = {
"off",
"on",
"full" };
560 if( op.domainSpace().gridPart().comm().rank() != 0 )
563 const std::string solverTypes[] = {
"cg",
"gcg",
"minres",
"bicgstab",
"gmres" };
567 solver_.reset(
new Dune::CGSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
573 solver_.reset(
new Dune::GeneralizedPCGSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity, restart ) );
578 solver_.reset(
new Dune::MINRESSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
582 solver_.reset(
new Dune::BiCGSTABSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, maxIterations_, verbosity ) );
588 solver_.reset(
new Dune::RestartedGMResSolver< vector_type >( *linearOperator_, *scalarProduct_, *preconditioner_, redEps_, restart, maxIterations_, verbosity ) );
597 preconditioner_.reset();
598 scalarProduct_.reset();
599 linearOperator_.reset();
600 communication_.reset();
603 int iterations ()
const {
return iterations_; }
605 void setMaxIterations (
int maxIterations ) { maxIterations_ = maxIterations; }
608 PreconditionerFactory preconditionerFactory_;
609 RealType redEps_, absLimit_;
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_;
618 mutable int iterations_;
double sqrt(const Dune::Fem::Double &v)
Definition: double.hh:977
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
static double sqrt(const Double &v)
Definition: double.hh:886
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
void buildCommunication(const DiscreteFunctionSpace &dfSpace, Dune::SolverCategory::Category solverCategory, std::shared_ptr< FemCommunication< DiscreteFunctionSpace > > &communication)
Definition: fem.hh:143
Definition: hierarchical/function.hh:106
static ParameterContainer & container()
Definition: io/parameter.hh:193
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
Definition: io/parameter/exceptions.hh:26
T getValue(const std::string &key) const
get mandatory parameter
Definition: reader.hh:159
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
const DiscreteFunctionSpace::GridPartType::CollectiveCommunicationType & communicator() const
Definition: fem.hh:84
Dune::FieldTraits< typenameT::field_type >::real_type norm(const T &x) const
Definition: fem.hh:122
void project(T &x) const
Definition: fem.hh:96
Dune::SolverCategory::Category getSolverCategory() const
Definition: fem.hh:130
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