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-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_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 PressureAdvectionDiffusionValidation
  {
    extern unsigned Flag;
 
    extern double Peclet;
 
    extern void wind_function(const Vector<double>& x, Vector<double>& wind);
 
    extern void get_exact_u(const Vector<double>& x, Vector<double>& u);
 
    extern void get_exact_u(const Vector<double>& x, double& u);
 
    extern double source_function(const Vector<double>& x_vect);
 
  } // namespace PressureAdvectionDiffusionValidation
 
 
 
 
  //=============================================================
  //===============================================================
  class FpPreconditionerAssemblyHandler : public AssemblyHandler
  {
  public:
    FpPreconditionerAssemblyHandler(const unsigned& ndim) // : Ndim(ndim)
    {
    }
 
    virtual ~FpPreconditionerAssemblyHandler() {}
 
    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);
    }
 
    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:
    // unsigned Ndim; // cgj: Ndim is never used.
  };
 
 
 
 
  //===========================================================================
  //===========================================================================
  template<class ELEMENT>
  class FpPressureAdvectionDiffusionProblem : public virtual Problem
  {
  public:
    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);
    }
 
 
    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;
    }
 
 
    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();
    }
 
 
    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]);
          }
        }
      }
    }
 
 
    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";
    }
 
 
    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:
    unsigned Ndim;
 
    Problem* Orig_problem_pt;
 
    std::map<Data*, std::vector<int>> Eqn_number_backup;
  };
 
 
 
 
  //===========================================================================
  //===========================================================================
  class NavierStokesSchurComplementPreconditioner
    : public BlockPreconditioner<CRDoubleMatrix>
  {
  public:
    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;
    }
 
    virtual ~NavierStokesSchurComplementPreconditioner()
    {
      clean_up_memory();
    }
 
    NavierStokesSchurComplementPreconditioner(
      const NavierStokesSchurComplementPreconditioner&) = delete;
 
    // Commented out broken assignment operator because this can lead to a
    // conflict warning when used in the virtual inheritence hierarchy.
    // Essentially the compiler doesn't realise that two separate
    // implementations of the broken function are the same and so, quite
    // rightly, it shouts.
    /*void operator=(const NavierStokesSchurComplementPreconditioner&) =
     * delete;*/
 
    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;
    }
 
    void enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
    {
      Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements = true;
    }
 
 
    void disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()
    {
      Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements = false;
    }
 
 
    void setup();
 
    using Preconditioner::setup;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    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;
    }
 
    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;
    }
 
    void set_p_superlu_preconditioner()
    {
      if (!Using_default_p_preconditioner)
      {
        P_preconditioner_pt = new SuperLUPreconditioner;
        Using_default_p_preconditioner = true;
      }
    }
 
    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;
    }
 
    void use_lsc()
    {
      Use_LSC = true;
    }
 
    void use_fp()
    {
      Use_LSC = false;
    }
 
    void set_f_superlu_preconditioner()
    {
      if (!Using_default_f_preconditioner)
      {
        F_preconditioner_pt = new SuperLUPreconditioner;
        Using_default_f_preconditioner = true;
      }
    }
 
    void enable_doc_time()
    {
      Doc_time = true;
    }
 
    void disable_doc_time()
    {
      Doc_time = false;
    }
 
    void clean_up_memory();
 
    void enable_robin_for_fp()
    {
      Use_robin_for_fp = true;
    }
 
    void disable_robin_for_fp()
    {
      Use_robin_for_fp = false;
    }
 
    void pin_first_pressure_dof_in_press_adv_diff()
    {
      Pin_first_pressure_dof_in_press_adv_diff = true;
    }
 
    void unpin_first_pressure_dof_in_press_adv_diff()
    {
      Pin_first_pressure_dof_in_press_adv_diff = false;
    }
 
    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);
    }
 
    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);
      }
    }
 
 
    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;
    }
 
 
    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
    // -----------------------------------------------------
    Preconditioner* P_preconditioner_pt;
 
    Preconditioner* F_preconditioner_pt;
 
    bool Using_default_f_preconditioner;
 
    bool Using_default_p_preconditioner;
 
    bool Preconditioner_has_been_setup;
 
    bool Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements;
 
    void assemble_inv_press_and_veloc_mass_matrix_diagonal(
      CRDoubleMatrix*& inv_p_mass_pt,
      CRDoubleMatrix*& inv_v_mass_pt,
      const bool& do_both);
 
    bool F_preconditioner_is_block_preconditioner;
 
    bool Doc_time;
 
    MatrixVectorProduct* QBt_mat_vec_pt;
 
    MatrixVectorProduct* Bt_mat_vec_pt;
 
    MatrixVectorProduct* F_mat_vec_pt;
 
    MatrixVectorProduct* E_mat_vec_pt;
 
    Mesh* Navier_stokes_mesh_pt;
 
    bool Allow_multiple_element_type_in_navier_stokes_mesh;
 
    bool Use_LSC;
 
    bool Use_robin_for_fp;
 
    std::map<Data*, std::vector<int>> Eqn_number_backup;
 
    bool Pin_first_pressure_dof_in_press_adv_diff;
 
    Problem* Problem_pt;
  };
 
 
 
 
  //============================================================================
  //=============================================================================
  template<typename MATRIX>
  class NavierStokesExactPreconditioner : public BlockPreconditioner<MATRIX>
  {
  public:
    NavierStokesExactPreconditioner() : BlockPreconditioner<MATRIX>() {}
 
 
    ~NavierStokesExactPreconditioner() {}
 
 
    NavierStokesExactPreconditioner(const NavierStokesExactPreconditioner&) =
      delete;
 
    /*void operator=(const NavierStokesExactPreconditioner&) = delete;*/
 
    void setup();
 
 
    void preconditioner_solve(const Vector<double>& r, Vector<double>& z);
 
  protected:
    MATRIX P_matrix;
  };
 
} // namespace oomph
#endif