mumps_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//====================================================================
// This is the header file for the C++ wrapper functions for the
// mumps solver
 
// Include guards to prevent multiple inclusions of the header
#ifndef NEW_MUMPS_SOLVER_HEADER
#define NEW_MUMPS_SOLVER_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
#include "linear_solver.h"
#include "preconditioner.h"
 
#include <mumps_c_types.h>
#include <dmumps_c.h>
 
// offset macros for the 1-indexed (FORTRAN) arrays in the
// MUMPS interface; makes the control integers, info etc. easier to
// read w.r.t. the MUMPS documentation
#define ICNTL(i) icntl[(i)-1]
#define INFOG(i) infog[(i)-1]
#define INFO(i) info[(i)-1]
 
namespace oomph
{
  //====================================================================
  //====================================================================
  class MumpsSolver : public LinearSolver
  {
  public:
    static bool Suppress_incorrect_rhs_distribution_warning_in_resolve;
 
    MumpsSolver();
 
    MumpsSolver(const MumpsSolver& dummy) = delete;
 
    void operator=(const MumpsSolver&) = delete;
 
    ~MumpsSolver();
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    void enable_suppress_warning_about_MPI_COMM_WORLD()
    {
      Suppress_warning_about_MPI_COMM_WORLD = true;
    }
 
    void disable_suppress_warning_about_MPI_COMM_WORLD()
    {
      Suppress_warning_about_MPI_COMM_WORLD = false;
    }
 
    void enable_suppress_mumps_info_during_solve()
    {
      Suppress_mumps_info_during_solve = true;
    }
 
    void disable_suppress_mumps_info_during_solve()
    {
      Suppress_mumps_info_during_solve = false;
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result);
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    void enable_doc_stats()
    {
      Doc_stats = true;
    }
 
    void disable_doc_stats()
    {
      Doc_stats = false;
    }
 
    double jacobian_setup_time()
    {
      return Jacobian_setup_time;
    }
 
    virtual double linear_solver_solution_time()
    {
      return Solution_time;
    }
 
    void enable_suppress_solve()
    {
      Suppress_solve = true;
    }
 
    void disable_suppress_solve()
    {
      Suppress_solve = false;
    }
 
    void enable_delete_matrix_data()
    {
      Delete_matrix_data = true;
    }
 
    void disable_delete_matrix_data()
    {
      Delete_matrix_data = false;
    }
 
    void factorise(DoubleMatrixBase* const& matrix_pt);
 
    void backsub(const DoubleVector& rhs, DoubleVector& result);
 
    void clean_up_memory();
 
    static int Default_workspace_scaling_factor;
 
    void declare_jacobian_is_unsymmetric()
    {
      Jacobian_symmetry_flag = Unsymmetric;
    }
 
    void declare_jacobian_is_symmetric()
    {
      Jacobian_symmetry_flag = Symmetric;
    }
 
    void declare_jacobian_is_symmetric_positive_definite()
    {
      Jacobian_symmetry_flag = Symmetric_positive_definite;
    }
 
    // tell MUMPS to use the PORD package for the ordering
    void use_pord_ordering()
    {
      Jacobian_ordering_flag = Pord_ordering;
    }
 
    // tell MUMPS to use the METIS package for the ordering
    void use_metis_ordering()
    {
      Jacobian_ordering_flag = Metis_ordering;
    }
 
    // tell MUMPS to use the SCOTCH package for the ordering
    void use_scotch_ordering()
    {
      Jacobian_ordering_flag = Scotch_ordering;
    }
 
  private:
    void initialise_mumps();
 
    void shutdown_mumps();
 
    double Jacobian_setup_time;
 
    double Solution_time;
 
    bool Suppress_solve;
 
    bool Doc_stats;
 
    bool Suppress_warning_about_MPI_COMM_WORLD;
 
    bool Suppress_mumps_info_during_solve;
 
    bool Mumps_is_initialised;
 
    // Work space scaling factor
    unsigned Workspace_scaling_factor;
 
    bool Delete_matrix_data;
 
    Vector<int> Irn_loc;
 
    // Vector for column numbers
    Vector<int> Jcn_loc;
 
    // Vector for entries
    Vector<double> A_loc;
 
    DMUMPS_STRUC_C* Mumps_struc_pt;
 
    enum MumpsJacobianSymmetryFlags
    {
      Unsymmetric = 0,
      Symmetric_positive_definite = 1,
      Symmetric = 2
    };
 
    enum MumpsJacobianOrderingFlags
    {
      Scotch_ordering = 3,
      Pord_ordering = 4,
      Metis_ordering = 5
    };
 
    unsigned Jacobian_symmetry_flag;
 
    unsigned Jacobian_ordering_flag;
  };
 
 
 
 
  //====================================================================
  //====================================================================
  class NewMumpsPreconditioner : public Preconditioner
  {
  public:
    NewMumpsPreconditioner() {}
 
    ~NewMumpsPreconditioner() {}
 
    NewMumpsPreconditioner(const NewMumpsPreconditioner&) = delete;
 
    void operator=(const NewMumpsPreconditioner&) = delete;
 
    void setup()
    {
      oomph_info << "Setting up Mumps (exact) preconditioner" << std::endl;
 
      DistributableLinearAlgebraObject* dist_matrix_pt =
        dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt());
      if (dist_matrix_pt != 0)
      {
        LinearAlgebraDistribution dist(dist_matrix_pt->distribution_pt());
        this->build_distribution(dist);
        Solver.factorise(matrix_pt());
      }
      else
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "NewMumpsPreconditioner can only be applied to matrices derived \n"
          << "DistributableLinearAlgebraObject.\n";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      Solver.resolve(r, z);
    }
 
 
    void clean_up_memory()
    {
      Solver.clean_up_memory();
    }
 
 
    void enable_doc_time()
    {
      Solver.enable_doc_time();
    }
 
    void disable_doc_time()
    {
      Solver.disable_doc_time();
    }
 
  private:
    MumpsSolver Solver;
  };
 
} // namespace oomph
 
#endif