trilinos_solver.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_TRILINOS_SOLVER_HEADER
27#define OOMPH_TRILINOS_SOLVER_HEADER
28
31
32#include "AztecOO.h"
33
34namespace oomph
35{
36 //=============================================================================
37 /// An Epetra_Operator class for oomph-lib preconditioners.
38 /// A helper class for TrilinosOomphLibPreconditioner to allow an oomph-lib
39 /// preconditioner (i.e. one derived from Preconditioner) to be used with
40 /// a trilinos solver (TrilinosAztecOOSolver)
41 //=============================================================================
42 class OomphLibPreconditionerEpetraOperator : public Epetra_Operator
43 {
44 public:
45 /// Constructor - takes the pointer to the oomph-lib
46 /// preconditioner and the distribution of the preconditioner
47 /// Note: the oomph-lib preconditioner must be setup.
48 /// If use_eptra_values is true then the epetra vector values is used
49 /// within the vectors passed to the oomph-lib
50 /// preconditioner. If this is true none of the vector rebuild methods can
51 /// be called.
53 bool use_epetra_values = false)
54#ifdef OOMPH_HAS_MPI
56 preconditioner_pt->distribution_pt()->communicator_pt()->mpi_comm()),
57 Use_epetra_values(use_epetra_values)
58#else
59 : Operator_comm(), Use_epetra_values(use_epetra_values)
60#endif
61 {
62 // set the ooomph-lib preconditioner
63 Oomph_lib_preconditioner_pt = preconditioner_pt;
64
65 // set the preconditioner label
66 Preconditioner_label = "oomph-lib Preconditioner";
67
68 // setup the Epetra_map
71 }
72
73 /// Destructor - deletes the Epetra_map and My_global_rows vector
74 /// (if MPI)
76 {
77 delete Operator_map_pt;
79 }
80
81 /// Broken copy constructor
84
85 /// Broken assignment operator.
87
88 /// Broken Epetra_Operator member - SetUseTranspose
90 {
91 std::ostringstream error_message;
92 error_message << "SetUseTranspose() is a pure virtual Epetra_Operator "
93 << "member that is not required for a Preconditioner"
94 << std::endl;
95 throw OomphLibError(
96 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
97 }
98
99 /// Broken Epetra_Operator member - Apply
100 int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
101 {
102 std::ostringstream error_message;
103 error_message << "Apply() is a pure virtual Epetra_Operator member"
104 << "that is not required for a Preconditioner" << std::endl;
105 throw OomphLibError(
106 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
107 }
108
109
110 /// applies the oomph-lib preconditioner. Converts the Epetra vector
111 /// applys the preconditioner by calling the oomph-lib preconditioner's
112 /// preconditioner_solve functionality.
113 /// NOTE : the oomph-lib preconditioner is setup prior to being passed to
114 /// this class
115 int ApplyInverse(const Epetra_MultiVector& epetra_r,
116 Epetra_MultiVector& epetra_z) const
117 {
118 // the oomph-lib vector for r
119 DoubleVector oomph_r;
120
121 // copy the Epetra_MultiVector r into an oomph-lib vector
122 double** r_pt;
123 epetra_r.ExtractView(&r_pt);
125 {
126 oomph_r.set_external_values(
128 }
129 else
130 {
132 unsigned nrow_local =
134 for (unsigned i = 0; i < nrow_local; i++)
135 {
136 oomph_r[i] = r_pt[0][i];
137 }
138 }
139
140 // oomph-lib vector for Y
141 DoubleVector oomph_z;
143 {
144 double** z_pt;
145 epetra_z.ExtractView(&z_pt);
146 DoubleVector oomph_z;
147 oomph_z.set_external_values(
149 }
150 else
151 {
153 }
154
155 // apply the preconditioner
157
158 // if not using external data copy the oomph-lib vector oomph_Y back to Y
160 {
161 unsigned nrow_local =
163 for (unsigned i = 0; i < nrow_local; i++)
164 {
165 epetra_z.ReplaceMyValue(i, 0, oomph_z[i]);
166 }
167 }
168
169 // return 0 to indicate success
170 return 0;
171 }
172
173 /// Broken Epetra_Operator member - NormInf
174 double NormInf() const
175 {
176 std::ostringstream error_message;
177 error_message << "NormInf() is a pure virtual Epetra_Operator member"
178 << "that is not required for a Preconditioner" << std::endl;
179 throw OomphLibError(
180 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
181 }
182
183 /// Epetra_Operator::Label - returns a string describing the operator
184 const char* Label() const
185 {
186 return Preconditioner_label.c_str();
187 }
188
189 /// Broken Epetra_Operator member - UseTranspose
190 bool UseTranspose() const
191 {
192 std::ostringstream error_message;
193 error_message
194 << "UseTranspose() is a pure virtual Epetra_Operator member "
195 << "that is not required for a Preconditioner" << std::endl;
196 throw OomphLibError(
197 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
198 }
199
200 /// Broken Epetra_Operator member - HasNormInf
201 bool HasNormInf() const
202 {
203 std::ostringstream error_message;
204 error_message << "HasNormInf() is a pure virtual Epetra_Operator member "
205 << "that is not required for a Preconditioner" << std::endl;
206 throw OomphLibError(
207 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
208 }
209
210 /// Returns the Epetra MPI_Comm object
211 const Epetra_Comm& Comm() const
212 {
213 return Operator_comm;
214 }
215
216 /// Epetra_Operator member - OperatorDomainMap
217 const Epetra_Map& OperatorDomainMap() const
218 {
219 return *Operator_map_pt;
220 }
221
222 /// Epetra_Operator member - OperatorRangeMap
223 const Epetra_Map& OperatorRangeMap() const
224 {
225 return *Operator_map_pt;
226 }
227
228
229 private:
230 /// A pointer to the oomph-lib preconditioner
232
233#ifdef OOMPH_HAS_MPI
234 /// An Epetra MPI Comm object
235 Epetra_MpiComm Operator_comm;
236#else
237 /// An Epetra Serial Comm object
238 Epetra_SerialComm Operator_comm;
239#endif
240
241 /// Use the epetra data within the vectors passed to the oomph-lib
242 /// preconditioner. If this is true none of the vector rebuild methods can
243 /// be called.
245
246 /// A pointer to an Epetra_Map object - describes distribution of the
247 /// preconditioner, in this instance it is primarily used to prescribe the
248 /// distribution
249 /// of the residual and solution vector
250 Epetra_Map* Operator_map_pt;
251
252 /// a label for the preconditioner ( for Epetra_Operator::Label() )
254 };
255
256
257 //=============================================================================
258 /// An interface to the Trilinos AztecOO classes allowing it to be used
259 /// as an Oomph-lib LinearSolver.
260 /// The AztecOO solver is a Krylov Subspace solver; the solver type (either
261 /// CG, GMRES or BiCGStab) can be set using solver_type(). This solver can be
262 /// preconditioned with Trilinos Preconditioners (derived from
263 /// TrilinosPreconditionerBase) or Oomph-lib preconditioners (derived from
264 /// Preconditioner). Preconditioners are set using preconditioner_pt().
265 //=============================================================================
267 {
268 public:
269 /// Constructor.
271 {
272 // By default use workaround for creating of epetra matrix that
273 // respects aztecoo's ordering requirements
275
276 // set pointer to Null
278
279 // initially assume not problem based solve
281
282 // if a problem based solve is performed then it should generate a
283 // serial matrix
285
286 // by default we do not use external values in the oomph-lib
287 // preconditioner
289
290 // null the pts
291 Problem_pt = 0;
294
295 // set solver defaults
297 Tolerance = 10e-10;
298
299 // don't delete the matrix
300 Delete_matrix = false;
301 }
302
303 /// Destructor - delete the solver and the matrices
305 {
306 // delete
308
309 // if Problem_based solve then the oomph matrix was generated by this
310 // class and should be deleted
312 {
313 delete Oomph_matrix_pt;
314 Oomph_matrix_pt = 0;
315 }
316 }
317
318 /// Broken copy constructor.
320
321 /// Broken assignment operator.
322 void operator=(const TrilinosAztecOOSolver&) = delete;
323
324 /// Enable workaround for creating of epetra matrix that respects
325 /// aztecoo's ordering requirements
327 {
329 }
330
331 /// Disable workaround for creating of epetra matrix that respects
332 /// aztecoo's ordering requirements
334 {
336 }
337
338 /// Is workaround for creating of epetra matrix that respects
339 /// aztecoo's ordering requirements enabled?
341 {
343 }
344
345 /// Clean up method - deletes the solver, the matrices and the
346 /// preconditioner
348 {
349 // delete the solver
350 delete AztecOO_solver_pt;
352
353 // delete the Epetra_Operator preconditioner (only if it is a wrapper to
354 // an oomphlib preconditioner in which case only the wrapper is deleted
355 // and not the actual preconditioner). if the preconditioner is Trilinos
356 // preconditioner then the Epetra_Operator is deleted when that
357 // preconditioner is deleted
359 {
360 if (dynamic_cast<OomphLibPreconditionerEpetraOperator*>(
362 {
365 }
366 }
367
368 // don't delete the preconditioner but call its clean up memory method
369 if (this->preconditioner_pt() != 0)
370 {
372 }
373
374
375 // delete the matrices
376 // This must now happen after the preconditioner delete because the
377 // ML preconditioner (and maybe others) use the communicator from the
378 // matrix in the destructor
379 delete Epetra_matrix_pt;
381 }
382
383 /// Function which uses problem_pt's get_jacobian(...) function to
384 /// generate a linear system which is then solved. This function deletes
385 /// any existing internal data and then generates a new AztecOO solver.
386 void solve(Problem* const& problem_pt, DoubleVector& solution);
387
388 /// Function to solve the linear system defined by matrix_pt and rhs.
389 /// \b NOTE 1. The matrix has to be of type CRDoubleMatrix or
390 /// DistributedCRDoubleMatrix.
391 /// \b NOTE 2. This function will delete any existing internal data and
392 /// generate a new AztecOO solver.
393 void solve(DoubleMatrixBase* const& matrix_pt,
394 const DoubleVector& rhs,
395 DoubleVector& solution);
396
397 /// Function to resolve a linear system using the existing solver
398 /// data, allowing a solve with a new right hand side vector. This
399 /// function must be used after a call to solve(...) with
400 /// enable_resolve set to true.
401 void resolve(const DoubleVector& rhs, DoubleVector& solution);
402
403 /// Disable resolve function (overloads the LinearSolver
404 /// disable_resolve function).
406 {
407 Enable_resolve = false;
409 }
410
411 /// Call if the matrix can be deleted
413 {
414 Delete_matrix = true;
415 }
416
417 /// Call if the matrix can not be deleted (default)
419 {
420 Delete_matrix = false;
421 }
422
423 /// Access function to Max_iter
424 unsigned& max_iter()
425 {
426 return Max_iter;
427 }
428
429 /// Acess function to Iterations
430 unsigned iterations() const
431 {
432 return Iterations;
433 }
434
435 /// Access function to Tolerance
436 double& tolerance()
437 {
438 return Tolerance;
439 }
440
441 /// Access function to Solver_type
442 unsigned& solver_type()
443 {
444 return Solver_type;
445 }
446
447 /// Function to return Jacobian_setup_time;
449 {
450 return Jacobian_setup_time;
451 }
452
453 /// Function to return Linear_solver_solution_time
455 {
457 }
458
459 /// Set the assembly of the serial jacobian
460 /// when performing a problem-based solve
462 {
464 }
465
466 /// Unset the assembly of the serial jacobian
468 {
470 }
471
472 /// if this solver is using an oomph-lib preconditioner then the vectors
473 /// passed to preconditioner_solve(...) should be using the values in the
474 /// epetra vector as external data. If the vectors are using external
475 /// data then rebuild(...) methods cannot be used be used in the
476 /// preconditioner.
477 // bool& if_oomphlib_preconditioner_use_epetra_values()
478 // {
479 // return If_oomphlib_preconditioner_use_epetra_values;
480 // }
481
482 /// Enumerated list to define which AztecOO solver is used
484 {
488 };
489
490 protected:
491 /// Helper function performs the actual solve once the AztecOO
492 /// solver is set up
493 void solve_using_AztecOO(Epetra_Vector*& rhs_pt, Epetra_Vector*& soln_pt);
494
495 /// Helper function for setting up the solver. Converts the oomph-lib
496 /// matrices to Epetra matrices, sets up the preconditioner, creates the
497 /// Trilinos Aztec00 solver and passes in the matrix, preconditioner and
498 /// parameters.
499 void solver_setup(DoubleMatrixBase* const& matrix_pt);
500
501 /// Use workaround for creating of epetra matrix that respects
502 /// aztecoo's ordering requirements
504
505 /// Stores number of iterations used
506 unsigned Iterations;
507
508 /// Pointer to the AztecOO solver
510
511 /// Stores set up time for Jacobian
513
514 /// Stores time for the solution (excludes time to set up preconditioner)
516
517 /// Trilinos copies matrix data from oomph-lib's own CRDoubleMatrix
518 /// or DistributedCRDoubleMatrix to Trilinos's Epetra format - the Trilinos
519 /// solver no longer requires the oomph-lib matrices and therefore they
520 /// could be deleted to save memory. This must be requested explicitly by
521 /// setting this flag to TRUE. \b NOTE: The matrix is deleted after the
522 /// preconditioner is setup.
524
525 /// If true, when performing a problem based solve a serial matrix
526 /// will be requested from Problem::get_jacobian(...). Defaults to true
528
529 /// Defines which solver is set up - available types are
530 /// defined in AztecOO_solver_types
531 unsigned Solver_type;
532
533 /// Helper flag keeping track of whether we called the
534 /// linear algebra or problem-based solve function.
536
537 /// A pointer for the linear system matrix in Epetra_CrsMatrix format
538 Epetra_CrsMatrix* Epetra_matrix_pt;
539
540 /// A pointer to the Epetra_Operator for the preconditioner. This is
541 /// only used if the preconditioner NOT a Trilinos preconditioner.
542 Epetra_Operator* Epetra_preconditioner_pt;
543
544 /// Oomph lib matrix pointer
546
547 /// A pointer to the underlying problem (NULL if MATRIX based solve)
548 /// The problem_pt is stored here in a problem based solve for the
549 /// preconditioner
551
552 /// if this solver is using an oomph-lib preconditioner then the vectors
553 /// passed to preconditioner_solve(...) should be using the values in the
554 /// epetra vector as external data. If the vectors are using external
555 /// data then rebuild(...) methods cannot be used.
557 };
558
559} // namespace oomph
560#endif
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Abstract base class for matrices of doubles – adds abstract interfaces for solving,...
Definition: matrices.h:261
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void set_external_values(const LinearAlgebraDistribution *const &dist_pt, double *external_values, bool delete_external_values)
Allows are external data to be used by this vector. WARNING: The size of the external data must corre...
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
double Tolerance
Convergence tolerance.
unsigned Max_iter
Maximum number of iterations.
Preconditioner *& preconditioner_pt()
Access function to preconditioner.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
bool Enable_resolve
Boolean that indicates whether the matrix (or its factors, in the case of direct solver) should be st...
Definition: linear_solver.h:73
An OomphLibError object which should be thrown when an run-time error is encountered....
An Epetra_Operator class for oomph-lib preconditioners. A helper class for TrilinosOomphLibPreconditi...
Epetra_SerialComm Operator_comm
An Epetra Serial Comm object.
Preconditioner * Oomph_lib_preconditioner_pt
A pointer to the oomph-lib preconditioner.
Epetra_MpiComm Operator_comm
An Epetra MPI Comm object.
OomphLibPreconditionerEpetraOperator(const OomphLibPreconditionerEpetraOperator &)=delete
Broken copy constructor.
bool Use_epetra_values
Use the epetra data within the vectors passed to the oomph-lib preconditioner. If this is true none o...
bool HasNormInf() const
Broken Epetra_Operator member - HasNormInf.
int ApplyInverse(const Epetra_MultiVector &epetra_r, Epetra_MultiVector &epetra_z) const
applies the oomph-lib preconditioner. Converts the Epetra vector applys the preconditioner by calling...
const char * Label() const
Epetra_Operator::Label - returns a string describing the operator.
void operator=(const OomphLibPreconditionerEpetraOperator &)=delete
Broken assignment operator.
std::string Preconditioner_label
a label for the preconditioner ( for Epetra_Operator::Label() )
bool UseTranspose() const
Broken Epetra_Operator member - UseTranspose.
const Epetra_Comm & Comm() const
Returns the Epetra MPI_Comm object.
double NormInf() const
Broken Epetra_Operator member - NormInf.
OomphLibPreconditionerEpetraOperator(Preconditioner *preconditioner_pt, bool use_epetra_values=false)
Constructor - takes the pointer to the oomph-lib preconditioner and the distribution of the precondit...
const Epetra_Map & OperatorDomainMap() const
Epetra_Operator member - OperatorDomainMap.
int SetUseTranspose(bool UseTranspose)
Broken Epetra_Operator member - SetUseTranspose.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Broken Epetra_Operator member - Apply.
const Epetra_Map & OperatorRangeMap() const
Epetra_Operator member - OperatorRangeMap.
Epetra_Map * Operator_map_pt
A pointer to an Epetra_Map object - describes distribution of the preconditioner, in this instance it...
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
virtual void clean_up_memory()
Clean up memory (empty). Generic interface function.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: problem.h:151
An interface to the Trilinos AztecOO classes allowing it to be used as an Oomph-lib LinearSolver....
bool Using_problem_based_solve
Helper flag keeping track of whether we called the linear algebra or problem-based solve function.
void resolve(const DoubleVector &rhs, DoubleVector &solution)
Function to resolve a linear system using the existing solver data, allowing a solve with a new right...
bool If_oomphlib_preconditioner_use_epetra_values
if this solver is using an oomph-lib preconditioner then the vectors passed to preconditioner_solve(....
double & tolerance()
Access function to Tolerance.
void solver_setup(DoubleMatrixBase *const &matrix_pt)
Helper function for setting up the solver. Converts the oomph-lib matrices to Epetra matrices,...
void disable_aztecoo_workaround_for_epetra_matrix_setup()
Disable workaround for creating of epetra matrix that respects aztecoo's ordering requirements.
unsigned Iterations
Stores number of iterations used.
Epetra_CrsMatrix * Epetra_matrix_pt
A pointer for the linear system matrix in Epetra_CrsMatrix format.
void clean_up_memory()
Clean up method - deletes the solver, the matrices and the preconditioner.
double jacobian_setup_time()
Function to return Jacobian_setup_time;.
bool Assemble_serial_jacobian
If true, when performing a problem based solve a serial matrix will be requested from Problem::get_ja...
void solve_using_AztecOO(Epetra_Vector *&rhs_pt, Epetra_Vector *&soln_pt)
Helper function performs the actual solve once the AztecOO solver is set up.
bool Delete_matrix
Trilinos copies matrix data from oomph-lib's own CRDoubleMatrix or DistributedCRDoubleMatrix to Trili...
void disable_delete_matrix()
Call if the matrix can not be deleted (default)
DoubleMatrixBase * Oomph_matrix_pt
Oomph lib matrix pointer.
Epetra_Operator * Epetra_preconditioner_pt
A pointer to the Epetra_Operator for the preconditioner. This is only used if the preconditioner NOT ...
~TrilinosAztecOOSolver()
Destructor - delete the solver and the matrices.
unsigned & solver_type()
Access function to Solver_type.
void disable_resolve()
Disable resolve function (overloads the LinearSolver disable_resolve function).
void solve(Problem *const &problem_pt, DoubleVector &solution)
Function which uses problem_pt's get_jacobian(...) function to generate a linear system which is then...
void enable_assemble_serial_jacobian()
Set the assembly of the serial jacobian when performing a problem-based solve.
double linear_solver_solution_time()
Function to return Linear_solver_solution_time.
AztecOO_solver_types
if this solver is using an oomph-lib preconditioner then the vectors passed to preconditioner_solve(....
unsigned & max_iter()
Access function to Max_iter.
void disable_assemble_serial_jacobian()
Unset the assembly of the serial jacobian.
TrilinosAztecOOSolver(const TrilinosAztecOOSolver &)=delete
Broken copy constructor.
void enable_aztecoo_workaround_for_epetra_matrix_setup()
Enable workaround for creating of epetra matrix that respects aztecoo's ordering requirements.
unsigned Solver_type
Defines which solver is set up - available types are defined in AztecOO_solver_types.
unsigned iterations() const
Acess function to Iterations.
double Linear_solver_solution_time
Stores time for the solution (excludes time to set up preconditioner)
Problem * Problem_pt
A pointer to the underlying problem (NULL if MATRIX based solve) The problem_pt is stored here in a p...
bool is_aztecoo_workaround_for_epetra_matrix_setup_enabled()
Is workaround for creating of epetra matrix that respects aztecoo's ordering requirements enabled?
void enable_delete_matrix()
Call if the matrix can be deleted.
AztecOO * AztecOO_solver_pt
Pointer to the AztecOO solver.
void operator=(const TrilinosAztecOOSolver &)=delete
Broken assignment operator.
bool Use_aztecoo_workaround_for_epetra_matrix_setup
Use workaround for creating of epetra matrix that respects aztecoo's ordering requirements.
double Jacobian_setup_time
Stores set up time for Jacobian.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
Epetra_Map * create_epetra_map(const LinearAlgebraDistribution *const dist)
create an Epetra_Map corresponding to the LinearAlgebraDistribution
//////////////////////////////////////////////////////////////////// ////////////////////////////////...