 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11
12 /// This HYPRE library interface has been taken originally from MFEM and modified
13 /// to suit the needs for the MOAB library.
14 /// Modified by: Vijay Mahadevan
15
16 #ifndef MOAB_HYPRESOLVER
17 #define MOAB_HYPRESOLVER
18
19 #include "moab/MOABConfig.h"
20 #include "moab/Core.hpp"
21
22 #ifdef MOAB_HAVE_MPI
23 #include <mpi.h>
24
25 // Enable internal hypre timing routines
26 #ifndef HYPRE_TIMING
27 #define HYPRE_TIMING
28 #endif
29
30 #include "HypreParVector.hpp"
31 #include "HypreParMatrix.hpp"
32
33 namespace moab
34 {
35
36 /// Base class for solvers
37 class AbstractSolver
38 {
39 protected:
40 const Eigen::SparseMatrix< double >* m_operator;
41
42 public:
43 /// If true, use the second argument of Mult as an initial guess
44 bool iterative_mode;
45
46 /// Initialize a Solver
47 explicit AbstractSolver( bool iter_mode )
48 {
49 iterative_mode = iter_mode;
50 }
51
52 /// Set/update the solver for the given operator
53 virtual void SetOperator( const HypreParMatrix& op ) = 0;
54 };
55
56 /// Abstract class for hypre's solvers and preconditioners
57 class HypreSolver : public AbstractSolver
58 {
59 public:
60 typedef HypreParVector Vector;
61
62 protected:
63 /// The linear system matrix
64 HypreParMatrix* A;
65
66 /// Right-hand side and solution vector
67 mutable HypreParVector *B, *X;
68
69 /// Was hypre's Setup function called already?
70 mutable int setup_called;
71
72 public:
73 HypreSolver( bool iterative = true );
74
75 HypreSolver( HypreParMatrix* _A, bool iterative = true );
76
77 /// Typecast to HYPRE_Solver -- return the solver
78 virtual operator HYPRE_Solver() const = 0;
79
80 virtual void Verbosity( int /*print_level*/ )
81 { /* nothing to do */
82 }
83
84 /// Set the hypre solver to be used as a preconditioner
85 virtual void SetPreconditioner( HypreSolver& /*precond*/ )
86 { /* Nothing to do */
87 }
88
89 /// hypre's internal Setup function
90 virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
91 /// hypre's internal Solve function
92 virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;
93
94 virtual void SetOperator( const HypreParMatrix& /*op*/ )
95 {
96 MB_SET_ERR_RET( "HypreSolvers do not support SetOperator!" );
97 }
98
99 virtual void GetFinalResidualNorm( double& normr ) = 0;
100
101 virtual void GetNumIterations( int& num_iterations ) = 0;
102
103 /// Solve the linear system Ax=b
104 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
105
106 virtual ~HypreSolver();
107 };
108
109 /// PCG solver in hypre
110 class HyprePCG : public HypreSolver
111 {
112 private:
113 HYPRE_Solver pcg_solver;
114
115 public:
116 HyprePCG( HypreParMatrix& _A );
117
118 void SetTol( double tol );
119 void SetMaxIter( int max_iter );
120 void SetLogging( int logging );
121 virtual void Verbosity( int print_lvl );
122
123 /// Set the hypre solver to be used as a preconditioner
124 virtual void SetPreconditioner( HypreSolver& precond );
125
126 /** Use the L2 norm of the residual for measuring PCG convergence, plus
127 (optionally) 1) periodically recompute true residuals from scratch; and
128 2) enable residual-based stopping criteria. */
129 void SetResidualConvergenceOptions( int res_frequency = -1, double rtol = 0.0 );
130
131 /// non-hypre setting
132 void SetZeroInintialIterate()
133 {
134 iterative_mode = false;
135 }
136
137 virtual void GetFinalResidualNorm( double& normr )
138 {
139 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( pcg_solver, &normr );
140 }
141
142 virtual void GetNumIterations( int& num_iterations )
143 {
144 HYPRE_Int num_it;
145 HYPRE_ParCSRPCGGetNumIterations( pcg_solver, &num_it );
146 num_iterations = internal::to_int( num_it );
147 }
148
149 /// The typecast to HYPRE_Solver returns the internal pcg_solver
150 virtual operator HYPRE_Solver() const
151 {
152 return pcg_solver;
153 }
154
155 /// PCG Setup function
156 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
157 {
158 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSetup;
159 }
160 /// PCG Solve function
161 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
162 {
163 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRPCGSolve;
164 }
165
166 /// Solve Ax=b with hypre's PCG
167 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
168 using HypreSolver::Solve;
169
170 virtual ~HyprePCG();
171 };
172
173 /// GMRES solver in hypre
174 class HypreGMRES : public HypreSolver
175 {
176 private:
177 HYPRE_Solver gmres_solver;
178
179 public:
180 HypreGMRES( HypreParMatrix& _A );
181
182 void SetTol( double atol, double rtol = 1e-20 );
183 void SetMaxIter( int max_iter );
184 void SetKDim( int dim );
185 void SetLogging( int logging );
186 virtual void Verbosity( int print_lvl );
187
188 /// Set the hypre solver to be used as a preconditioner
189 virtual void SetPreconditioner( HypreSolver& precond );
190
191 /** Use the L2 norm of the residual for measuring PCG convergence, plus
192 (optionally) 1) periodically recompute true residuals from scratch; and
193 2) enable residual-based stopping criteria. */
194 void SkipRealResidualCheck()
195 {
196 HYPRE_GMRESSetSkipRealResidualCheck( gmres_solver, 1 );
197 }
198
199 /// non-hypre setting
200 void SetZeroInintialIterate()
201 {
202 iterative_mode = false;
203 }
204
205 virtual void GetFinalResidualNorm( double& normr )
206 {
207 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( gmres_solver, &normr );
208 }
209
210 virtual void GetNumIterations( int& num_iterations )
211 {
212 HYPRE_Int num_it;
213 HYPRE_GMRESGetNumIterations( gmres_solver, &num_it );
214 num_iterations = internal::to_int( num_it );
215 }
216
217 /// The typecast to HYPRE_Solver returns the internal gmres_solver
218 virtual operator HYPRE_Solver() const
219 {
220 return gmres_solver;
221 }
222
223 /// GMRES Setup function
224 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
225 {
226 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSetup;
227 }
228 /// GMRES Solve function
229 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
230 {
231 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRGMRESSolve;
232 }
233
234 /// Solve Ax=b with hypre's GMRES
235 virtual HYPRE_Int Solve( const HypreParVector& b, HypreParVector& x ) const;
236 using HypreSolver::Solve;
237
238 virtual ~HypreGMRES();
239 };
240
241 /// The identity operator as a hypre solver
242 class HypreIdentity : public HypreSolver
243 {
244 public:
245 virtual operator HYPRE_Solver() const
246 {
247 return NULL;
248 }
249
250 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
251 {
252 return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentitySetup;
253 }
254 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
255 {
256 return (HYPRE_PtrToParSolverFcn)hypre_ParKrylovIdentity;
257 }
258
259 virtual void GetNumIterations( int& num_iterations )
260 {
261 num_iterations = 1;
262 }
263
264 virtual void GetFinalResidualNorm( double& normr )
265 {
266 normr = 0.0;
267 }
268
269 virtual ~HypreIdentity() {}
270 };
271
272 /// Jacobi preconditioner in hypre
273 class HypreDiagScale : public HypreSolver
274 {
275 public:
276 HypreDiagScale() : HypreSolver() {}
277 explicit HypreDiagScale( HypreParMatrix& A ) : HypreSolver( &A ) {}
278 virtual operator HYPRE_Solver() const
279 {
280 return NULL;
281 }
282
283 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
284 {
285 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScaleSetup;
286 }
287 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
288 {
289 return (HYPRE_PtrToParSolverFcn)HYPRE_ParCSRDiagScale;
290 }
291
292 virtual void GetNumIterations( int& num_iterations )
293 {
294 num_iterations = 1;
295 }
296
297 virtual void GetFinalResidualNorm( double& normr )
298 {
299 normr = 0.0;
300 }
301
302 HypreParMatrix* GetData()
303 {
304 return A;
305 }
306 virtual ~HypreDiagScale() {}
307 };
308
309 /// The ParaSails preconditioner in hypre
310 class HypreParaSails : public HypreSolver
311 {
312 private:
313 HYPRE_Solver sai_precond;
314
315 public:
316 HypreParaSails( HypreParMatrix& A );
317
318 void SetSymmetry( int sym );
319
320 virtual void GetNumIterations( int& num_iterations )
321 {
322 num_iterations = 1;
323 }
324
325 virtual void GetFinalResidualNorm( double& normr )
326 {
327 normr = 0.0;
328 }
329
330 /// The typecast to HYPRE_Solver returns the internal sai_precond
331 virtual operator HYPRE_Solver() const
332 {
333 return sai_precond;
334 }
335
336 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
337 {
338 return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSetup;
339 }
340 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
341 {
342 return (HYPRE_PtrToParSolverFcn)HYPRE_ParaSailsSolve;
343 }
344
345 virtual ~HypreParaSails();
346 };
347
348 /// The BoomerAMG solver in hypre
349 class HypreBoomerAMG : public HypreSolver
350 {
351 private:
352 HYPRE_Solver amg_precond;
353
354 /// Default, generally robust, BoomerAMG options
355 void SetDefaultOptions();
356
357 // If amg_precond is NULL, allocates it and sets default options.
358 // Otherwise saves the options from amg_precond, destroys it, allocates a new
359 // one, and sets its options to the saved values.
360 void ResetAMGPrecond();
361
362 public:
363 HypreBoomerAMG();
364
365 HypreBoomerAMG( HypreParMatrix& A );
366
367 virtual void SetOperator( const HypreParMatrix& op );
368
369 virtual void GetNumIterations( int& num_iterations )
370 {
371 num_iterations = 1;
372 }
373
374 virtual void GetFinalResidualNorm( double& normr )
375 {
376 normr = 0.0;
377 }
378
379 /** More robust options for systems, such as elasticity. Note that BoomerAMG
380 assumes Ordering::byVDIM in the finite element space used to generate the
381 matrix A. */
382 void SetSystemsOptions( int dim );
383
384 virtual void Verbosity( int print_level )
385 {
386 HYPRE_BoomerAMGSetPrintLevel( amg_precond, print_level );
387 }
388
389 /// The typecast to HYPRE_Solver returns the internal amg_precond
390 virtual operator HYPRE_Solver() const
391 {
392 return amg_precond;
393 }
394
395 virtual HYPRE_PtrToParSolverFcn SetupFcn() const
396 {
397 return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSetup;
398 }
399 virtual HYPRE_PtrToParSolverFcn SolveFcn() const
400 {
401 return (HYPRE_PtrToParSolverFcn)HYPRE_BoomerAMGSolve;
402 }
403
404 virtual ~HypreBoomerAMG();
405 };
406
407 } // namespace moab
408
409 #endif // MOAB_HAVE_MPI
410
411 #endif
1.9.1.