trilinos_preconditioners.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_OPERATORS_HEADER
#define OOMPH_TRILINOS_OPERATORS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// Needed in trilinos (as of gcc 4.6.* or so...)
#include <cstddef>
 
// trilinos headers
#include "ml_include.h"
#include "ml_MultiLevelPreconditioner.h"
#include "Ifpack.h"
 
// oomph-lib headers
#include "trilinos_helpers.h"
#include "preconditioner.h"
 
namespace oomph
{
  //=============================================================================
  //=============================================================================
  class TrilinosPreconditionerBase : public Preconditioner
  {
  public:
    TrilinosPreconditionerBase()
    {
      // Initialise pointers
      Epetra_preconditioner_pt = 0;
      Epetra_matrix_pt = 0;
    }
 
    virtual ~TrilinosPreconditionerBase()
    {
      clean_up_memory();
    }
 
    static double Cumulative_preconditioner_solve_time;
 
    void clean_up_memory()
    {
      // delete the Epetra preconditioner
      delete Epetra_preconditioner_pt;
      Epetra_preconditioner_pt = 0;
 
      // delete the epetra matrix
      delete Epetra_matrix_pt;
      Epetra_matrix_pt = 0;
    }
 
    TrilinosPreconditionerBase(const TrilinosPreconditionerBase&) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const TrilinosPreconditionerBase&) = delete;*/
 
    void setup();
 
    void setup(Epetra_CrsMatrix* epetra_matrix_pt);
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    Epetra_Operator*& epetra_operator_pt()
    {
      return Epetra_preconditioner_pt;
    }
 
    Epetra_Operator* epetra_operator_pt() const
    {
      return Epetra_preconditioner_pt;
    }
 
  protected:
    virtual void setup_trilinos_preconditioner(
      Epetra_CrsMatrix* epetra_matrix_pt) = 0;
 
    Epetra_Operator* Epetra_preconditioner_pt;
 
    Epetra_CrsMatrix* Epetra_matrix_pt;
  };
 
 
  //============================================================================
  //============================================================================
  class TrilinosMLPreconditioner : public TrilinosPreconditionerBase
  {
  public:
    TrilinosMLPreconditioner()
    {
      // set default values
      ML_Epetra::SetDefaults("SA", ML_parameters);
 
      // Set number of MG cycles performed in preconditioner
      ML_parameters.set("cycle applications", Default_n_cycles);
    }
 
    virtual ~TrilinosMLPreconditioner() {}
 
    TrilinosMLPreconditioner(const TrilinosMLPreconditioner&) = delete;
 
    /*void operator=(const TrilinosMLPreconditioner&) = delete;*/
 
    void set_NSSA_default_values()
    {
      ML_Epetra::SetDefaults("NSSA", ML_parameters);
    }
 
 
    void set_DD_default_values()
    {
      ML_Epetra::SetDefaults("DD", ML_parameters);
    }
 
 
    void set_DDML_default_values()
    {
      ML_Epetra::SetDefaults("DD-ML", ML_parameters);
    }
 
    void set_SA_default_values()
    {
      ML_Epetra::SetDefaults("SA", ML_parameters);
    }
 
    void set_max_levels(int max_levels)
    {
      ML_parameters.set("max levels", max_levels);
    }
 
    void set_n_cycles(int n_cycles)
    {
      ML_parameters.set("cycle applications", n_cycles);
    }
 
    void set_smoother_damping(double smoother_damping)
    {
      ML_parameters.set("smoother: damping factor", smoother_damping);
    }
 
    void set_smoother_sweeps(int smoother_sweeps)
    {
      ML_parameters.set("smoother: sweeps", smoother_sweeps);
    }
 
    void set_smoother_jacobi()
    {
      ML_parameters.set("smoother: type", "Jacobi");
    }
 
    void set_smoother_gauss_seidel()
    {
      ML_parameters.set("smoother: type", "symmetric Gauss-Seidel");
    }
 
    void set_output(int output)
    {
      ML_parameters.set("output", output);
    }
 
 
    static int Default_n_cycles;
 
  protected:
    void setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt);
 
    // Parameter list of control flags for the preconditioner
    Teuchos::ParameterList ML_parameters;
  };
 
 
  //============================================================================
  //============================================================================
  class TrilinosIFPACKPreconditioner : public TrilinosPreconditionerBase
  {
  public:
    TrilinosIFPACKPreconditioner()
    {
      // set default values
      Preconditioner_type = "ILU";
      ILU_fill_level = 0;
      ILUT_fill_level = 1.0;
      Overlap = 0;
      Absolute_threshold = 0.0;
      Relative_threshold = 1.0;
    }
 
    virtual ~TrilinosIFPACKPreconditioner() {}
 
    TrilinosIFPACKPreconditioner(const TrilinosIFPACKPreconditioner&) = delete;
 
    /*void operator=(const TrilinosIFPACKPreconditioner&) = delete;*/
 
    void set_preconditioner_ILU()
    {
      Preconditioner_type = "ILU";
    }
 
    void set_preconditioner_ILUT()
    {
      Preconditioner_type = "ILUT";
    }
 
    int& ilu_fill_level()
    {
      return ILU_fill_level;
    }
 
    double& ilut_fill_level()
    {
      return ILUT_fill_level;
    }
 
    double& absolute_threshold()
    {
      return Absolute_threshold;
    }
 
    double& relative_threshold()
    {
      return Relative_threshold;
    }
 
  protected:
    void setup_trilinos_preconditioner(Epetra_CrsMatrix* epetra_matrix_pt);
 
    string Preconditioner_type;
 
    int ILU_fill_level;
 
    double ILUT_fill_level;
 
    int Overlap;
 
    double Absolute_threshold;
 
    double Relative_threshold;
  };
} // namespace oomph
#endif