general_purpose_block_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//====================================================================
 
// Include guards
#ifndef OOMPH_GENERAL_BLOCK_PRECONDITIONERS
#define OOMPH_GENERAL_BLOCK_PRECONDITIONERS
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// c++ include
// #include<list>
 
// oomph-lib includes
#include "matrices.h"
// #include "mesh.h"
// #include "problem.h"
#include "block_preconditioner.h"
#include "SuperLU_preconditioner.h"
#include "preconditioner_array.h"
#include "matrix_vector_product.h"
 
 
namespace oomph
{
  namespace PreconditionerCreationFunctions
  {
    inline Preconditioner* create_super_lu_preconditioner()
    {
      return new SuperLUPreconditioner;
    }
  } // namespace PreconditionerCreationFunctions
 
 
  //============================================================================
  //============================================================================
  template<typename MATRIX>
  class GeneralPurposeBlockPreconditioner : public BlockPreconditioner<MATRIX>
  {
  public:
    typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
 
    GeneralPurposeBlockPreconditioner()
      : BlockPreconditioner<MATRIX>(),
        Subsidiary_preconditioner_creation_function_pt(
          &PreconditionerCreationFunctions::create_super_lu_preconditioner)
    {
      // Make sure that the Gp_mesh_pt container is size zero.
      Gp_mesh_pt.resize(0);
    }
 
    virtual ~GeneralPurposeBlockPreconditioner()
    {
      this->clean_up_memory();
 
      for (unsigned j = 0, nj = Subsidiary_preconditioner_pt.size(); j < nj;
           j++)
      {
        delete Subsidiary_preconditioner_pt[j];
      }
    }
 
    virtual void clean_up_memory()
    {
      // Call clean up in any subsidiary precondtioners that are set.
      for (unsigned j = 0, nj = Subsidiary_preconditioner_pt.size(); j < nj;
           j++)
      {
        if (Subsidiary_preconditioner_pt[j] != 0)
        {
          Subsidiary_preconditioner_pt[j]->clean_up_memory();
        }
      }
 
      // Clean up the block preconditioner base class stuff
      this->clear_block_preconditioner_base();
    }
 
    GeneralPurposeBlockPreconditioner(
      const GeneralPurposeBlockPreconditioner&) = delete;
 
    void operator=(const GeneralPurposeBlockPreconditioner&) = delete;
 
    void set_subsidiary_preconditioner_function(
      SubsidiaryPreconditionerFctPt sub_prec_fn)
    {
      Subsidiary_preconditioner_creation_function_pt = sub_prec_fn;
    }
 
    void reset_subsidiary_preconditioner_function_to_default()
    {
      Subsidiary_preconditioner_creation_function_pt =
        &PreconditionerCreationFunctions::create_super_lu_preconditioner;
    }
    void set_subsidiary_preconditioner_pt(Preconditioner* prec,
                                          const unsigned& i)
    {
      // If the vector is currently too small to hold that many
      // preconditioners then expand it and fill with nulls.
      if (Subsidiary_preconditioner_pt.size() < i + 1)
      {
        Subsidiary_preconditioner_pt.resize(i + 1, 0);
      }
      // Note: the size of the vector is checked by
      // fill_in_subsidiary_preconditioners(..)  when we know what size it
      // should be.
 
      // I'm assuming that the number of preconditioners is always "small"
      // compared to Jacobian size, so a resize doesn't waste much time.
 
      // Put the pointer in the vector
      Subsidiary_preconditioner_pt[i] = prec;
    }
 
    Preconditioner* subsidiary_preconditioner_pt(const unsigned& i) const
    {
      return Subsidiary_preconditioner_pt[i];
    }
 
    void set_dof_to_block_map(Vector<unsigned>& dof_to_block_map)
    {
      Dof_to_block_map = dof_to_block_map;
    }
 
