pseudo_elastic_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//====================================================================
#ifndef OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER
#define OOMPH_PSEUDO_ELASTIC_SUBSIDIARY_PRECONDITIONER
 
// includes
#include "../generic/problem.h"
#include "../generic/block_preconditioner.h"
#include "../generic/preconditioner.h"
#include "../generic/SuperLU_preconditioner.h"
#include "../generic/matrix_vector_product.h"
#include "../generic/general_purpose_preconditioners.h"
#include "../generic/general_purpose_block_preconditioners.h"
#ifdef OOMPH_HAS_HYPRE
#include "../generic/hypre_solver.h"
#endif
#ifdef OOMPH_HAS_TRILINOS
#include "../generic/trilinos_solver.h"
#endif
namespace oomph
{
  //=============================================================================
  //=============================================================================
  namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
  {
#ifdef OOMPH_HAS_HYPRE
 
    Preconditioner* get_elastic_preconditioner_hypre();
 
    Preconditioner* get_elastic_preconditioner();
 
#endif
 
#ifdef OOMPH_HAS_TRILINOS
 
    Preconditioner* get_elastic_preconditioner_trilinos_ml();
 
    Preconditioner* get_lagrange_multiplier_preconditioner();
#endif
  } // namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
 
 
 
 
  //=============================================================================
  //=============================================================================
  class PseudoElasticPreconditioner : public BlockPreconditioner<CRDoubleMatrix>
  {
    friend class PseudoElasticFSIPreconditioner;
 
  public:
    typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
 
    enum Elastic_preconditioner_type
    {
      Exact_block_preconditioner,
      Block_diagonal_preconditioner,
      Block_lower_triangular_preconditioner,
      Block_upper_triangular_preconditioner
    };
 
    PseudoElasticPreconditioner()
    {
      // null pointers
      Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0;
      Elastic_subsidiary_preconditioner_function_pt = 0;
      Elastic_preconditioner_pt = 0;
 
      // set defaults
      Use_inf_norm_of_s_scaling = true;
      E_preconditioner_type = Exact_block_preconditioner;
 
      // resize the Mesh_pt
      Lagrange_multiplier_mesh_pt = 0;
      Elastic_mesh_pt = 0;
    }
 
    virtual ~PseudoElasticPreconditioner()
    {
      this->clean_up_memory();
    }
 
