navier_stokes_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_NAVIER_STOKES_PRECONDITIONERS_HEADER
#define OOMPH_NAVIER_STOKES_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"
#include "navier_stokes_elements.h"
#include "refineable_navier_stokes_elements.h"


namespace oomph
{




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


//======start_of_namespace============================================
/// Namespace for exact solution for pressure advection diffusion 
/// problem
//====================================================================
 namespace PressureAdvectionDiffusionValidation
 {
  
  /// Flag for solution
  extern unsigned Flag;
    
  /// Peclet number -- overwrite with actual Reynolds number
  extern double Peclet;
    
  /// Wind
  extern void wind_function(const Vector<double>& x, Vector<double>& wind);
  
  /// Exact solution as a Vector
  extern void get_exact_u(const Vector<double>& x, Vector<double>& u);
  
  /// Exact solution as a scalar
  extern void get_exact_u(const Vector<double>& x, double& u);
  
  /// Source function required to make the solution above an exact solution 
  extern double source_function(const Vector<double>& x_vect);
  
 } // end of namespace
 
 
 
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
 
 
 
//=============================================================
/// A class that is used to define the functions used to
/// assemble the elemental contributions to the pressure
/// advection diffusion problem used by the Fp preconditioner.
//===============================================================
 class FpPreconditionerAssemblyHandler : public AssemblyHandler
 {
  
   public:
  
  /// Constructor. Pass spatial dimension
  FpPreconditionerAssemblyHandler(const unsigned& ndim) // : Ndim(ndim)
   {}
   
   /// \short Empty virtual destructor
   virtual ~FpPreconditionerAssemblyHandler() {}
   
   /// \short Return the contribution to the residuals of the element elem_pt
   void get_residuals(GeneralisedElement* const &elem_pt,
                      Vector<double> &residuals)
   {
    unsigned n_dof=elem_pt->ndof();
    for (unsigned i=0;i<n_dof;i++)
     {
      residuals[i]=0.0;
     }

    dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
     elem_pt)->fill_in_pressure_advection_diffusion_residuals(residuals);    
   }
   
