trilinos_solver.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_TRILINOS_SOLVER_HEADER
#define OOMPH_TRILINOS_SOLVER_HEADER

#include "trilinos_preconditioners.h"
#include "iterative_linear_solver.h"

#include "AztecOO.h"

namespace oomph 
{

//=============================================================================
/// \short An Epetra_Operator class for oomph-lib preconditioners.
/// A helper class for TrilinosOomphLibPreconditioner to allow an oomph-lib
/// preconditioner (i.e. one derived from Preconditioner) to be used with 
/// a trilinos solver (TrilinosAztecOOSolver)
//=============================================================================
 class OomphLibPreconditionerEpetraOperator : public Epetra_Operator 
 {
  public:

  /// \short Constructor - takes the pointer to the oomph-lib 
  /// preconditioner and the distribution of the preconditioner
  /// Note: the oomph-lib preconditioner must be setup.
  /// If use_eptra_values is true then the epetra vector values is used
  /// within the vectors passed to the oomph-lib
  /// preconditioner. If this is true none of the vector rebuild methods can
  /// be called.
  OomphLibPreconditionerEpetraOperator(Preconditioner* preconditioner_pt,
                                       bool use_epetra_values = false)
#ifdef OOMPH_HAS_MPI
   : Operator_comm
   (preconditioner_pt->distribution_pt()->communicator_pt()->mpi_comm()), 
   Use_epetra_values(use_epetra_values)
#else
   : Operator_comm(), Use_epetra_values(use_epetra_values)
#endif
   {
    // set the ooomph-lib preconditioner
    Oomph_lib_preconditioner_pt = preconditioner_pt;
    
    // set the preconditioner label
    Preconditioner_label = "oomph-lib Preconditioner";

    // setup the Epetra_map
    Operator_map_pt = TrilinosEpetraHelpers::create_epetra_map
     (Oomph_lib_preconditioner_pt->distribution_pt());
   }

  /// \short Destructor - deletes the Epetra_map and My_global_rows vector
  /// (if MPI)
  ~OomphLibPreconditionerEpetraOperator()
   {
    delete Operator_map_pt;
    Operator_map_pt = 0;
   }

  /// Broken copy constructor
  OomphLibPreconditionerEpetraOperator
   (const OomphLibPreconditionerEpetraOperator&)
#ifdef OOMPH_HAS_MPI
   : Operator_comm(MPI_Helpers::communicator_pt()->mpi_comm()) 
#else
   : Operator_comm()
#endif
  {
   BrokenCopy::broken_copy("OomphLibPreconditionerEpetraOperator");
  }

  /// Broken assignment operator.
  void operator=(const OomphLibPreconditionerEpetraOperator&)
  {
   BrokenCopy::broken_assign("OomphLibPreconditionerEpetraOperator");
  }

  /// Broken Epetra_Operator member - SetUseTranspose
  int SetUseTranspose(bool UseTranspose)
   {
    std::ostringstream error_message;
    error_message << "SetUseTranspose() is a pure virtual Epetra_Operator "
                  << "member that is not required for a Preconditioner" 
                  << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }

  /// Broken Epetra_Operator member - Apply
  int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
   {
    std::ostringstream error_message;
    error_message << "Apply() is a pure virtual Epetra_Operator member"
                  << "that is not required for a Preconditioner" << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }


  /// \short applies the oomph-lib preconditioner. Converts the Epetra vector
  /// applys the preconditioner by calling the oomph-lib preconditioner's
  /// preconditioner_solve functionality. 
  /// NOTE : the oomph-lib preconditioner is setup prior to being passed to 
  /// this class
  int ApplyInverse(const Epetra_MultiVector &epetra_r, 
                   Epetra_MultiVector &epetra_z) const
   {
    // the oomph-lib vector for r
    DoubleVector oomph_r;

    // copy the Epetra_MultiVector r into an oomph-lib vector
    double** r_pt;
    epetra_r.ExtractView(&r_pt);
    if (Use_epetra_values)
     {
      oomph_r.set_external_values
       (Oomph_lib_preconditioner_pt->distribution_pt(),*r_pt,false);
     }
    else
     {
      oomph_r.build(Oomph_lib_preconditioner_pt->distribution_pt(),0.0);
      unsigned nrow_local = 
       Oomph_lib_preconditioner_pt->distribution_pt()->nrow_local();
      for (unsigned i = 0; i < nrow_local; i++)
       {
        oomph_r[i] = r_pt[0][i];
       }
     }

    // oomph-lib vector for Y
    DoubleVector oomph_z;
    if (Use_epetra_values)
     {
      double** z_pt;
      epetra_z.ExtractView(&z_pt);
      DoubleVector oomph_z;
      oomph_z.set_external_values(
       Oomph_lib_preconditioner_pt->distribution_pt(),
       *z_pt,false);
     }
    else
     {
      oomph_z.build(Oomph_lib_preconditioner_pt->distribution_pt(),0.0);
     }

    // apply the preconditioner
    Oomph_lib_preconditioner_pt->preconditioner_solve(oomph_r,oomph_z);

    // if not using external data copy the oomph-lib vector oomph_Y back to Y
    if (!Use_epetra_values)
     {
      unsigned nrow_local = 
       Oomph_lib_preconditioner_pt->distribution_pt()->nrow_local();
      for (unsigned i = 0; i < nrow_local; i++)
       {
        epetra_z.ReplaceMyValue(i,0,oomph_z[i]);
       }
     }

    // return 0 to indicate success
    return 0;
   }

