1 #ifndef DUNE_FEM_AMGXSOLVER_HH 2 #define DUNE_FEM_AMGXSOLVER_HH 3 4 #include <limits> 5 6 #include <dune/fem/function/common/scalarproducts.hh> 7 #include <dune/fem/operator/common/operator.hh> 8 #include <dune/fem/io/parameter.hh> 9 10 #include <dune/fem/solver/inverseoperatorinterface.hh> 11 #include <dune/fem/solver/parameter.hh> 12 #include <dune/fem/operator/linear/petscoperator.hh> 13 #include <dune/fem/misc/petsc/petsccommon.hh> 14 #include <dune/fem/function/petscdiscretefunction.hh> 15 16 #if HAVE_AMGXSOLVER 17 // AMGX solver wrapper based on Petsc data structures 18 #include <AmgXSolver.hpp> 19 #endif 20 21 namespace Dune 22 { 23 24 namespace Fem 25 { 26 27 //===================================================================== 28 // Implementation of AMGX solver wrapper using PETSc matrix 29 //===================================================================== 30 31 /** @ingroup OEMSolver 32 @{ 33 **/ 34 35 struct AMGXSolverParameter : public LocalParameter< SolverParameter, AMGXSolverParameter > 36 { 37 typedef LocalParameter< SolverParameter, AMGXSolverParameter > BaseType; 38 39 public: 40 using BaseType :: parameter; 41 using BaseType :: keyPrefix; 42 AMGXSolverParameterDune::Fem::AMGXSolverParameter43 AMGXSolverParameter( const ParameterReader& parameter = Parameter::container() ) 44 : BaseType( parameter ) 45 {} 46 47 AMGXSolverParameterDune::Fem::AMGXSolverParameter48 AMGXSolverParameter( const std::string &keyPrefix, const ParameterReader ¶meter = Parameter::container() ) 49 : BaseType( keyPrefix, parameter ) 50 {} 51 AMGXSolverParameterDune::Fem::AMGXSolverParameter52 AMGXSolverParameter( const SolverParameter& sp ) 53 : AMGXSolverParameter( sp.keyPrefix(), sp.parameter() ) 54 {} 55 solvermodeDune::Fem::AMGXSolverParameter56 virtual std::string solvermode () const 57 { 58 const std::string modes [] = { "dDDI" , "dDFI", "dFFI", "hDDI", "hDFI", "hFFI" }; 59 int mode = parameter().getEnum(keyPrefix() + "amgx.mode", modes, 0 ); 60 return modes[ mode ]; 61 } 62 solverconfigDune::Fem::AMGXSolverParameter63 virtual std::string solverconfig () const 64 { 65 return parameter().template getValue< std::string >( keyPrefix() + "amgx.config", "amgxconfig.json"); 66 } 67 }; 68 69 70 71 // AMGXSolver 72 // -------------- 73 74 template< class DiscreteFunction > 75 class AMGXInverseOperator; 76 77 template< class DiscreteFunction > 78 struct AMGXInverseOperatorTraits 79 { 80 typedef DiscreteFunction DiscreteFunctionType; 81 typedef PetscDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType; 82 83 typedef Dune::Fem::Operator< DiscreteFunction, DiscreteFunction > OperatorType; 84 typedef OperatorType PreconditionerType; 85 86 #if HAVE_AMGXSOLVER 87 typedef Fem::PetscLinearOperator< DiscreteFunction, DiscreteFunction > AssembledOperatorType; 88 #else 89 typedef OperatorType AssembledOperatorType; 90 #endif 91 92 typedef AMGXInverseOperator< DiscreteFunction > InverseOperatorType; 93 94 typedef AMGXSolverParameter SolverParameterType; 95 }; 96 97 98 99 100 /** \brief AMGX solver context for PETSc Mat and PETSc Vec */ 101 template< class DF > 102 class AMGXInverseOperator 103 : public InverseOperatorInterface< AMGXInverseOperatorTraits< DF > > 104 { 105 typedef AMGXInverseOperatorTraits< DF > Traits; 106 typedef InverseOperatorInterface< Traits > BaseType; 107 friend class InverseOperatorInterface< Traits >; 108 public: 109 using BaseType :: parameter; 110 111 typedef typename BaseType::SolverDiscreteFunctionType SolverDiscreteFunctionType; 112 typedef typename BaseType::OperatorType OperatorType; 113 typedef typename BaseType::PreconditionerType PreconditionerType; 114 typedef typename BaseType::AssembledOperatorType AssembledOperatorType; 115 116 /** \brief constructor 117 * 118 * \param[in] parameter parameter for the solver 119 */ AMGXInverseOperator(const AMGXSolverParameter & parameter=AMGXSolverParameter ())120 AMGXInverseOperator ( const AMGXSolverParameter& parameter = AMGXSolverParameter() ) 121 : BaseType( parameter ) 122 { 123 } 124 AMGXInverseOperator(const OperatorType & op,const AMGXSolverParameter & parameter=AMGXSolverParameter ())125 AMGXInverseOperator ( const OperatorType &op, 126 const AMGXSolverParameter & parameter = AMGXSolverParameter() ) 127 : AMGXInverseOperator ( parameter ) 128 { 129 bind( op ); 130 } 131 AMGXInverseOperator(const OperatorType & op,PreconditionerType & preconditioner,const AMGXSolverParameter & parameter=AMGXSolverParameter ())132 AMGXInverseOperator ( const OperatorType &op, PreconditionerType &preconditioner, 133 const AMGXSolverParameter & parameter = AMGXSolverParameter() ) 134 : AMGXInverseOperator( parameter ) 135 { 136 bind( op, preconditioner ); 137 } 138 AMGXInverseOperator(const AMGXInverseOperator & other)139 AMGXInverseOperator ( const AMGXInverseOperator& other ) 140 : AMGXInverseOperator( other.parameter() ) 141 { 142 if( other.operator_ ) 143 bind( *(other.operator_) ); 144 } 145 bind(const OperatorType & op)146 void bind( const OperatorType& op ) 147 { 148 BaseType::bind( op ); 149 init( parameter() ); 150 } 151 unbind()152 void unbind() 153 { 154 #if HAVE_AMGXSOLVER 155 amgXSolver_->finalize(); 156 amgXSolver_.reset(); 157 #endif 158 BaseType :: unbind(); 159 } 160 161 protected: init(const AMGXSolverParameter & parameter)162 void init( const AMGXSolverParameter& parameter ) 163 { 164 if( assembledOperator_ ) 165 { 166 std::string mode = parameter.solvermode(); 167 std::string config = parameter.solverconfig(); 168 #if HAVE_AMGXSOLVER 169 amgXSolver_.reset( new AmgXSolver() ); 170 amgXSolver_->initialize(PETSC_COMM_WORLD, mode, config ); 171 172 // check that PetscMat was assembled not in block mode 173 if( assembledOperator_->blockedMode() ) 174 DUNE_THROW(InvalidStateException, "AMGXInverseOperator only works with PetscLinearOperator in non-blocked mode!"); 175 176 // attach Matrix to linear solver context 177 Mat& A = const_cast<Mat &> (assembledOperator_->exportMatrix()); 178 179 // set matrix 180 amgXSolver_->setA( A ); 181 #else 182 DUNE_THROW(InvalidStateException,"AMGX solver or PETSc not found during cmake config. Please reconfigure!"); 183 #endif 184 } 185 } 186 apply(const SolverDiscreteFunctionType & arg,SolverDiscreteFunctionType & dest) const187 int apply( const SolverDiscreteFunctionType& arg, SolverDiscreteFunctionType& dest ) const 188 { 189 if( !assembledOperator_ ) 190 DUNE_THROW(NotImplemented,"AMGX solver with matrix free implementations is not supported!"); 191 192 193 int iterations = -1; 194 #if HAVE_AMGXSOLVER 195 assert( amgXSolver_ ); 196 197 // need to have a 'distributed' destination vector for continuous spaces 198 if( dest.space().continuous() ) 199 dest.dofVector().clearGhost(); 200 201 // call PETSc solvers, dest = x, arg = rhs 202 Vec& x = *dest.petscVec(); 203 Vec& rhs = *(const_cast< SolverDiscreteFunctionType& > (arg).petscVec()); 204 amgXSolver_->solve( x, rhs ); 205 206 // a continuous solution is 'distributed' so need a communication here 207 if( dest.space().continuous() ) 208 { 209 dest.communicate(); 210 } 211 212 // get number of iterations 213 amgXSolver_->getIters( iterations ); 214 #else 215 DUNE_THROW(InvalidStateException,"AMGX solver or PETSc not found during cmake config. Please reconfigure!"); 216 #endif 217 return iterations; 218 } 219 220 protected: 221 #if HAVE_AMGXSOLVER 222 mutable std::unique_ptr< AmgXSolver > amgXSolver_; 223 #endif 224 using BaseType :: assembledOperator_; 225 using BaseType :: parameter_; 226 }; 227 228 ///@} 229 230 } // namespace Fem 231 232 } // namespace Dune 233 234 #endif // #ifndef DUNE_FEM_PETSCSOLVER_HH 235