   /// \short Calculate the elemental Jacobian matrix "d equation 
   /// / d variable" for elem_pt. 
   void get_jacobian(GeneralisedElement* const &elem_pt,
                     Vector<double> &residuals, 
                     DenseMatrix<double> &jacobian)
   {
    // Initialise
    unsigned n_dof=elem_pt->ndof();
    for (unsigned i=0;i<n_dof;i++)
     {
      residuals[i]=0.0;
      for (unsigned j=0;j<n_dof;j++)
       {
        jacobian(i,j)=0.0;
       }
     }
    
    dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
     elem_pt)->fill_in_pressure_advection_diffusion_jacobian(residuals,
                                                             jacobian); 
   }
   
   private:
   
   /// Spatial dimension of problem
   //unsigned Ndim; // cgj: Ndim is never used.
   
 };
 
 
 
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
 
 
//===========================================================================
/// \short Auxiliary Problem that can be used to assemble the
/// pressure advection diffusion matrix needed by the FpPreconditoner
//===========================================================================
 template<class ELEMENT>
  class FpPressureAdvectionDiffusionProblem : public virtual Problem
 {
  
   public:
  
  /// Constructor: Pass Navier-Stokes mesh and pointer to orig problem
  FpPressureAdvectionDiffusionProblem(Mesh* navier_stokes_mesh_pt,
   Problem* orig_problem_pt)
   {
    // Pass across mesh -- boundary conditions have already
    // been applied and the equations have been enumerated.
    mesh_pt()=navier_stokes_mesh_pt;
    
    // Store pointer to orig problem so we can re-assign the
    // orig eqn numbers when we're done
    Orig_problem_pt=orig_problem_pt;

    // Get the spatial dimension
    Ndim=mesh_pt()->finite_element_pt(0)->dim();
    
    // Set new assembly handler
    this->assembly_handler_pt()=new FpPreconditionerAssemblyHandler(Ndim);
   }
  
  
  /// Get the pressure advection diffusion Jacobian
  void get_pressure_advection_diffusion_jacobian(CRDoubleMatrix& jacobian)
  {
   // Pin all non-pressure dofs 
   pin_all_non_pressure_dofs(); 
   
   // Re-assign the equation numbers for all elements in Navier-Stokes mesh
   unsigned n_dof=assign_eqn_numbers();
   oomph_info << "Eqn numbers after pinning veloc dofs: " 
              << n_dof << std::endl;
   
   // Get "Jacobian" of the modified system
   DoubleVector dummy_residuals;
   this->get_jacobian(dummy_residuals,jacobian);
   
   // Reset pin status
   reset_pin_status();
   
   // Reassign equation numbering of original problem
   oomph_info << "Eqn numbers in orig problem after re-assignment: " <<
    Orig_problem_pt->assign_eqn_numbers() << std::endl;
  }
  
  
  
  /// \short Reset pin status of all values
  void reset_pin_status()
  {
   // Reset pin status for nodes
   unsigned nnod=mesh_pt()->nnode();
   for (unsigned j=0;j<nnod;j++)
    {
     Node* nod_pt=mesh_pt()->node_pt(j);
     unsigned nval=nod_pt->nvalue();
     for (unsigned i=0;i<nval;i++)
      {
       nod_pt->eqn_number(i)=Eqn_number_backup[nod_pt][i];
      }
    }
   
   // ... and internal data
   unsigned nelem=mesh_pt()->nelement();
   for (unsigned e=0;e<nelem;e++)
    {
     // Get actual element
     ELEMENT* el_pt=dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(e));
     

#ifdef PARANOID
     if (el_pt == 0)
      {
       std::ostringstream error_message;
       error_message << "Navier Stokes mesh must contain only Navier Stokes"
                     << "bulk elements\n";
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
#endif
     
     unsigned nint=el_pt->ninternal_data();
     for (unsigned j=0;j<nint;j++)
      {
       Data* data_pt=el_pt->internal_data_pt(j);
       unsigned nvalue=data_pt->nvalue();
       for (unsigned i=0;i<nvalue;i++)
        {
         data_pt->eqn_number(i)=Eqn_number_backup[data_pt][i];
        }
      }
    }

   // Free up storage
   Eqn_number_backup.clear();
  }
  
  

  /// \short Pin all non-pressure dofs and (if boolean flag is set to true
  /// also all pressure dofs along domain boundaries -- only used 
  /// for validation)
  void pin_all_non_pressure_dofs(const bool& 
                                 set_pressure_bc_for_validation=false)
  {

   // Backup pin status and then pin all non-pressure degrees of freedom
   unsigned nelem=mesh_pt()->nelement();
   for (unsigned e=0;e<nelem;e++)
    {
     // Get actual element (needed because different elements
     // store nodal pressure in different places)
     ELEMENT* el_pt=dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(e));
     

#ifdef PARANOID
     if (el_pt == 0)
      {
       std::ostringstream error_message;
       error_message << "Navier Stokes mesh must contain only Navier Stokes"
                     << "bulk elements\n";
       throw OomphLibError(
        error_message.str(),
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      }
#endif

     // Check if element has internal pressure representation
     // usually discontinuous -- preconditioner doesn't work for that case!
     if (el_pt->p_nodal_index_nst()<0)
      {
       std::ostringstream error_message;
       error_message << "Cannot use Fp preconditioner with discontinuous\n"
                     << "pressures.\n";
       throw OomphLibError(error_message.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }

     // Loop over internal data and pin the values (having established that
     // pressure dofs aren't amongst those)
     unsigned nint=el_pt->ninternal_data();
     for (unsigned j=0;j<nint;j++)
      {
       Data* data_pt=el_pt->internal_data_pt(j);
       if (Eqn_number_backup[data_pt].size()==0)
        {
         unsigned nvalue=data_pt->nvalue();
         Eqn_number_backup[data_pt].resize(nvalue);
         for (unsigned i=0;i<nvalue;i++)
          {
           // Backup
           Eqn_number_backup[data_pt][i]=data_pt->eqn_number(i);
           
           // Pin everything
           data_pt->pin(i);
          }
        }
      }
     
     // Now deal with nodal values
     unsigned nnod=el_pt->nnode();
     for (unsigned j=0;j<nnod;j++)
      {
       Node* nod_pt=el_pt->node_pt(j);
       if (Eqn_number_backup[nod_pt].size()==0)
        {
         unsigned nvalue=nod_pt->nvalue();
         Eqn_number_backup[nod_pt].resize(nvalue);
         for (unsigned i=0;i<nvalue;i++)
          {
           // Pin everything apart from the nodal pressure 
           // value
           if (int(i)!=el_pt->p_nodal_index_nst())
            {
             // Backup and pin
             Eqn_number_backup[nod_pt][i]=nod_pt->eqn_number(i);
             nod_pt->pin(i);
            }
           // Else it's a pressure value
           else 
            {
             // Exclude non-nodal pressure based elements
             if (el_pt->p_nodal_index_nst()>=0)
              {
               // Backup
               Eqn_number_backup[nod_pt][i]=nod_pt->eqn_number(i);
              }
            }
          }
        }
   
       // Set wind
       if (set_pressure_bc_for_validation)
        {
         Vector<double> x(2);
         x[0]=nod_pt->x(0);
         x[1]=nod_pt->x(1);
         Vector<double> u(2);
         PressureAdvectionDiffusionValidation::wind_function(x,u);
         nod_pt->set_value(el_pt->u_index_nst(0),u[0]);
         nod_pt->set_value(el_pt->u_index_nst(1),u[1]);
        }
      }
    }
  }


  /// \short Validate pressure advection diffusion problem and doc exact
  /// and computed solution in directory specified in DocInfo object
  void validate(DocInfo& doc_info)
  {
   
   oomph_info << "\n\n==============================================\n\n";
   oomph_info << "Doing validation for pressure adv diff problem\n";
   
   // Choose exact solution
   PressureAdvectionDiffusionValidation::Flag=1;

   // Pin all non-pressure dofs and pressure dofs along boundaries for
   // validation
   bool set_pressure_bc_for_validation=true; 
   pin_all_non_pressure_dofs(set_pressure_bc_for_validation); 
      

   // Set Peclet number
   PressureAdvectionDiffusionValidation::Peclet=
    dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(0))
    ->re();
   
   // Loop over all elements and set source function pointer
   unsigned nel=mesh_pt()->nelement();
   for (unsigned e=0;e<nel;e++)
    {
     dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(e))
      ->source_fct_for_pressure_adv_diff()=
      &PressureAdvectionDiffusionValidation::source_function;
    }
   
   // Re-assign the equation numbers for all elements in Navier-Stokes mesh
   oomph_info << "Eqn numbers after pinning veloc dofs: " 
              << assign_eqn_numbers() << std::endl;
   
   // Attach Robin BC elements
   unsigned nbound=mesh_pt()->nboundary(); 
   if (PressureAdvectionDiffusionValidation::Flag==1)
    {
     // Loop over all boundaries of Navier Stokes mesh
     for (unsigned b=0;b<nbound;b++)
      {
       // How many bulk elements are adjacent to boundary b?
       unsigned n_element = mesh_pt()->nboundary_element(b);
       
       // Loop over the bulk elements adjacent to boundary b?
       for(unsigned e=0;e<n_element;e++)
        {
         TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
          dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
           mesh_pt()->boundary_element_pt(b,e));
         
         //What is the index of the face of the bulk element e on bondary b
         int face_index=mesh_pt()->face_index_at_boundary(b,e);
         
         // Build face element
         bulk_elem_pt->build_fp_press_adv_diff_robin_bc_element(face_index);
         
        } //end of loop over bulk elements adjacent to boundary b
      }
    }

   
   // Loop over all elements and set source function pointer
      std::ofstream outfile;
   outfile.open("robin_elements.dat");
   for (unsigned e=0;e<nel;e++)
    {
     dynamic_cast<NavierStokesEquations<2>*>(mesh_pt()->element_pt(e))->
      output_pressure_advection_diffusion_robin_elements(outfile);   
    }
   outfile.close();

   // Solve it
   newton_solve();
   
   // And output the solution...
   doc_solution(doc_info);

   // Kill Robin BC elements
   if (PressureAdvectionDiffusionValidation::Flag==1)
    {
     // Loop over all boundaries of Navier Stokes mesh
     for (unsigned b=0;b<nbound;b++)
      {
       // How many bulk elements are adjacent to boundary b?
       unsigned n_element = mesh_pt()->nboundary_element(b);
       
       // Loop over the bulk elements adjacent to boundary b?
       for(unsigned e=0;e<n_element;e++)
        {
         
         TemplateFreeNavierStokesEquationsBase* bulk_elem_pt = 
          dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
           mesh_pt()->boundary_element_pt(b,e));
         
         // Kill associated Robin elements
         bulk_elem_pt->delete_pressure_advection_diffusion_robin_elements();
         
        } //end of loop over bulk elements adjacent to boundary b
      }
    }
   
   // Reset pin status
   reset_pin_status();

   // Reassign equation numbering of original problem
   oomph_info << "Eqn numbers in orig problem after re-assignment: " <<
    Orig_problem_pt->assign_eqn_numbers() << std::endl;


   oomph_info << "Done validation for pressure adv diff problem\n";
   oomph_info << "\n\n==============================================\n\n";
  }
 


  /// \short Doc solution (only needed during development -- kept alive
  /// in case validation is required during code maintenance)
  void doc_solution(DocInfo& doc_info)
  { 
   std::ofstream some_file;
   std::ostringstream filename;
   
   // Number of plot points
   unsigned npts;
   npts=5; 
   
   // Check value of FE solution in first element
   Vector<double> s(Ndim,0.0);
   Vector<double> x(Ndim);   
   double p=dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(0))
    ->interpolated_p_nst(s);   
   dynamic_cast<ELEMENT*>(mesh_pt()->finite_element_pt(0))
    ->interpolated_x(s,x);

   // Get offset-free exact solution
   double u_exact=0;
   PressureAdvectionDiffusionValidation::get_exact_u(x,u_exact);

   // Adjust offset
   unsigned nnode=mesh_pt()->nnode();
   for (unsigned j=0;j<nnode;j++)
    {
     if (mesh_pt()->node_pt(j)->nvalue()==3)
      {
       *(mesh_pt()->node_pt(j)->value_pt(2))-=p-u_exact;
      }
    }

   // Output solution 
   filename << doc_info.directory() << "/fp_soln" << doc_info.number() 
            << ".dat";
   some_file.open(filename.str().c_str());
   some_file.precision(20);
   mesh_pt()->output(some_file,npts);
   some_file.close();
   
   filename.str("");
   filename << doc_info.directory() << "/fp_exact_soln" 
            << doc_info.number() << ".dat";
   some_file.open(filename.str().c_str());
   some_file.precision(20);
   mesh_pt()->output_fct(some_file,npts,
                         PressureAdvectionDiffusionValidation::get_exact_u); 
   some_file.close();
   
  }
  
  private:

  /// The spatial dimension
  unsigned Ndim;
  
  /// \short Pointer to orig problem (required so we can re-assign 
  /// eqn numbers)
  Problem* Orig_problem_pt;

  /// Map to store original equation numbers
  std::map<Data*,std::vector<int> > Eqn_number_backup;

};



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