  /// Broken Epetra_Operator member - NormInf
  double NormInf() const
   {
    std::ostringstream error_message;
    error_message << "NormInf() is a pure virtual Epetra_Operator member"
                  << "that is not required for a Preconditioner" << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }

  /// Epetra_Operator::Label - returns a string describing the operator
  const char* Label() const
   {
    return Preconditioner_label.c_str();
   }

  /// Broken Epetra_Operator member - UseTranspose
  bool UseTranspose() const
   {
    std::ostringstream error_message;
    error_message << "UseTranspose() is a pure virtual Epetra_Operator member "
                  << "that is not required for a Preconditioner" << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }

  /// Broken Epetra_Operator member - HasNormInf
  bool HasNormInf() const
   {
    std::ostringstream error_message;
    error_message << "HasNormInf() is a pure virtual Epetra_Operator member "
                  << "that is not required for a Preconditioner" << std::endl;
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }

  /// Returns the Epetra MPI_Comm object
  const Epetra_Comm& Comm() const
   {
    return Operator_comm;
   }

  /// Epetra_Operator member - OperatorDomainMap
  const Epetra_Map& OperatorDomainMap() const
   {
    return *Operator_map_pt;
   }

  /// Epetra_Operator member - OperatorRangeMap
  const Epetra_Map& OperatorRangeMap() const
   {
    return *Operator_map_pt;
   }


   private:

  /// A pointer to the oomph-lib preconditioner
  Preconditioner* Oomph_lib_preconditioner_pt;

#ifdef OOMPH_HAS_MPI
  /// An Epetra MPI Comm object 
  Epetra_MpiComm Operator_comm;
#else
  /// An Epetra Serial Comm object 
  Epetra_SerialComm Operator_comm;
#endif

  /// \short Use the epetra data within the vectors passed to the oomph-lib
  /// preconditioner. If this is true none of the vector rebuild methods can
  /// be called.
  bool Use_epetra_values;

  /// \short A pointer to an Epetra_Map object - describes distribution of the
  /// preconditioner, in this instance it is primarily used to prescribe the 
  /// distribution
  /// of the residual and solution vector
  Epetra_Map* Operator_map_pt;

  /// a label for the preconditioner ( for Epetra_Operator::Label() )
  std::string Preconditioner_label;
 };


//=============================================================================
/// \short An interface to the Trilinos AztecOO classes allowing it to be used 
/// as an Oomph-lib LinearSolver. 
/// The AztecOO solver is a Krylov Subspace solver; the solver type (either CG,
/// GMRES or BiCGStab) can be set using solver_type(). 
/// This solver can be preconditioned with Trilinos Preconditioners (derived
/// from TrilinosPreconditionerBase) or Oomph-lib preconditioners (derived
/// from Preconditioner). Preconditioners are set using preconditioner_pt().
//=============================================================================
class TrilinosAztecOOSolver : public IterativeLinearSolver
{
 public:

