biharmonic_preconditioner.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_BIHARMONIC_PRECONDITIONER_HEADER
#define OOMPH_BIHARMONIC_PRECONDITIONER_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/preconditioner.h"
#include "../generic/block_preconditioner.h"
#include "../generic/hijacked_elements.h"
#include "biharmonic_elements.h"
#include "../meshes/hermite_element_quad_mesh.template.h"
#include "../generic/SuperLU_preconditioner.h"
#include "../generic/general_purpose_preconditioners.h"
 
#ifdef OOMPH_HAS_HYPRE
#include "../generic/hypre_solver.h"
#endif
 
namespace oomph
{
#ifdef OOMPH_HAS_HYPRE
 
  //=============================================================================
  // defaults settings for the Hypre solver (AMG) when used as the approximate
  // linear solver for the Schur complement (non-compound) linear subsidiary
  // linear systems
  //=============================================================================
  namespace Biharmonic_schur_complement_Hypre_defaults
  {
    extern unsigned AMG_smoother;
 
    extern unsigned AMG_coarsening;
 
    extern unsigned N_cycle;
 
    extern double AMG_strength;
 
    extern double AMG_jacobi_damping;
 
    extern unsigned AMG_smoother_iterations;
 
    extern void set_defaults(HyprePreconditioner* hypre_prec_pt);
 
  } // namespace Biharmonic_schur_complement_Hypre_defaults
#endif
 
 
  //=============================================================================
  //=============================================================================
  class BiharmonicPreconditioner : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    BiharmonicPreconditioner()
    {
      // initialise both preconditioners to zero
      Sub_preconditioner_1_pt = 0;
      Sub_preconditioner_2_pt = 0;
      Hijacked_sub_block_preconditioner_pt = 0;
 
      // by default we use the inexact biharmonic preconditioner
      // and if possible use Hypre preconditioner
#ifdef OOMPH_HAS_HYPRE
      Preconditioner_type = 2;
#else
      Preconditioner_type = 1;
#endif
 
      // size mesh pt correctly
      this->set_nmesh(1);
      Bulk_element_mesh_pt = 0;
    }
 
    ~BiharmonicPreconditioner()
    {
      this->clean_up_memory();
    }
 
    // delete the subsidiary preconditioners and memory
    void clean_up_memory()
    {
      // delete the sub preconditioners
      delete Sub_preconditioner_1_pt;
      delete Sub_preconditioner_2_pt;
      delete Hijacked_sub_block_preconditioner_pt;
    }
 
    BiharmonicPreconditioner(const BiharmonicPreconditioner&) = delete;
 
    void operator=(const BiharmonicPreconditioner&) = delete;
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    unsigned& preconditioner_type()
    {
      return Preconditioner_type;
    }
 
    Mesh*& bulk_element_mesh_pt()
    {
      return Bulk_element_mesh_pt;
    }
 
  private:
    unsigned Preconditioner_type;
 
    Preconditioner* Sub_preconditioner_1_pt;
 
    Preconditioner* Sub_preconditioner_2_pt;
 
    Preconditioner* Hijacked_sub_block_preconditioner_pt;
 
    Mesh* Bulk_element_mesh_pt;
  };
 
 
  //=============================================================================
  //=============================================================================
  class ExactSubBiharmonicPreconditioner
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    ExactSubBiharmonicPreconditioner(BiharmonicPreconditioner* master_prec_pt,
                                     const bool& retain_all_blocks = false)
      : Retain_all_blocks(retain_all_blocks)
    {
      // Block mapping for ExactSubBiharmonicPreconditioner
      Vector<unsigned> block_lookup(3);
      block_lookup[0] = 0;
      block_lookup[1] = 1;
      block_lookup[2] = 2;
 
      // set as subsidiary block preconditioner
      this->turn_into_subsidiary_block_preconditioner(master_prec_pt,
                                                      block_lookup);
 
      // null the Sub preconditioner pt
      Sub_preconditioner_pt = 0;
    }
 
    ~ExactSubBiharmonicPreconditioner()
    {
      this->clean_up_memory();
    }
 
    void clean_up_memory()
    {
      delete Sub_preconditioner_pt;
      Sub_preconditioner_pt = 0;
    }
 
    ExactSubBiharmonicPreconditioner(const ExactSubBiharmonicPreconditioner&) =
      delete;
 
    void operator=(const ExactSubBiharmonicPreconditioner&) = delete;
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    //   private:
 
    // Pointer to the sub preconditioner
    Preconditioner* Sub_preconditioner_pt;
 
 
    bool Retain_all_blocks;
  };
 
  //=============================================================================
  //=============================================================================
  class InexactSubBiharmonicPreconditioner
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    InexactSubBiharmonicPreconditioner(BiharmonicPreconditioner* master_prec_pt,
                                       const bool use_amg)
    {
      // Block mapping for ExactSubBiharmonicPreconditioner
      Vector<unsigned> block_lookup(3);
      block_lookup[0] = 0;
      block_lookup[1] = 1;
      block_lookup[2] = 2;
 
      // set as subsidiary block preconditioner
      this->turn_into_subsidiary_block_preconditioner(master_prec_pt,
                                                      block_lookup);
 
      // zero all the preconditioner pt
      S_00_preconditioner_pt = 0;
      Lumped_J_11_preconditioner_pt = 0;
      Lumped_J_22_preconditioner_pt = 0;
 
      // store the preconditioner type
      Use_amg = use_amg;
    }
 
    ~InexactSubBiharmonicPreconditioner()
    {
      this->clean_up_memory();
    }
 
    void clean_up_memory()
    {
      // delete the sub preconditioner pt
      delete S_00_preconditioner_pt;
      S_00_preconditioner_pt = 0;
 
      // delete the lumped preconditioners
      delete Lumped_J_11_preconditioner_pt;
      Lumped_J_11_preconditioner_pt = 0;
      delete Lumped_J_22_preconditioner_pt;
      Lumped_J_22_preconditioner_pt = 0;
 
      // Number of block types
      unsigned n = Matrix_of_block_pointers.nrow();
 
      // delete the block matrices
      for (unsigned i = 0; i < n; i++)
      {
        for (unsigned j = 0; j < n; j++)
        {
          if (Matrix_of_block_pointers(i, j) != 0)
          {
            delete Matrix_of_block_pointers(i, j);
            Matrix_of_block_pointers(i, j) = 0;
          }
        }
      }
    }
 
    InexactSubBiharmonicPreconditioner(
      const InexactSubBiharmonicPreconditioner&) = delete;
 
    void operator=(const InexactSubBiharmonicPreconditioner&) = delete;
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    //    private:
 
    void compute_inexact_schur_complement();
 
    Preconditioner* S_00_preconditioner_pt;
 
    MatrixBasedLumpedPreconditioner<CRDoubleMatrix>*
      Lumped_J_11_preconditioner_pt;
 
    MatrixBasedLumpedPreconditioner<CRDoubleMatrix>*
      Lumped_J_22_preconditioner_pt;
 
    DenseMatrix<CRDoubleMatrix*> Matrix_of_block_pointers;
 
    CRDoubleMatrix* S_00_pt;
 
    unsigned Use_amg;
  };
} // namespace oomph
#endif