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// Copyright (C) 2006-2024 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);
 
  } // namespace PressureAdvectionDiffusionValidation
 
 
  /// /////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////
 
 
  //=============================================================
  /// 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)
    {
    }
 
    /// Empty virtual destructor
    virtual ~FpPreconditionerAssemblyHandler() {}
 
    /// 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);
    }
 
    /// 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.
  };
 
 
  /// ////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// 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;
    }
 
 
    /// 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();
    }
 
 
    /// 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]);
          }
        }
      }
    }
 
 
    /// 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";
    }
 
 
    /// 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;
 
    /// 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;
  };
 
 
  /// /////////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////////
  /// /////////////////////////////////////////////////////////////////////////
 
 
  //===========================================================================
  /// 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&) = delete;
 
    /// 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&) =
     * delete;*/
 
    /// 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;
    }
 
    /// 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;
    }
 
 
    /// 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();
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// Helper function to delete preconditioner data.
    void clean_up_memory();
 
    /// Use  Robin BC elements for the Fp preconditioner
    void enable_robin_for_fp()
    {
      Use_robin_for_fp = true;
    }
 
    /// Don't use Robin BC elements for the Fp preconditioenr
    void disable_robin_for_fp()
    {
      Use_robin_for_fp = false;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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);
    }
 
    /// 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;
    }
 
 
    /// 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;
 
    /// 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;
 
    /// 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;
 
    /// 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);
 
    /// 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;
 
    /// the pointer to the mesh of block preconditionable Navier
    /// Stokes elements.
    Mesh* Navier_stokes_mesh_pt;
 
    /// 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;
 
    /// 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;
 
    /// 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;
  };
 
 
  /// ////////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////////
  /// ////////////////////////////////////////////////////////////////////////////
 
 
  //============================================================================
  /// 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&) =
      delete;
 
    /// Broken assignment operator
    /*void operator=(const NavierStokesExactPreconditioner&) = delete;*/
 
    /// 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;
  };
 
} // namespace oomph
#endif