    void add_mesh(const Mesh* mesh_pt,
                  const bool& allow_multiple_element_type_in_mesh = false)
    {
#ifdef PARANOID
      // Check that the mesh pointer is not null.
      if (mesh_pt == 0)
      {
        std::ostringstream err_msg;
        err_msg << "The mesh_pt is null, please point it to a mesh.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // Push back the mesh pointer and the boolean in a pair.
      Gp_mesh_pt.push_back(
        std::make_pair(mesh_pt, allow_multiple_element_type_in_mesh));
    }
 
    unsigned gp_nmesh()
    {
      return Gp_mesh_pt.size();
    }
 
  protected:
    void gp_preconditioner_set_all_meshes()
    {
      const unsigned nmesh = gp_nmesh();
#ifdef PARANOID
      if (nmesh == 0)
      {
        std::ostringstream err_msg;
        err_msg << "There are no meshes set.\n"
                << "Have you remembered to call add_mesh(...)?\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      this->set_nmesh(nmesh);
      for (unsigned mesh_i = 0; mesh_i < nmesh; mesh_i++)
      {
        this->set_mesh(
          mesh_i, Gp_mesh_pt[mesh_i].first, Gp_mesh_pt[mesh_i].second);
      }
    }
 
    void gp_preconditioner_block_setup()
    {
      if (Dof_to_block_map.size() > 0)
      {
        BlockPreconditioner<MATRIX>::block_setup(Dof_to_block_map);
      }
      else
      {
        BlockPreconditioner<MATRIX>::block_setup();
      }
    }
 
    void fill_in_subsidiary_preconditioners(const unsigned& nprec_needed)
    {
      // If it's empty then fill it in with null pointers.
      if (Subsidiary_preconditioner_pt.empty())
      {
        Subsidiary_preconditioner_pt.assign(nprec_needed, 0);
      }
      else
      {
        // Otherwise check we have the right number of them
#ifdef PARANOID
        if (Subsidiary_preconditioner_pt.size() != nprec_needed)
        {
          using namespace StringConversion;
          std::string error_msg = "Wrong number of precondtioners in";
          error_msg += "Subsidiary_preconditioner_pt, should have ";
          error_msg += to_string(nprec_needed) + " but we actually have ";
          error_msg += to_string(Subsidiary_preconditioner_pt.size());
          throw OomphLibError(
            error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
      }
 
 
      // Now replace any null pointers with new preconditioners
      for (unsigned j = 0, nj = Subsidiary_preconditioner_pt.size(); j < nj;
           j++)
      {
        if (Subsidiary_preconditioner_pt[j] == 0)
        {
          Subsidiary_preconditioner_pt[j] =
            (*Subsidiary_preconditioner_creation_function_pt)();
        }
      }
    }
 
    Vector<Preconditioner*> Subsidiary_preconditioner_pt;
 
    SubsidiaryPreconditionerFctPt
      Subsidiary_preconditioner_creation_function_pt;
 
  private:
    Vector<unsigned> Dof_to_block_map;
 
    Vector<std::pair<const Mesh*, bool>> Gp_mesh_pt;
  };
 
 
  //=============================================================================
  //=============================================================================
  template<typename MATRIX>
  class BlockDiagonalPreconditioner
    : public GeneralPurposeBlockPreconditioner<MATRIX>
  {
  public:
    BlockDiagonalPreconditioner() : GeneralPurposeBlockPreconditioner<MATRIX>()
    {
      // by default we do not use two level parallelism
      Use_two_level_parallelisation = false;
 
      // null the Preconditioner array pt
      Preconditioner_array_pt = 0;
 
      // Don't doc by default
      Doc_time_during_preconditioner_solve = false;
    }
 
    virtual ~BlockDiagonalPreconditioner()
    {
      this->clean_up_memory();
    }
 
    virtual void clean_up_memory()
    {
      if (Use_two_level_parallelisation)
      {
        delete Preconditioner_array_pt;
        Preconditioner_array_pt = 0;
      }
 
      // Clean up the base class too
      GeneralPurposeBlockPreconditioner<MATRIX>::clean_up_memory();
    }
 
    BlockDiagonalPreconditioner(const BlockDiagonalPreconditioner&) = delete;
 
    void operator=(const BlockDiagonalPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    virtual void setup();
 
    void enable_two_level_parallelisation()
    {
#ifndef OOMPH_HAS_MPI
      throw OomphLibError("Cannot do any parallelism since we don't have MPI.",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
#endif
      Use_two_level_parallelisation = true;
    }
 
    void disable_two_level_parallelisation()
    {
      Use_two_level_parallelisation = false;
    }
 
    void enable_doc_time_during_preconditioner_solve()
    {
      Doc_time_during_preconditioner_solve = true;
    }
 
    void disable_doc_time_during_preconditioner_solve()
    {
      Doc_time_during_preconditioner_solve = false;
    }
 
    void fill_in_subsidiary_preconditioners(const unsigned& nprec_needed)
    {
#ifdef PARANOID
      if ((Use_two_level_parallelisation) &&
          !this->Subsidiary_preconditioner_pt.empty())
      {
        std::string err_msg =
          "Two level parallelism diagonal block preconditioners cannot have";
        err_msg +=
          " any preset preconditioners (due to weird memory management";
        err_msg += " in the PreconditionerArray, you could try fixing it).";
        throw OomphLibError(
          err_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Now call the real function
      GeneralPurposeBlockPreconditioner<
        MATRIX>::fill_in_subsidiary_preconditioners(nprec_needed);
    }
 
  protected:
    virtual unsigned get_other_diag_ds(const unsigned& i,
                                       const unsigned& nblock) const
    {
      return i;
    }
 
 
  private:
    PreconditionerArray* Preconditioner_array_pt;
 
    bool Use_two_level_parallelisation;
 
    bool Doc_time_during_preconditioner_solve;
  };
 
 
 
 
  //=============================================================================
  //=============================================================================
  template<typename MATRIX>
  class BlockTriangularPreconditioner
    : public GeneralPurposeBlockPreconditioner<MATRIX>
  {
  public:
    BlockTriangularPreconditioner()
      : GeneralPurposeBlockPreconditioner<MATRIX>()
    {
      // default to upper triangular
      Upper_triangular = true;
    }
 
    virtual ~BlockTriangularPreconditioner()
    {
      this->clean_up_memory();
    }
 
    virtual void clean_up_memory()
    {
      // Delete anything in Off_diagonal_matrix_vector_products
      for (unsigned i = 0, ni = Off_diagonal_matrix_vector_products.nrow();
           i < ni;
           i++)
      {
        for (unsigned j = 0, nj = Off_diagonal_matrix_vector_products.ncol();
             j < nj;
             j++)
        {
          delete Off_diagonal_matrix_vector_products(i, j);
          Off_diagonal_matrix_vector_products(i, j) = 0;
        }
      }
 
      // Clean up the base class too
      GeneralPurposeBlockPreconditioner<MATRIX>::clean_up_memory();
    }
 
    BlockTriangularPreconditioner(const BlockTriangularPreconditioner&) =
      delete;
 
    void operator=(const BlockTriangularPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
    void upper_triangular()
    {
      Upper_triangular = true;
    }
 
    void lower_triangular()
    {
      Upper_triangular = false;
    }
 
  private:
    DenseMatrix<MatrixVectorProduct*> Off_diagonal_matrix_vector_products;
 
    bool Upper_triangular;
  };
 
 
 
 
  //=============================================================================
  //=============================================================================
  template<typename MATRIX>
  class ExactBlockPreconditioner
    : public GeneralPurposeBlockPreconditioner<MATRIX>
  {
  public:
    ExactBlockPreconditioner() : GeneralPurposeBlockPreconditioner<MATRIX>() {}
 
    virtual ~ExactBlockPreconditioner() {}
 
    ExactBlockPreconditioner(const ExactBlockPreconditioner&) = delete;
 
    void operator=(const ExactBlockPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
    Preconditioner*& preconditioner_pt()
    {
      return this->Subsidiary_preconditioner_pt[0];
    }
  };
 
 
  // =================================================================
  // =================================================================
  template<typename MATRIX>
  class BlockAntiDiagonalPreconditioner
    : public BlockDiagonalPreconditioner<MATRIX>
  {
  protected:
    unsigned get_other_diag_ds(const unsigned& i, const unsigned& nblock) const
    {
      return nblock - i - 1;
    }
  };
 
 
  // =================================================================
  // =================================================================
  template<typename MATRIX>
  class DummyBlockPreconditioner
    : public GeneralPurposeBlockPreconditioner<MATRIX>
  {
  public:
    DummyBlockPreconditioner() : GeneralPurposeBlockPreconditioner<MATRIX>() {}
 
    ~DummyBlockPreconditioner() {}
 
    DummyBlockPreconditioner(const DummyBlockPreconditioner&) = delete;
 
    void operator=(const DummyBlockPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      z.build(r);
    }
 
    void setup()
    {
      // Set up the block look up schemes
      this->block_setup();
    }
  };
 
} // namespace oomph
#endif