 /// Constructor.
 TrilinosAztecOOSolver()
 {
  
  // By default use workaround for creating of epetra matrix that 
  // respects aztecoo's ordering requirements
  Use_aztecoo_workaround_for_epetra_matrix_setup=true;
  
  // set pointer to Null
  AztecOO_solver_pt = 0;

  // initially assume not problem based solve
  Using_problem_based_solve = false;

  // if a problem based solve is performed then it should generate a 
  // serial matrix
  Assemble_serial_jacobian = false;
  
  // by default we do not use external values in the oomph-lib preconditioner
  If_oomphlib_preconditioner_use_epetra_values = false;

  // null the pts
  Problem_pt = 0;
  Epetra_matrix_pt = 0;
  Epetra_preconditioner_pt = 0;

  // set solver defaults
  Solver_type = GMRES;
  Tolerance = 10e-10;

  // don't delete the matrix
  Delete_matrix = false;
 }

 /// Destructor - delete the solver and the matrices
 ~TrilinosAztecOOSolver()
 {
  // delete 
  clean_up_memory();

  // if Problem_based solve then the oomph matrix was generated by this class
  // and should be deleted
  if (Using_problem_based_solve)
   {
    delete Oomph_matrix_pt;
    Oomph_matrix_pt = 0; 
   }
 }

 /// Broken copy constructor.
 TrilinosAztecOOSolver(const TrilinosAztecOOSolver&)
 {
  BrokenCopy::broken_copy("TrilinosAztecOOSolver");
 }

 /// Broken assignment operator.
 void operator=(const TrilinosAztecOOSolver&)
 {
  BrokenCopy::broken_assign("TrilinosAztecOOSolver");
 }

 /// \short Enable workaround for creating of epetra matrix that respects
 /// aztecoo's ordering requirements
 void enable_aztecoo_workaround_for_epetra_matrix_setup()
  {Use_aztecoo_workaround_for_epetra_matrix_setup=true;}

 /// \short Disable workaround for creating of epetra matrix that respects
 /// aztecoo's ordering requirements
 void disable_aztecoo_workaround_for_epetra_matrix_setup()
  {Use_aztecoo_workaround_for_epetra_matrix_setup=false;}

 /// \short Is workaround for creating of epetra matrix that respects
 /// aztecoo's ordering requirements enabled?
 bool is_aztecoo_workaround_for_epetra_matrix_setup_enabled()
  {return Use_aztecoo_workaround_for_epetra_matrix_setup;}

 /// Clean up method - deletes the solver, the matrices and the preconditioner
 void clean_up_memory()
  {
   // delete the solver
   delete AztecOO_solver_pt;
   AztecOO_solver_pt = 0;

   // delete the Epetra_Operator preconditioner (only if it is a wrapper to an 
   // oomphlib preconditioner in which case only the wrapper is deleted and 
   // not the actual preconditioner).
   // if the preconditioner is Trilinos preconditioner then the
   // Epetra_Operator is deleted when that preconditioner is deleted
   if (Epetra_preconditioner_pt != 0)
    {
     if (dynamic_cast<OomphLibPreconditionerEpetraOperator*>
         (Epetra_preconditioner_pt) != 0)
      {
       delete Epetra_preconditioner_pt;
       Epetra_preconditioner_pt = 0;
      }
    }

   // don't delete the preconditioner but call its clean up memory method
   if (this->preconditioner_pt() != 0)
    {
     this->preconditioner_pt()->clean_up_memory();
    }


   // delete the matrices
   // This must now happen after the preconditioner delete because the 
   // ML preconditioner (and maybe others) use the communicator from the
   // matrix in the destructor
   delete Epetra_matrix_pt;
   Epetra_matrix_pt = 0;
  }

 /// \short Function which uses problem_pt's get_jacobian(...) function to
 /// generate a linear system which is then solved. This function deletes
 /// any existing internal data and then generates a new AztecOO solver.
 void solve(Problem* const &problem_pt,DoubleVector &solution);

 /// \short Function to solve the linear system defined by matrix_pt and rhs.
 /// \b NOTE 1. The matrix has to be of type CRDoubleMatrix or 
 /// DistributedCRDoubleMatrix.
 /// \b NOTE 2. This function will delete any existing internal data and 
 /// generate a new AztecOO solver.
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &solution);

 /// \short Function to resolve a linear system using the existing solver
 /// data, allowing a solve with a new right hand side vector. This
 /// function must be used after a call to solve(...) with
 /// enable_resolve set to true.
 void resolve(const DoubleVector &rhs,
              DoubleVector &solution);

 /// \short Disable resolve function (overloads the LinearSolver
 /// disable_resolve function).
 void disable_resolve()
 {
  Enable_resolve=false;
  clean_up_memory();
 }

 /// \short Call if the matrix can be deleted
 void enable_delete_matrix() {Delete_matrix=true;}

