solid_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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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_SOLID_PRECONDITIONERS_HEADER
#define OOMPH_SOLID_PRECONDITIONERS_HEADER


// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif


// oomphlib headers
#include "../generic/matrices.h"
#include "../generic/assembly_handler.h"
#include "../generic/problem.h"
#include "../generic/block_preconditioner.h"
#include "../generic/preconditioner.h"
#include "../generic/SuperLU_preconditioner.h"
#include "../generic/matrix_vector_product.h"


namespace oomph
{



//===========================================================================
/// \short The least-squares commutator (LSC; formerly BFBT) 
/// preconditioner. It uses blocks corresponding to the displacement/position 
/// and pressure unknowns, i.e. there are a total of 2x2 blocks, 
/// and all displacement/position components are treated as a 
/// single block of unknowns.
/// 
/// Here are the details: An "ideal" preconditioner
/// would solve the saddle point system
/// \f[
/// \left( 
/// \begin{array}{cc}
/// {\bf F} & {\bf G} \\ {\bf D} & {\bf 0} 
/// \end{array} 
/// \right)
/// \left( 
/// \begin{array}{c}
/// {\bf z}_u \\ {\bf z}_p
/// \end{array} 
/// \right) =
/// \left( 
/// \begin{array}{c}
/// {\bf r}_u \\ {\bf r}_p
/// \end{array} 
/// \right)
/// \f]
/// where \f$ {\bf F}\f$,  \f$ {\bf G} \f$, and \f$ {\bf D}\f$ are 
/// the blocks that arise in the Jacobian of the pressure-based
/// equations of linear and nonlinear elasticity (with dofs in order
/// of displacement/position and pressure).
/// The use of this preconditioner would ensure the convergence
/// of any iterative linear solver in a single iteration but its
/// application is, of course, exactly as expensive as a direct solve.
/// The LSC/BFBT preconditioner replaces the exact Jacobian by 
/// a block-triangular approximation
/// \f[
/// \left( 
/// \begin{array}{cc}
/// {\bf F} & {\bf G} \\ {\bf 0} & -{\bf M}_s 
/// \end{array} 
/// \right) 
/// \left( 
/// \begin{array}{c}
/// {\bf z}_u \\ {\bf z}_p
/// \end{array} 
/// \right) =
/// \left( 
/// \begin{array}{c}
/// {\bf r}_u \\ {\bf r}_p
/// \end{array} 
/// \right),
/// \f]
/// where \f${\bf M}_s\f$ is an approximation to the pressure 
/// Schur-complement \f$ {\bf S} = {\bf D} {\bf F}^{-1}{\bf G}. \f$
/// This system can be solved in two steps:
/// -# Solve the second row for \f$ {\bf z}_p\f$ via
///    \f[ 
///    {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p
///    \f]
/// -# Given \f$ {\bf z}_p \f$ , solve the first row for \f$ {\bf z}_u\f$ via
///    \f[ 
///    {\bf z}_u = {\bf F}^{-1} \big( {\bf r}_u - {\bf G} {\bf z}_p \big)
///    \f]
/// .
/// In the LSC/BFBT preconditioner, the action of the inverse pressure
/// Schur complement 
/// \f[
/// {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p
/// \f]
/// is approximated by
/// \f[
/// {\bf z}_p = - 
/// \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1}
/// \big({\bf D} \widehat{\bf Q}^{-1}{\bf F} \widehat{\bf Q}^{-1}{\bf G}\big) 
/// \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1}
/// {\bf r}_p,
/// \f]
/// where  \f$ \widehat{\bf Q} \f$ is the diagonal of the displacement/position
/// mass matrix. The evaluation of this expression involves
/// two linear solves involving the matrix
/// \f[
/// {\bf P} = \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)
/// \f]
/// which has the character of a matrix arising from the discretisation 
/// of a Poisson problem on the pressure space. We also have
/// to evaluate matrix-vector products with the matrix 
/// \f[
/// {\bf E}={\bf D}\widehat{\bf Q}^{-1}{\bf F}\widehat{\bf Q}^{-1}{\bf G}
/// \f]
/// Details of the theory can be found in "Finite Elements and 
/// Fast Iterative Solvers with Applications in Incompressible Fluid 
/// Dynamics" by Howard C. Elman, David J. Silvester, and Andrew J. Wathen,
/// published by Oxford University Press, 2006.
/// 
/// In our implementation of the preconditioner, the linear systems
/// can either be solved "exactly", using SuperLU (in its incarnation
/// as an exact preconditioner; this is the default) or by any 
/// other Preconditioner (inexact solver) specified via the access functions
/// \code
/// PressureBasedSolidLSCPreconditioner::set_f_preconditioner(...)
/// \endcode
/// or 
/// \code
/// PressureBasedSolidLSCPreconditioner::set_p_preconditioner(...)
/// \endcode
//===========================================================================
  class PressureBasedSolidLSCPreconditioner :
  public BlockPreconditioner<CRDoubleMatrix>
  {
   
   public :
    
    /// Constructor - sets defaults for control flags
    PressureBasedSolidLSCPreconditioner() : BlockPreconditioner<CRDoubleMatrix>()
    {
     // Flag to indicate that the preconditioner has been setup
     // previously -- if setup() is called again, data can
     // be wiped.
     Preconditioner_has_been_setup = false;

     // By default we use SuperLU for both p and f blocks
     Using_default_p_preconditioner=true;
     Using_default_f_preconditioner=true;

     // resize the mesh pt
     // note: meaningless if subsidiary preconditioner
     this->set_nmesh(1);
     Solid_mesh_pt = 0;

     // Set default preconditioners (inexact solvers) -- they are 
     // members of this class!
     P_preconditioner_pt = 0;
     F_preconditioner_pt = 0;

     // Flag to determine if mass matrix diagonal Q^{-1}
     // is used for scaling.
     P_matrix_using_scaling = true;

     // set Doc_time to false
     Doc_time = false;

     // null the off diagonal Block matrix pt
     Bt_mat_vec_pt = 0;

     // null the F matrix vector product helper
     F_mat_vec_pt = 0;

     // null the QBt matrix vector product pt
     QBt_mat_vec_pt = 0;

     // null the E matrix vector product helper
     E_mat_vec_pt = 0;

     // by default we do not form the E matrix (BQFQBt)
     Form_BFBt_product = false;
    }

   /// Destructor
   ~PressureBasedSolidLSCPreconditioner()
    {
     clean_up_memory();
    }

   /// Broken copy constructor
   PressureBasedSolidLSCPreconditioner(const PressureBasedSolidLSCPreconditioner&)
    {
     BrokenCopy::broken_copy("PressureBasedSolidLSCPreconditioner");
    }

   /// Broken assignment operator
//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 PressureBasedSolidLSCPreconditioner&)
    {
     BrokenCopy::broken_assign("PressureBasedSolidLSCPreconditioner");
     }*/

   /// Setup the preconditioner
   void setup();

   /// Apply preconditioner to Vector r
   void preconditioner_solve(const DoubleVector&r, DoubleVector &z);
   
   /// specify the mesh containing the mesh containing the 
   /// block-preconditionable solid elements. The dimension of the
   /// problem must also be specified.
   void set_solid_mesh(Mesh* mesh_pt)
    { 
     Solid_mesh_pt = mesh_pt;
    }

   /// \short Enable mass matrix diagonal scaling in the 
   /// Schur complement approximation
   void enable_p_matrix_scaling() {P_matrix_using_scaling=true;}

    /// \short Enable mass matrix diagonal scaling in the 
   /// Schur complement approximation
   void disable_p_matrix_scaling() {P_matrix_using_scaling=false;}

   /// \short Return whether the mass matrix is using diagonal
   /// scaling or not
   bool is_p_matrix_using_scaling() const {return P_matrix_using_scaling;}

   /// Function to set a new pressure matrix preconditioner (inexact solver)
   void set_p_preconditioner(Preconditioner* new_p_preconditioner_pt)
   {
    // If the default preconditioner has been used
    // clean it up now...
    if (Using_default_p_preconditioner)
     {
      delete P_preconditioner_pt;
     }
    P_preconditioner_pt = new_p_preconditioner_pt;
    Using_default_p_preconditioner = false;
   }

   /// \short Function to (re-)set pressure matrix preconditioner  (inexact 
   /// solver) to SuperLU
   void set_p_superlu_preconditioner()
   {
    if (!Using_default_p_preconditioner)
     {
      P_preconditioner_pt = new SuperLUPreconditioner;
      Using_default_p_preconditioner = true;
     }
   }

   /// Function to set a new momentum matrix preconditioner (inexact solver)
   void set_f_preconditioner(Preconditioner* new_f_preconditioner_pt)
   {
    // If the default preconditioner has been used
    // clean it up now...
    if (Using_default_f_preconditioner)
     {
      delete F_preconditioner_pt;
     }
    F_preconditioner_pt = new_f_preconditioner_pt;
    Using_default_f_preconditioner = false;
   }

   ///\short Function to (re-)set momentum matrix preconditioner (inexact 
   /// solver) to SuperLU
   void set_f_superlu_preconditioner()
   {
    if (!Using_default_f_preconditioner)
     {
      F_preconditioner_pt = new SuperLUPreconditioner;
      Using_default_f_preconditioner = true;
     }
   }

 
   ///Enable documentation of time
   void enable_doc_time() {Doc_time = true;}

   ///Disable documentation of time
   void disable_doc_time() {Doc_time = false;}

   /// \short If this function is called then:
   /// in setup(...) : BFBt is computed.
   /// in preconditioner_solve(...) : a single matrix vector product with 
   /// BFBt is performed.
   void enable_form_BFBt_product() {Form_BFBt_product = true; }

   /// \short if this function is called  then:
   /// in setup(...) : the matrices B, F are assembled and stored 
   /// (the default behaviour) .
   /// in preconditioner_solve(...) : a sequence of matrix vector products
   /// with B, F, and Bt is performed.
   /// (Note: in this discussion no scaling was considered but B and Bt 
   ///  are replaced with BQ and QBt with scaling)
   void disable_form_BFBt_product() {Form_BFBt_product = false;}

   /// \short Helper function to delete preconditioner data.
   void clean_up_memory();

    private:

   // oomph-lib objects
   // -----------------

   // Pointers to preconditioner (=inexact solver) objects
   // -----------------------------------------------------
   /// Pointer to the 'preconditioner' for the pressure matrix
   Preconditioner* P_preconditioner_pt;

   /// Pointer to the 'preconditioner' for the F matrix
   Preconditioner* F_preconditioner_pt;

   /// flag indicating whether the default F preconditioner is used
   bool Using_default_f_preconditioner;

   /// flag indicating whether the default P preconditioner is used
   bool Using_default_p_preconditioner;

   /// \short Control flag is true if the preconditioner has been setup
   /// (used so we can wipe the data when the preconditioner is
   /// called again)
   bool Preconditioner_has_been_setup;

   /// \short Control flag is true if mass matrix diagonal scaling
   /// is used in the Schur complement approximation
   bool P_matrix_using_scaling;

   /// Helper function to assemble the diagonal of the 
   /// mass matrix from the elemental contributions defined in
   /// PressureBasedSolidEquations<DIM>::get_mass_matrix_diagonal(...).
   CRDoubleMatrix* assemble_mass_matrix_diagonal();

   /// \short Boolean indicating whether the momentum system preconditioner 
   /// is a block preconditioner
   bool F_preconditioner_is_block_preconditioner;

   /// Set Doc_time to true for outputting results of timings
   bool Doc_time;

   /// MatrixVectorProduct operator for F if BFBt is not to be formed.
   MatrixVectorProduct* F_mat_vec_pt;

   /// MatrixVectorProduct operator for QBt if BFBt is not to be formed.
   MatrixVectorProduct* QBt_mat_vec_pt;

   /// MatrixVectorProduct operator for Bt;
   MatrixVectorProduct* Bt_mat_vec_pt;

   /// MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.
   MatrixVectorProduct* E_mat_vec_pt;

   /// \short indicates whether BFBt should be formed or the component matrices
   /// should be retained.
   /// If true then:
   /// in setup(...) : BFBt is computed.
   /// in preconditioner_solve(...) : a single matrix vector product with 
   /// BFBt is performed.
   /// if false then:
   /// in setup(...) : the matrices B, F are assembled and stored.
   /// in preconditioner_solve(...) : a sequence of matrix vector products
   /// with B, F, and Bt is performed.
   /// (Note: in this discussion no scaling was considered but B and Bt 
   ///  are replaced with BQ and QBt with scaling)
   bool Form_BFBt_product;

   /// \short the pointer to the mesh of block preconditionable solid
   /// elements.
   Mesh* Solid_mesh_pt;
  };



///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////



//============================================================================
/// \short The exact solid preconditioner. This extracts 2x2 blocks
/// (corresponding to the displacement/position and pressure unknowns) 
/// and uses these to
/// build a single preconditioner matrix for testing purposes.
/// Iterative solvers should converge in a single step if this is used.
/// If it doesn't something is wrong in the setup of the block matrices.
//=============================================================================
 template<typename MATRIX>
  class PressureBasedSolidExactPreconditioner : public BlockPreconditioner<MATRIX>
  {

   public :
     
    /// Constructor - do nothing
    PressureBasedSolidExactPreconditioner() : BlockPreconditioner<MATRIX>(){}
   
   
   /// Destructor - do nothing
   ~PressureBasedSolidExactPreconditioner(){}

   
   /// Broken copy constructor
   PressureBasedSolidExactPreconditioner(const PressureBasedSolidExactPreconditioner&)
    {
     BrokenCopy::broken_copy("PressureBasedSolidExactPreconditioner");
    }
   

   /// Broken assignment operator   
   /*void operator=(const PressureBasedSolidExactPreconditioner&)
    {
     BrokenCopy::broken_assign("PressureBasedSolidExactPreconditioner");
     }*/
   
   
   /// Setup the preconditioner
   void setup();
   
   /// Apply preconditioner to r
   void preconditioner_solve(const Vector<double>&r,
                             Vector<double> &z);

   protected :
    
    /// Preconditioner matrix
    MATRIX P_matrix;

  };
 
}
#endif