1 #ifndef DUNE_PYTHON_ISTL_SOLVER_HH 2 #define DUNE_PYTHON_ISTL_SOLVER_HH 3 4 #include <dune/common/typeutilities.hh> 5 #include <dune/common/version.hh> 6 7 #include <dune/istl/solver.hh> 8 #include <dune/istl/solvers.hh> 9 #include <dune/istl/preconditioners.hh> 10 11 #include <dune/python/istl/preconditioners.hh> 12 13 #include <dune/python/pybind11/pybind11.h> 14 15 namespace Dune 16 { 17 18 namespace Python 19 { 20 21 // registerInverseOperator 22 // ----------------------- 23 24 template< class Solver, class... options > registerInverseOperator(pybind11::class_<Solver,options...> cls)25 inline void registerInverseOperator ( pybind11::class_< Solver, options... > cls ) 26 { 27 typedef typename Solver::domain_type Domain; 28 typedef typename Solver::range_type Range; 29 30 using pybind11::operator""_a; 31 32 cls.def( "__call__", [] ( Solver &self, Domain &x, Range &b, double reduction ) { 33 InverseOperatorResult result; 34 self.apply( x, b, reduction, result ); 35 return std::make_tuple( result.iterations, result.reduction, result.converged, result.conv_rate, result.elapsed ); 36 }, "x"_a, "b"_a, "reduction"_a, 37 R"doc( 38 Solve linear system 39 40 Args: 41 x: solution of linear system 42 b: right hand side of the system 43 reduction: factor to reduce the defect by 44 45 Returns: (iterations, reduction, converged, conv_rate, elapsed) 46 iterations: number of iterations performed 47 reduction: actual factor, the error has been reduced by 48 converged: True, if the solver has achieved its reduction requirements 49 conv_rate: rate of convergence 50 elapsed: time in seconds used to solve the linear system 51 52 Note: 53 - If the reduction is omitted, the default value of the solver is used. 54 - For iterative solvers, the solution must be initialized to the starting point. 55 - The right hand side b will be replaced by the residual. 56 )doc" ); 57 58 cls.def( "__call__", [] ( Solver &self, Domain &x, Range &b ) { 59 InverseOperatorResult result; 60 self.apply( x, b, result ); 61 return std::make_tuple( result.iterations, result.reduction, result.converged, result.conv_rate, result.elapsed ); 62 }, "x"_a, "b"_a ); 63 64 cls.def_property_readonly( "category", [] ( const Solver &self ) { return self.category(); }, 65 R"doc( 66 Obtain category of the linear solver 67 )doc" ); 68 69 cls.def( "asPreconditioner", [] ( Solver &self ) { 70 return new InverseOperator2Preconditioner< Solver >( self ); 71 }, pybind11::keep_alive< 0, 1 >(), 72 R"doc( 73 Convert linear solver into preconditioner 74 )doc" ); 75 } 76 77 78 79 namespace detail 80 { 81 82 // registerEndomorphismSolvers 83 // --------------------------- 84 85 template< class X, class Y, class... options > 86 inline std::enable_if_t< std::is_same< X, Y >::value > registerEndomorphismSolvers(pybind11::module module,pybind11::class_<LinearOperator<X,Y>,options...>,PriorityTag<1>)87 registerEndomorphismSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... >, PriorityTag< 1 > ) 88 { 89 typedef Dune::InverseOperator< X, Y > Solver; 90 91 using pybind11::operator""_a; 92 93 pybind11::options opts; 94 opts.disable_function_signatures(); 95 96 module.def( "LoopSolver", [] ( LinearOperator< X, X > &op, Preconditioner< X, X > &prec, double reduction, int maxit, int verbose ) { 97 return static_cast< Solver * >( new Dune::LoopSolver< X >( op, prec, reduction, maxit, verbose ) ); 98 }, "operator"_a, "preconditioner"_a, "reduction"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 99 R"doc( 100 Loop solver 101 102 Args: 103 operator: operator to invert 104 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 105 reduction: factor to reduce the defect by 106 maxIterations: maximum number of iterations to perform 107 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 108 109 Returns: 110 ISTL Loop solver 111 112 Note: 113 The loop solver will apply the preconditioner once in each step. 114 )doc" ); 115 116 module.def( "GradientSolver", [] ( LinearOperator< X, X > &op, Preconditioner< X, X > &prec, double reduction, int maxit, int verbose ) { 117 return static_cast< Solver * >( new Dune::GradientSolver< X >( op, prec, reduction, maxit, verbose ) ); 118 }, "operator"_a, "preconditioner"_a, "reduction"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 119 R"doc( 120 Gradient iterative solver 121 122 Args: 123 operator: operator to invert 124 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 125 reduction: factor to reduce the defect by 126 maxIterations: maximum number of iterations to perform 127 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 128 129 Returns: 130 ISTL Gradient solver 131 132 Note: 133 This method is also know as steepest descend method. 134 )doc" ); 135 136 module.def( "CGSolver", [] ( LinearOperator< X, X > &op, Preconditioner< X, X > &prec, double reduction, int maxit, int verbose ) { 137 return static_cast< Solver * >( new Dune::CGSolver< X >( op, prec, reduction, maxit, verbose ) ); 138 }, "operator"_a, "preconditioner"_a, "reduction"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 139 R"doc( 140 Conjugate gradient iterative solver 141 142 Args: 143 operator: operator to invert 144 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 145 reduction: factor to reduce the defect by 146 maxIterations: maximum number of iterations to perform 147 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 148 149 Returns: 150 ISTL Conjugate gradient solver 151 152 Note: 153 The conjucate gradient method can only be applied if the operator and the preconditioner are both symmetric and positive definite. 154 )doc" ); 155 156 module.def( "BiCGSTABSolver", [] ( LinearOperator< X, X > &op, Preconditioner< X, X > &prec, double reduction, int maxit, int verbose ) { 157 return static_cast< Solver * >( new Dune::CGSolver< X >( op, prec, reduction, maxit, verbose ) ); 158 }, "operator"_a, "preconditioner"_a, "reduction"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 159 R"doc( 160 Biconjugate gradient stabilized iterative solver 161 162 Args: 163 operator: operator to invert 164 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 165 reduction: factor to reduce the defect by 166 maxIterations: maximum number of iterations to perform 167 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 168 169 Returns: 170 ISTL Biconjugate gradient stabilized solver 171 )doc" ); 172 173 module.def( "MinResSolver", [] ( LinearOperator< X, X > &op, Preconditioner< X, X > &prec, double reduction, int maxit, int verbose ) { 174 return static_cast< Solver * >( new Dune::CGSolver< X >( op, prec, reduction, maxit, verbose ) ); 175 }, "operator"_a, "preconditioner"_a, "reduction"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 176 R"doc( 177 Minimal residual iterative solver 178 179 Args: 180 operator: operator to invert 181 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 182 reduction: factor to reduce the defect by 183 maxIterations: maximum number of iterations to perform 184 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 185 186 Returns: 187 ISTL Minimal residual solver 188 189 Note: 190 The minimal residual method can only be applied if the operator and the preconditioner are both symmetric. 191 )doc" ); 192 } 193 194 template< class X, class Y, class... options > 195 inline std::enable_if_t< std::is_same< X, Y >::value > registerEndomorphismSolvers(pybind11::module module,pybind11::class_<LinearOperator<X,Y>,options...>,PriorityTag<0>)196 registerEndomorphismSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... >, PriorityTag< 0 > ) 197 {} 198 199 template< class X, class Y, class... options > registerEndomorphismSolvers(pybind11::module module,pybind11::class_<LinearOperator<X,Y>,options...> cls)200 inline void registerEndomorphismSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... > cls ) 201 { 202 registerEndomorphismSolvers( module, cls, PriorityTag< 42 >() ); 203 } 204 205 } // namespace detail 206 207 208 209 // registerSolvers 210 // --------------- 211 212 template< class X, class Y, class... options > registerSolvers(pybind11::module module,pybind11::class_<LinearOperator<X,Y>,options...> cls)213 inline void registerSolvers ( pybind11::module module, pybind11::class_< LinearOperator< X, Y >, options... > cls ) 214 { 215 typedef Dune::InverseOperator< X, Y > Solver; 216 217 using pybind11::operator""_a; 218 219 pybind11::options opts; 220 opts.disable_function_signatures(); 221 222 pybind11::class_< Solver > clsSolver( module, "InverseOperator" ); 223 registerInverseOperator( clsSolver ); 224 225 detail::registerEndomorphismSolvers( module, cls ); 226 227 module.def( "RestartedGMResSolver", [] ( LinearOperator< X, Y > &op, Preconditioner< X, Y > &prec, double reduction, int restart, int maxit, int verbose ) { 228 return static_cast< Solver * >( new Dune::RestartedGMResSolver< X, Y >( op, prec, reduction, restart, maxit, verbose ) ); 229 }, "operator"_a, "preconditioner"_a, "reduction"_a, "restart"_a, "maxIterations"_a = std::numeric_limits< int >::max(), "verbose"_a = 0, pybind11::keep_alive< 0, 1 >(), pybind11::keep_alive< 0, 2 >(), 230 R"doc( 231 Restarted generalized minimal residual iterative solver 232 233 Args: 234 operator: operator to invert 235 preconditioner: preconditioner to use (i.e., apprixmate inverse of the operator) 236 reduction: factor to reduce the defect by 237 restart: number of iterations before restart 238 maxIterations: maximum number of iterations to perform 239 verbose: verbosity level (0 = quiet, 1 = summary, 2 = verbose) 240 241 Returns: 242 ISTL Restarted generalized minimal residual solver 243 244 Note: 245 The restarted generalized minimal residual method holds restart many vectors in memory during application. 246 This can lead to a large memory consumption. 247 )doc" ); 248 } 249 250 } // namespace Python 251 252 } // namespace Dune 253 254 #endif // #ifndef DUNE_PYTHON_ISTL_SOLVER_HH 255