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 &parameter = 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