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// Copyright (C) 2006-2022 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
{
  //=============================================================================
  //=============================================================================
  class OomphLibPreconditionerEpetraOperator : public Epetra_Operator
  {
  public:
    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());
    }
 
    ~OomphLibPreconditionerEpetraOperator()
    {
      delete Operator_map_pt;
      Operator_map_pt = 0;
    }
 
    OomphLibPreconditionerEpetraOperator(
      const OomphLibPreconditionerEpetraOperator&) = delete;
 
    void operator=(const OomphLibPreconditionerEpetraOperator&) = delete;
 
    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);
    }
 
    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);
    }
 
 
    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;
    }
 
    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);
    }
 
    const char* Label() const
    {
      return Preconditioner_label.c_str();
    }
 
    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);
    }
 
    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);
    }
 
    const Epetra_Comm& Comm() const
    {
      return Operator_comm;
    }
 
    const Epetra_Map& OperatorDomainMap() const
    {
      return *Operator_map_pt;
    }
 
    const Epetra_Map& OperatorRangeMap() const
    {
      return *Operator_map_pt;
    }
 
 
  private:
    Preconditioner* Oomph_lib_preconditioner_pt;
 
#ifdef OOMPH_HAS_MPI
    Epetra_MpiComm Operator_comm;
#else
    Epetra_SerialComm Operator_comm;
#endif
 
    bool Use_epetra_values;
 
    Epetra_Map* Operator_map_pt;
 
    std::string Preconditioner_label;
  };
 
 
  //=============================================================================
  //=============================================================================
  class TrilinosAztecOOSolver : public IterativeLinearSolver
  {
  public:
    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;
    }
 
    ~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;
      }
    }
 
    TrilinosAztecOOSolver(const TrilinosAztecOOSolver&) = delete;
 
    void operator=(const TrilinosAztecOOSolver&) = delete;
 
    void enable_aztecoo_workaround_for_epetra_matrix_setup()
    {
      Use_aztecoo_workaround_for_epetra_matrix_setup = true;
    }
 
    void disable_aztecoo_workaround_for_epetra_matrix_setup()
    {
      Use_aztecoo_workaround_for_epetra_matrix_setup = false;
    }
 
    bool is_aztecoo_workaround_for_epetra_matrix_setup_enabled()
    {
      return Use_aztecoo_workaround_for_epetra_matrix_setup;
    }
 
    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;
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& solution);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution);
 
    void resolve(const DoubleVector& rhs, DoubleVector& solution);
 
    void disable_resolve()
    {
      Enable_resolve = false;
      clean_up_memory();
    }
 
    void enable_delete_matrix()
    {
      Delete_matrix = true;
    }
 
    void disable_delete_matrix()
    {
      Delete_matrix = false;
    }
 
    unsigned& max_iter()
    {
      return Max_iter;
    }
 
    unsigned iterations() const
    {
      return Iterations;
    }
 
    double& tolerance()
    {
      return Tolerance;
    }
 
    unsigned& solver_type()
    {
      return Solver_type;
    }
 
    double jacobian_setup_time()
    {
      return Jacobian_setup_time;
    }
 
    double linear_solver_solution_time()
    {
      return Linear_solver_solution_time;
    }
 
    void enable_assemble_serial_jacobian()
    {
      Assemble_serial_jacobian = true;
    }
 
    void disable_assemble_serial_jacobian()
    {
      Assemble_serial_jacobian = false;
    }
 
    // bool& if_oomphlib_preconditioner_use_epetra_values()
    // {
    //  return If_oomphlib_preconditioner_use_epetra_values;
    // }
 
    enum AztecOO_solver_types
    {
      CG,
      GMRES,
      BiCGStab
    };
 
  protected:
    void solve_using_AztecOO(Epetra_Vector*& rhs_pt, Epetra_Vector*& soln_pt);
 
    void solver_setup(DoubleMatrixBase* const& matrix_pt);
 
    bool Use_aztecoo_workaround_for_epetra_matrix_setup;
 
    unsigned Iterations;
 
    AztecOO* AztecOO_solver_pt;
 
    double Jacobian_setup_time;
 
    double Linear_solver_solution_time;
 
    bool Delete_matrix;
 
    bool Assemble_serial_jacobian;
 
    unsigned Solver_type;
 
    bool Using_problem_based_solve;
 
    Epetra_CrsMatrix* Epetra_matrix_pt;
 
    Epetra_Operator* Epetra_preconditioner_pt;
 
    DoubleMatrixBase* Oomph_matrix_pt;
 
    Problem* Problem_pt;
 
    bool If_oomphlib_preconditioner_use_epetra_values;
  };
 
} // namespace oomph
#endif