 ///\ short Call if the matrix can not be deleted (default)
 void disable_delete_matrix() {Delete_matrix=false;}

 /// Access function to Max_iter
 unsigned& max_iter() {return Max_iter;}

 /// Acess function to Iterations
 unsigned iterations() const {return Iterations;}

 /// Access function to Tolerance
 double& tolerance() {return Tolerance;}

 /// Access function to Solver_type
 unsigned& solver_type() {return Solver_type;}

 /// Function to return Jacobian_setup_time;
 double jacobian_setup_time() {return Jacobian_setup_time;}

 /// Function to return Linear_solver_solution_time
 double linear_solver_solution_time() {return Linear_solver_solution_time;}

 /// \short Set the assembly of the serial jacobian 
 /// when performing a problem-based solve
 void enable_assemble_serial_jacobian() {Assemble_serial_jacobian=true;}

 /// Unset the assembly of the serial jacobian
 void disable_assemble_serial_jacobian() {Assemble_serial_jacobian=false;}

 /// if this solver is using an oomph-lib preconditioner then the vectors
 /// passed to preconditioner_solve(...) should be using the values in the
 /// epetra vector as external data. If the vectors are using external
 /// data then rebuild(...) methods cannot be used be used in the 
 /// preconditioner.
 //bool& if_oomphlib_preconditioner_use_epetra_values()
 // {
 //  return If_oomphlib_preconditioner_use_epetra_values;
 // }

 /// \short Enumerated list to define which AztecOO solver is used
 enum AztecOO_solver_types{CG,
                           GMRES,
                           BiCGStab};

 protected:

 /// \short Helper function performs the actual solve once the AztecOO
 /// solver is set up
 void solve_using_AztecOO(Epetra_Vector* &rhs_pt, Epetra_Vector* &soln_pt);

 /// \short Helper function for setting up the solver. Converts the oomph-lib 
 /// matrices to Epetra matrices, sets up the preconditioner, creates the 
 /// Trilinos Aztec00 solver and passes in the matrix, preconditioner and 
 /// parameters.
 void solver_setup(DoubleMatrixBase* const& matrix_pt);

 /// \short Use workaround for creating of epetra matrix that respects
 /// aztecoo's ordering requirements
 bool Use_aztecoo_workaround_for_epetra_matrix_setup;
 
 /// Stores number of iterations used
 unsigned Iterations;

 /// Pointer to the AztecOO solver
 AztecOO* AztecOO_solver_pt;

 /// Stores set up time for Jacobian
 double Jacobian_setup_time;

 /// Stores time for the solution (excludes time to set up preconditioner)
 double Linear_solver_solution_time;

 /// \short Trilinos copies matrix data from oomph-lib's own CRDoubleMatrix
 /// or DistributedCRDoubleMatrix to Trilinos's Epetra format - the Trilinos
 /// solver no longer requires the oomph-lib matrices and therefore they could
 /// be deleted to save memory. This must be requested explicitly by setting
 /// this flag to TRUE. 
 /// \b NOTE: The matrix is deleted after the preconditioner is setup.
 bool Delete_matrix;

 /// \short If true, when performing a problem based solve a serial matrix
 /// will be requested from Problem::get_jacobian(...). Defaults to true
 bool Assemble_serial_jacobian;

 /// \short Defines which solver is set up - available types are
 /// defined in AztecOO_solver_types
 unsigned Solver_type;

 /// \short Helper flag keeping track of whether we called the 
 /// linear algebra or problem-based solve function.
 bool Using_problem_based_solve;
 
 /// \short A pointer for the linear system matrix in Epetra_CrsMatrix format 
 Epetra_CrsMatrix* Epetra_matrix_pt;

 /// \short A pointer to the Epetra_Operator for the preconditioner. This is 
 /// only used if the preconditioner NOT a Trilinos preconditioner.
 Epetra_Operator* Epetra_preconditioner_pt;

 /// Oomph lib matrix pointer
 DoubleMatrixBase* Oomph_matrix_pt;

 /// \short A pointer to the underlying problem (NULL if MATRIX based solve)
 /// The problem_pt is stored here in a problem based solve for the
 /// preconditioner
 Problem* Problem_pt;

 /// if this solver is using an oomph-lib preconditioner then the vectors
 /// passed to preconditioner_solve(...) should be using the values in the
 /// epetra vector as external data. If the vectors are using external
 /// data then rebuild(...) methods cannot be used.
 bool If_oomphlib_preconditioner_use_epetra_values;
};

}
#endif