    PseudoElasticPreconditioner(const PseudoElasticPreconditioner&) = 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 PseudoElasticPreconditioner&) = delete;*/
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      this->elastic_preconditioner_solve(r, z);
      this->lagrange_multiplier_preconditioner_solve(r, z);
    }
 
    void set_elastic_mesh(Mesh* mesh_pt)
    {
      Elastic_mesh_pt = mesh_pt;
    }
 
    void set_lagrange_multiplier_mesh(Mesh* mesh_pt)
    {
      Lagrange_multiplier_mesh_pt = mesh_pt;
    }
 
    void enable_inf_norm_of_s_scaling()
    {
      Use_inf_norm_of_s_scaling = true;
    }
 
    void disable_inf_norm_of_s_scaling()
    {
      Use_inf_norm_of_s_scaling = false;
    }
 
    void set_lagrange_multiplier_subsidiary_preconditioner(
      SubsidiaryPreconditionerFctPt prec_fn)
    {
      Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn;
    }
 
    void set_elastic_subsidiary_preconditioner(
      SubsidiaryPreconditionerFctPt prec_fn)
    {
      Elastic_subsidiary_preconditioner_function_pt = prec_fn;
    }
 
    Elastic_preconditioner_type& elastic_preconditioner_type()
    {
      return E_preconditioner_type;
    }
 
    void clean_up_memory();
 
  private:
    void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                                  DoubleVector& z);
 
    double Scaling;
 
    bool Use_inf_norm_of_s_scaling;
 
    Elastic_preconditioner_type E_preconditioner_type;
 
    unsigned Dim;
 
    Preconditioner* Elastic_preconditioner_pt;
 
    Vector<Preconditioner*> Lagrange_multiplier_preconditioner_pt;
 
    SubsidiaryPreconditionerFctPt
      Lagrange_multiplier_subsidiary_preconditioner_function_pt;
 
    SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt;
 
    Mesh* Elastic_mesh_pt;
 
    Mesh* Lagrange_multiplier_mesh_pt;
  };
 
 
 
  //=============================================================================
  //=============================================================================
  class PseudoElasticPreconditionerOld
    : public BlockPreconditioner<CRDoubleMatrix>
  {
    friend class PseudoElasticFSIPreconditioner;
 
  public:
    typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
 
    enum Elastic_preconditioner_type
    {
      Exact_block_preconditioner,
      Block_diagonal_preconditioner,
      Block_lower_triangular_preconditioner,
      Block_upper_triangular_preconditioner
    };
 
    PseudoElasticPreconditionerOld()
    {
      // null pointers
      Lagrange_multiplier_subsidiary_preconditioner_function_pt = 0;
      Elastic_subsidiary_preconditioner_function_pt = 0;
      Elastic_preconditioner_pt = 0;
 
      // set defaults
      Use_inf_norm_of_s_scaling = true;
      E_preconditioner_type = Exact_block_preconditioner;
 
      // resize the Mesh_pt
      Lagrange_multiplier_mesh_pt = 0;
      Elastic_mesh_pt = 0;
    }
 
    virtual ~PseudoElasticPreconditionerOld()
    {
      this->clean_up_memory();
    }
 
    PseudoElasticPreconditionerOld(const PseudoElasticPreconditionerOld&) =
      delete;
 
    /*void operator=(const PseudoElasticPreconditionerOld&) = delete;*/
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      this->elastic_preconditioner_solve(r, z);
      this->lagrange_multiplier_preconditioner_solve(r, z);
    }
 
    void set_elastic_mesh(Mesh* mesh_pt)
    {
      Elastic_mesh_pt = mesh_pt;
    }
 
    void set_lagrange_multiplier_mesh(Mesh* mesh_pt)
    {
      Lagrange_multiplier_mesh_pt = mesh_pt;
    }
 
    void enable_inf_norm_of_s_scaling()
    {
      Use_inf_norm_of_s_scaling = true;
    }
 
    void disable_inf_norm_of_s_scaling()
    {
      Use_inf_norm_of_s_scaling = false;
    }
 
    void set_lagrange_multiplier_subsidiary_preconditioner(
      SubsidiaryPreconditionerFctPt prec_fn)
    {
      Lagrange_multiplier_subsidiary_preconditioner_function_pt = prec_fn;
    }
 
    void set_elastic_subsidiary_preconditioner(
      SubsidiaryPreconditionerFctPt prec_fn)
    {
      Elastic_subsidiary_preconditioner_function_pt = prec_fn;
    }
 
    Elastic_preconditioner_type& elastic_preconditioner_type()
    {
      return E_preconditioner_type;
    }
 
    void clean_up_memory();
 
  private:
    void elastic_preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                                  DoubleVector& z);
 
    double Scaling;
 
    bool Use_inf_norm_of_s_scaling;
 
    Elastic_preconditioner_type E_preconditioner_type;
 
    unsigned Dim;
 
    Preconditioner* Elastic_preconditioner_pt;
 
    Vector<Preconditioner*> Lagrange_multiplier_preconditioner_pt;
 
    SubsidiaryPreconditionerFctPt
      Lagrange_multiplier_subsidiary_preconditioner_function_pt;
 
    SubsidiaryPreconditionerFctPt Elastic_subsidiary_preconditioner_function_pt;
 
    Mesh* Elastic_mesh_pt;
 
    Mesh* Lagrange_multiplier_mesh_pt;
  };
 
 
 
 
  //=============================================================================
  //=============================================================================
  class PseudoElasticPreconditionerSubsidiaryPreconditionerOld
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
 
    PseudoElasticPreconditionerSubsidiaryPreconditionerOld()
    {
      Scaling = 1.0;
      Preconditioner_pt = 0;
      Subsidiary_preconditioner_function_pt = 0;
    }
 
    ~PseudoElasticPreconditionerSubsidiaryPreconditionerOld()
    {
      this->clean_up_memory();
    }
 
    PseudoElasticPreconditionerSubsidiaryPreconditionerOld(
      const PseudoElasticPreconditionerSubsidiaryPreconditionerOld&) = delete;
 
    /*void operator=(const
     PseudoElasticPreconditionerSubsidiaryPreconditionerOld&) = delete;*/
 
    // Setup the preconditioner
    void setup();
 
    // Apply the preconditioner
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    double& scaling()
    {
      return Scaling;
    }
 
    void set_subsidiary_preconditioner_function(
      SubsidiaryPreconditionerFctPt sub_prec_fn)
    {
      Subsidiary_preconditioner_function_pt = sub_prec_fn;
    }
 
  private:
    void clean_up_memory()
    {
      delete Preconditioner_pt;
      Preconditioner_pt = 0;
    }
 
    // the augmentation scaling
    double Scaling;
 
    Preconditioner* Preconditioner_pt;
 
    SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt;
  }; // end of PseudoElasticPreconditionerSubsidiaryPreconditionerOld
 
 
 
 
  //=============================================================================
  //=============================================================================
  class PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();
 
    PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld()
      : BlockPreconditioner<CRDoubleMatrix>()
    {
      // null the subsidiary preconditioner function pointer
      Subsidiary_preconditioner_function_pt = 0;
 
      // default to block diagonal
      Method = 0;
 
      // default scaling = 1.0
      Scaling = 1.0;
    };
 
    ~PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld()
    {
      this->clean_up_memory();
    }
 
    PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld(
      const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld&) =
      delete;
 
    /*void operator=
     (const PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld&) =
     delete;*/
 
    void clean_up_memory();
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& res, DoubleVector& z);
 
    void set_subsidiary_preconditioner_function(
      SubsidiaryPreconditionerFctPt sub_prec_fn)
    {
      Subsidiary_preconditioner_function_pt = sub_prec_fn;
    };
 
    void use_block_diagonal_approximation()
    {
      Method = 0;
    }
 
    void use_upper_triangular_approximation()
    {
      Method = 1;
    }
 
    void use_lower_triangular_approximation()
    {
      Method = 2;
    }
 
    double& scaling()
    {
      return Scaling;
    }
 
  private:
    Vector<PseudoElasticPreconditionerSubsidiaryPreconditionerOld*>
      Diagonal_block_preconditioner_pt;
 
    DenseMatrix<MatrixVectorProduct*> Off_diagonal_matrix_vector_products;
 
    unsigned Method;
 
    SubsidiaryPreconditionerFctPt Subsidiary_preconditioner_function_pt;
 
    double Scaling;
  };
 
 
 
 
  // /*
  //=============================================================================
  //=============================================================================
  class PseudoElasticPreconditionerScalingHelperOld
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    PseudoElasticPreconditionerScalingHelperOld(
      BlockPreconditioner<CRDoubleMatrix>* master_prec_pt,
      CRDoubleMatrix* matrix_pt,
      Vector<unsigned>& dof_list,
      const Mesh* const solid_mesh_pt,
      const OomphCommunicator* comm_pt)
    {
      // turn into a subisiary preconditioner
      this->turn_into_subsidiary_block_preconditioner(master_prec_pt, dof_list);
 
      // all dofs are of the same block type
      Vector<unsigned> dof_to_block_map(dof_list.size(), 0);
 
      // store the matrix_pt
      set_matrix_pt(matrix_pt);
 
      // set the mesh
      this->set_nmesh(1);
      this->set_mesh(0, solid_mesh_pt);
 
      // set the communicator pointer
      this->set_comm_pt(comm_pt);
 
      // call block_setup(...)
      this->block_setup(dof_to_block_map);
    }
 
    ~PseudoElasticPreconditionerScalingHelperOld()
    {
      this->clear_block_preconditioner_base();
    }
 
    PseudoElasticPreconditionerScalingHelperOld(
      const PseudoElasticPreconditionerScalingHelperOld&) = delete;
 
    /*void operator=(const PseudoElasticPreconditionerScalingHelperOld&) =
     * delete;*/
 
    double s_inf_norm()
    {
      CRDoubleMatrix* m_pt = new CRDoubleMatrix;
      this->get_block(0, 0, *m_pt);
      double s_inf_norm = m_pt->inf_norm();
      delete m_pt;
      return s_inf_norm;
    }
 
    // broken preconditioner setup
    void setup()
    {
      std::ostringstream error_message;
      error_message << "This method is intentionally broken. This class is not "
                       "a functioning "
                    << "preconditioner.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // broken preconditioner solve
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      std::ostringstream error_message;
      error_message << "This method is intentionally broken. This class is not "
                       "a functioning "
                    << "preconditioner.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
  }; // end of PseudoElasticPreconditionerScalingHelperOld
  // */
 
} // namespace oomph
#endif