//===========================================================================
/// \short The least-squares commutator (LSC; formerly BFBT) Navier Stokes 
/// preconditioner. It uses blocks corresponding to the velocity 
/// and pressure unknowns, i.e. there are a total of 2x2 blocks, 
/// and all velocity components are treated as a single block of unknowns.
/// 
/// Here are the details: An "ideal" Navier-Stokes preconditioner
/// would solve the 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$ is the momentum block,  \f$ {\bf G} \f$ the
/// discrete gradient operator, and \f$ {\bf D}\f$ the discrete
/// divergence operator. (For unstabilised elements, we have 
/// \f$ {\bf D} = {\bf G}^T \f$ and in much of the literature
/// the divergence matrix is denoted by \f$ {\bf B} \f$ .)
/// 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 velocity
/// 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
/// NavierStokesSchurComplementPreconditioner::set_f_preconditioner(...)
/// \endcode
/// or 
/// \code
/// NavierStokesSchurComplementPreconditioner::set_p_preconditioner(...)
/// \endcode
//===========================================================================
 class NavierStokesSchurComplementPreconditioner :
 public BlockPreconditioner<CRDoubleMatrix>
 {
  
   public :
  
    /// Constructor - sets defaults for control flags
    NavierStokesSchurComplementPreconditioner(Problem* problem_pt) :
     BlockPreconditioner<CRDoubleMatrix>(), Problem_pt(problem_pt)
    {
     // Insist that all elements are of type 
     // NavierStokesElementWithDiagonalMassMatrices and issue a warning
     // if one is not.
     Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements=false;

     // Use Robin BC elements for Fp preconditioner -- yes by default
     Use_robin_for_fp=true;

     // 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 the LSC version
     Use_LSC=true;

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

     // Pin pressure dof in press adv diff problem for Fp precond
     Pin_first_pressure_dof_in_press_adv_diff=true;

     Navier_stokes_mesh_pt = 0;

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

     // 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 in Fp
     E_mat_vec_pt = 0;

     // Initially assume that there are no multiple element types in the 
     // Navier-Stokes mesh
     Allow_multiple_element_type_in_navier_stokes_mesh = false;
    }

   /// Destructor
   virtual ~NavierStokesSchurComplementPreconditioner()
    {
     clean_up_memory();
    }

   /// Broken copy constructor
   NavierStokesSchurComplementPreconditioner(
    const NavierStokesSchurComplementPreconditioner&)
    {
     BrokenCopy::broken_copy("NavierStokesSchurComplementPreconditioner");
    }
   
   /// 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 NavierStokesSchurComplementPreconditioner&)
    {
     BrokenCopy::broken_assign("NavierStokesSchurComplementPreconditioner");
     }*/

   /// Set the problem pointer (non-const pointer, the problem WILL be
   /// changed) for use in get_pressure_advection_diffusion_matrix().
   void set_problem_pt(Problem* problem_pt)
   {Problem_pt = problem_pt;}

  Problem* problem_pt() const
  {
#ifdef PARANOID
   if(Problem_pt == 0)
    {
     std::ostringstream error_msg;
     error_msg << "Problem pointer is null, did you set it yet?";
     throw OomphLibError(error_msg.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return Problem_pt;
  }
  
  /// \short Accept presence of elements that are not of type
  /// NavierStokesElementWithDiagonalMassMatrices without issuing a warning.
  void enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
  {
   Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements=true;
  }


  /// \short Do not accept presence of elements that are not of type
  /// NavierStokesElementWithDiagonalMassMatrices without issuing a warning.
  void disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
  {
   Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements=false;
  }


  /// Setup the preconditioner
  void setup();
  
  /// \short for some reason we have to remind the compiler that there is a
  /// setup() function in Preconditioner base class.
  using Preconditioner::setup;
  
   /// Apply preconditioner to Vector r
   void preconditioner_solve(const DoubleVector&r, DoubleVector &z);
   
   /// Specify the mesh containing the block-preconditionable Navier-Stokes
   /// elements. The optional argument indicates if there are multiple types
   /// of elements in the same mesh.
   void set_navier_stokes_mesh
     (Mesh* mesh_pt, 
      const bool& allow_multiple_element_type_in_navier_stokes_mesh = false)
    { 
     // Store the mesh pointer.
     Navier_stokes_mesh_pt = mesh_pt;
     
     // Are there multiple element types in this mesh?
     Allow_multiple_element_type_in_navier_stokes_mesh 
       = allow_multiple_element_type_in_navier_stokes_mesh;
    }

   /// 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;
   }

   /// Use LSC version of the preconditioner
   void use_lsc(){Use_LSC=true;}
   
   /// Use Fp version of the preconditioner
   void use_fp(){Use_LSC=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 Helper function to delete preconditioner data.
   void clean_up_memory();

   /// \short Use  Robin BC elements for the Fp preconditioner
   void enable_robin_for_fp() {Use_robin_for_fp = true;}

   /// \short Don't use Robin BC elements for the Fp preconditioenr
   void disable_robin_for_fp() {Use_robin_for_fp = false;}

   /// \short Set boolean indicating that we want to pin first pressure 
   /// dof in Navier Stokes mesh when
   /// assembling the pressure advection diffusion system for
   /// Fp preconditoner -- needed at zero Reynolds number
   /// for non-enclosed flows but seems harmless in any case
   void pin_first_pressure_dof_in_press_adv_diff()
   {Pin_first_pressure_dof_in_press_adv_diff=true;}

    /// \short Set boolean indicating that we do not want to pin first pressure 
   /// dof in Navier Stokes mesh when
   /// assembling the pressure advection diffusion system for
   /// Fp preconditoner -- needed at zero Reynolds number
   /// for non-enclosed flows but seems harmless in any case
   void unpin_first_pressure_dof_in_press_adv_diff()
   {Pin_first_pressure_dof_in_press_adv_diff=false;}

   /// \short Validate auxiliary pressure advection diffusion problem
   /// in 2D
   template<class ELEMENT>
   void validate(DocInfo& doc_info, Problem* orig_problem_pt)
   {
    FpPressureAdvectionDiffusionProblem<ELEMENT> 
     aux_problem(Navier_stokes_mesh_pt,orig_problem_pt);
    aux_problem.validate(doc_info); 
   }
   
   /// \short Pin all non-pressure dofs
   void pin_all_non_pressure_dofs()
   {
    // Backup pin status and then pin all non-pressure degrees of freedom
    unsigned nelem=Navier_stokes_mesh_pt->nelement();
    for (unsigned e=0;e<nelem;e++)
     {
      // Get pointer to the bulk element that is adjacent to boundary b
      TemplateFreeNavierStokesEquationsBase* el_pt = 
       dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
        Navier_stokes_mesh_pt->element_pt(e));
      
      // Pin 
      el_pt->pin_all_non_pressure_dofs(Eqn_number_backup);
     }
   }


   /// Get the pressure advection diffusion matrix
   void get_pressure_advection_diffusion_matrix(CRDoubleMatrix& fp_matrix)
   {
    // Pin all non-pressure dofs 
    pin_all_non_pressure_dofs(); 
    
    // Get the spatial dimension
    unsigned ndim=Navier_stokes_mesh_pt->finite_element_pt(0)->dim(); 

    // Backup assembly handler and set new one
    AssemblyHandler* backed_up_assembly_handler_pt=
     problem_pt()->assembly_handler_pt();
    problem_pt()->assembly_handler_pt()=
     new FpPreconditionerAssemblyHandler(ndim);

    // Backup submeshes (if any)
    unsigned n_sub_mesh=problem_pt()->nsub_mesh();
    Vector<Mesh*> backed_up_sub_mesh_pt(n_sub_mesh);
    for (unsigned i=0;i<n_sub_mesh;i++)
     {
      backed_up_sub_mesh_pt[i]=problem_pt()->mesh_pt(i);
     }
    // Flush sub-meshes but don't call rebuild_global_mesh()
    // so we can simply re-instate the sub-meshes below.
    problem_pt()->flush_sub_meshes();
    
    // Backup the problem's own mesh pointer 
    Mesh* backed_up_mesh_pt=problem_pt()->mesh_pt();
    problem_pt()->mesh_pt()=Navier_stokes_mesh_pt;

#ifdef OOMPH_HAS_MPI

    // Backup the start and end elements for the distributed
    // assembly process
    Vector<unsigned> backed_up_first_el_for_assembly;
    Vector<unsigned> backed_up_last_el_for_assembly;
    problem_pt()->get_first_and_last_element_for_assembly(
     backed_up_first_el_for_assembly,
     backed_up_last_el_for_assembly);

    // Now re-assign
    problem_pt()->set_default_first_and_last_element_for_assembly();

#endif

    // Identify pinned pressure dof
    int pinned_pressure_eqn=-2;
    if (Pin_first_pressure_dof_in_press_adv_diff)
     {
      // Loop over all Navier Stokes elements to find first non-pinned 
      // pressure dof
      unsigned nel=Navier_stokes_mesh_pt->nelement();
      for (unsigned e=0;e<nel;e++)
       {
        TemplateFreeNavierStokesEquationsBase* bulk_elem_pt = 
         dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
          Navier_stokes_mesh_pt->element_pt(e));
        int local_eqn=bulk_elem_pt->p_local_eqn(0);
        if (local_eqn>=0)
         {
          pinned_pressure_eqn=bulk_elem_pt->eqn_number(local_eqn);
          break;
         }
       }



#ifdef OOMPH_HAS_MPI
      if (problem_pt()->distributed())
       {
        int pinned_pressure_eqn_local=pinned_pressure_eqn;
        int pinned_pressure_eqn_global=pinned_pressure_eqn;
        MPI_Allreduce(&pinned_pressure_eqn_local,
                      &pinned_pressure_eqn_global,1,MPI_INT,MPI_MIN,
                      this->comm_pt()->mpi_comm());
        pinned_pressure_eqn=pinned_pressure_eqn_global;
       }    
      
#endif
      
     }
    
    
    // Loop over all Navier Stokes elements
    unsigned nel=Navier_stokes_mesh_pt->nelement();
    for (unsigned e=0;e<nel;e++)
     {
      TemplateFreeNavierStokesEquationsBase* bulk_elem_pt = 
       dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
        Navier_stokes_mesh_pt->element_pt(e));
      
      // Set pinned pressure equation
      bulk_elem_pt->pinned_fp_pressure_eqn()=pinned_pressure_eqn;      
     }
    
 
    // Attach Robin BC elements
    unsigned nbound=Navier_stokes_mesh_pt->nboundary(); 
    if (Use_robin_for_fp)
     {
      // Loop over all boundaries of Navier Stokes mesh
      for (unsigned b=0;b<nbound;b++)
       {
        // How many bulk elements are adjacent to boundary b?
        unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
        
        // Loop over the bulk elements adjacent to boundary b?
        for(unsigned e=0;e<n_element;e++)
         {
          TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
           dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
            Navier_stokes_mesh_pt->boundary_element_pt(b,e));
          
          //What is the index of the face of the bulk element e on bondary b
          int face_index=Navier_stokes_mesh_pt->face_index_at_boundary(b,e);
          
          // Build face element
          bulk_elem_pt->build_fp_press_adv_diff_robin_bc_element(face_index);
          
         } //end of loop over bulk elements adjacent to boundary b
       }
     }
        
    // Get "Jacobian" of the modified system
    DoubleVector dummy_residuals;
    problem_pt()->get_jacobian(dummy_residuals,fp_matrix);

    // Kill Robin BC elements
    if (Use_robin_for_fp)
     {
      // Loop over all boundaries of Navier Stokes mesh
      for (unsigned b=0;b<nbound;b++)
       {
        // How many bulk elements are adjacent to boundary b?
        unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
        
        // Loop over the bulk elements adjacent to boundary b?
        for(unsigned e=0;e<n_element;e++)
         {
          
          TemplateFreeNavierStokesEquationsBase* bulk_elem_pt = 
           dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
            Navier_stokes_mesh_pt->boundary_element_pt(b,e));
          
          // Kill associated Robin elements
          bulk_elem_pt->delete_pressure_advection_diffusion_robin_elements();

         } //end of loop over bulk elements adjacent to boundary b
       }
     }

    // Reset pin status
    reset_pin_status();

#ifdef OOMPH_HAS_MPI

    // Reset start and end elements for the distributed
    // assembly process
    problem_pt()->set_first_and_last_element_for_assembly(
     backed_up_first_el_for_assembly,
     backed_up_last_el_for_assembly);
    
#endif

    // Cleanup and reset assembly handler
    delete problem_pt()->assembly_handler_pt();
    problem_pt()->assembly_handler_pt()=backed_up_assembly_handler_pt;

    // Re-instate submeshes. (No need to call rebuild_global_mesh()
    // as it was never unbuilt).
    for (unsigned i=0;i<n_sub_mesh;i++)
     {
      problem_pt()->add_sub_mesh(backed_up_sub_mesh_pt[i]);
     }

    
    // Reset the problem's mesh pointer 
    problem_pt()->mesh_pt()= backed_up_mesh_pt;    
   }

   
   /// \short Reset pin status of all values
   void reset_pin_status()
   {

    // Reset pin status for nodes
    unsigned nnod=Navier_stokes_mesh_pt->nnode();
    for (unsigned j=0;j<nnod;j++)
     {
      Node* nod_pt=Navier_stokes_mesh_pt->node_pt(j);
      unsigned nval=nod_pt->nvalue();
      for (unsigned i=0;i<nval;i++)
       {
        nod_pt->eqn_number(i)=Eqn_number_backup[nod_pt][i];
       }

      // If it's a solid node deal with its positional data too
      SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
      if (solid_nod_pt!=0)
       {
        Data* solid_posn_data_pt=solid_nod_pt->variable_position_pt();
        unsigned nval=solid_posn_data_pt->nvalue();
        for (unsigned i=0;i<nval;i++)
         {
          solid_posn_data_pt->eqn_number(i)=
           Eqn_number_backup[solid_posn_data_pt][i];
         }
       }
      
     }
    
    // ... and internal data
    unsigned nelem=Navier_stokes_mesh_pt->nelement();
    for (unsigned e=0;e<nelem;e++)
     {
      // Pointer to element
      GeneralisedElement* el_pt=Navier_stokes_mesh_pt->element_pt(e);
      
      // Pin/unpin internal data
      unsigned nint=el_pt->ninternal_data();
      for (unsigned j=0;j<nint;j++)
       {
        Data* data_pt=el_pt->internal_data_pt(j);
        unsigned nvalue=data_pt->nvalue();
        for (unsigned i=0;i<nvalue;i++)
         {
          data_pt->eqn_number(i)=Eqn_number_backup[data_pt][i];
         }
       }
     }
    
    // Free up storage
    Eqn_number_backup.clear();
   }
  
    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 Boolean to indicate that presence of elements that are not of type
   /// NavierStokesElementWithDiagonalMassMatrices is acceptable (suppresses
   /// warning that issued otherwise).
   bool Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements;

   /// \short Helper function to assemble the inverse diagonals of the pressure
   /// and velocity mass matrices from the elemental contributions defined in
   /// NavierStokesEquations<DIM>.
   /// If do_both=true, both are computed, otherwise only the velocity
   /// mass matrix (the LSC version of the preconditioner only needs
   /// that one)
   void assemble_inv_press_and_veloc_mass_matrix_diagonal(
    CRDoubleMatrix*& inv_p_mass_pt, 
    CRDoubleMatrix*& inv_v_mass_pt, 
    const bool& do_both);
   
   /// \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 Qv^{-1} Bt 
   MatrixVectorProduct* QBt_mat_vec_pt;

   /// MatrixVectorProduct operator for Bt
   MatrixVectorProduct* Bt_mat_vec_pt;

   /// MatrixVectorProduct operator for F
   MatrixVectorProduct* F_mat_vec_pt;

   /// MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
   MatrixVectorProduct* E_mat_vec_pt;

   /// \short the pointer to the mesh of block preconditionable Navier
   /// Stokes elements.
   Mesh* Navier_stokes_mesh_pt;
   
   /// \short Flag to indicate if there are multiple element types in the
   /// Navier-Stokes mesh.
   bool Allow_multiple_element_type_in_navier_stokes_mesh;

   /// Boolean to indicate use of LSC (true) or Fp (false) variant
   bool Use_LSC;
   
   /// Use Robin BC elements for Fp preconditoner?
   bool Use_robin_for_fp;

   /// \short Map to store original eqn numbers of various Data values when 
   /// assembling pressure advection diffusion matrix
   std::map<Data*,std::vector<int> > Eqn_number_backup;

   /// \short Boolean indicating if we want to pin first pressure 
   /// dof in Navier Stokes mesh when
   /// assembling the pressure advection diffusion system for
   /// Fp preconditoner -- needed at zero Reynolds number
   /// for non-enclosed flows but seems harmless in any case
   bool Pin_first_pressure_dof_in_press_adv_diff;
   
   /// Storage for the (non-const!) problem pointer for use in
   /// get_pressure_advection_diffusion_matrix().
   Problem* Problem_pt;

   };




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



//============================================================================
/// \short The exact Navier Stokes preconditioner. This extracts 2x2 blocks
/// (corresponding to the velocity 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 NavierStokesExactPreconditioner : public BlockPreconditioner<MATRIX>
  {

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

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

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