lagrange_enforced_flow_preconditioner.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    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_LAGRANGE_ENFORCED_FLOW_PRECONDITIONERS_HEADER
#define OOMPH_LAGRANGE_ENFORCED_FLOW_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 "../generic/general_purpose_preconditioners.h"
#include "../generic/general_purpose_block_preconditioners.h"
#ifdef OOMPH_HAS_HYPRE
#include "../generic/hypre_solver.h"
#endif
#ifdef OOMPH_HAS_TRILINOS
#include "../generic/trilinos_solver.h"
#endif

namespace oomph
{
//==========================================================================
/// Namespace for subsidiary preconditioner creation helper functions
//==========================================================================
namespace Lagrange_Enforced_Flow_Preconditioner_Subsidiary_Operator_Helper
{
  /// \short CG with diagonal preconditioner for W-block subsidiary linear 
  /// systems.
  extern Preconditioner* get_w_cg_preconditioner();

  /// \short Hypre Boomer AMG setting for the augmented momentum block
  /// of a 2D Navier-Stokes problem using the simple form of the viscous
  /// term (for serial code).
  extern Preconditioner* boomer_amg_for_2D_momentum_simple_visc();

  /// \short Hypre Boomer AMG setting for the augmented momentum block
  /// of a 2D Navier-Stokes problem using the stress divergence form of the 
  /// viscous term (for serial code).
  extern Preconditioner* boomer_amg_for_2D_momentum_stressdiv_visc();

  /// \short Hypre Boomer AMG setting for the augmented momentum block
  /// of a 3D Navier-Stokes problem (for serial code).
  extern Preconditioner* boomer_amg_for_3D_momentum();

  /// \short Hypre Boomer AMG setting for the augmented momentum block
  /// of a 3D Navier-Stokes problem (for serial code).
  extern Preconditioner* boomer_amg2v22_for_3D_momentum();

  /// \short Hypre Boomer AMG setting for the 2D Poisson problem 
  /// (for serial code).  
  extern Preconditioner* boomer_amg_for_2D_poisson_problem();

  /// \short Hypre Boomer AMG setting for the 3D Poisson problem 
  /// (for serial code).
  extern Preconditioner* boomer_amg_for_3D_poisson_problem();

} // namespace


//==========================================================================
/// \short The preconditioner for the Lagrange multiplier constrained 
/// Navier-Stokes equations. The velocity components are constrained by 
/// Lagrange multiplier, which are applied via OOMPH-LIB's FACE elements.
///
///
/// The linearised Jacobian takes the block form:
/// 
/// | F_ns | L^T |
/// |------------|
/// |   L  | 0   |
///
/// where L correspond to the constrained block,
/// F_ns is the Navier-Stokes block with the following block structure
///
/// |  F | G^T |
/// |----------|
/// |  D |  0  |
///
/// Here F is the momentum block, G the discrete gradient operator,
/// and D the discrete divergence operator. For unstabilised elements, 
/// we have D = G^T and in much of the literature the divergence matrix is 
/// denoted by B.
///
/// The Lagrange enforced flow preconditioner takes the form:
/// | F_aug |    |
/// |-------|----|
/// |       | Wd |
///
/// where  F_aug = F_ns + L^T*inv(Wd)*L is an augmented Navier-Stokes block 
/// and Wd=(1/Scaling_sigma)*diag(LL^T).
///
/// In our implementation of the preconditioner, the linear systems 
/// associated with the (1,1) block 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
///
/// LagrangeEnforcedFlowPreconditioner::set_navier_stokes_preconditioner(...)
///
/// Access to the elements is provided via meshes. However, a Vector of 
/// meshes is taken, each mesh contains a different type of block 
/// preconditionable element. This allows the (re-)classification of the 
/// constrained velocity DOF types.
/// 
/// The first mesh in the Vector Mesh_pt must be the 'bulk' mesh.
/// The rest are assumed to contain FACEELMENTS. 
/// 
/// Thus, the most general block structure (in 3D) is: 
/// 
///  0 1 2 3   4 5 6 7  8  ..x   x+0 x+1 x+2 x+3 x+4 
/// [u v w p] [u v w l1 l2 ...] [u   v   w   l1  l2 ...] ... 
///   Bulk       Surface 1             Surface 2         ... 
/// 
/// where the DOF types in [] are the DOF types associated with each mesh.
/// 
/// For example, consider a unit cube domain [0,1]^3 with parallel outflow 
/// imposed (in mesh 0) and tangential flow imposed (in mesh 1), then there 
/// are 13 DOF types and our implementation respects the following 
/// (natural) DOF type order:
///
///   bulk          mesh 0           mesh 1
/// [0 1 2 3] [4  5  6   7   8 ] [9  10 11 12 ]
/// [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1]
///
/// Via the appropriate mapping, the block_setup(...) function will
/// re-order the DOF types into the following block types:
/// 
///  0 1  2  3 4  5   6 7  8     9   10  11  12   <- Block type
///  0 4  9  1 5  10  2 6  11    3    7   8  12   <- DOF type
/// [u up ut v vp vt  w wp wt ] [p] [Lp1 Lp2 Lt1] 
/// 
//==========================================================================
class LagrangeEnforcedFlowPreconditioner 
: public BlockPreconditioner<CRDoubleMatrix>
{
public:

  /// \short This preconditioner includes the option to use subsidiary 
  /// operators other than SuperLUPreconditioner for this problem. 
  /// This is the typedef of a function that should return an instance
  /// of a subsidiary preconditioning operator.  This preconditioner is 
  /// responsible for the destruction of the subsidiary preconditioners.
  typedef Preconditioner* (*SubsidiaryPreconditionerFctPt)();

  /// Constructor - initialise variables.
  LagrangeEnforcedFlowPreconditioner()
    : BlockPreconditioner<CRDoubleMatrix>()
  {
    // The null pointer.
    Navier_stokes_preconditioner_pt = 0;

    // By default, the linear systems associated with the diagonal blocks 
    // are solved "exactly" using SuperLU (in its incarnation as an exact 
    // preconditioner. This is not a block preconditioner. 
    Navier_stokes_preconditioner_is_block_preconditioner = false;

    // Flag to indicate to use SuperLU or not.
    Using_superlu_ns_preconditioner = true;

    // Empty vector of meshes and set the number of meshes to zero.
    My_mesh_pt.resize(0,0);
    My_nmesh = 0;

    // The number of DOF types within the meshes.
    My_ndof_types_in_mesh.resize(0,0);

    // Initialise other variables.
    Use_norm_f_for_scaling_sigma = true;
    Scaling_sigma = 0.0;
    N_lagrange_doftypes = 0;
    N_fluid_doftypes = 0;
    N_velocity_doftypes = 0;
    Preconditioner_has_been_setup = false;
  } // constructor

  /// \short Destructor
  virtual ~LagrangeEnforcedFlowPreconditioner()
  {
    this->clean_up_memory();
  }

  /// \short Broken copy constructor
  LagrangeEnforcedFlowPreconditioner 
    (const LagrangeEnforcedFlowPreconditioner&)
    {
      BrokenCopy::broken_copy("LagrangeEnforcedFlowPreconditioner");
    }

  /// \short Broken assignment operator
  void operator=(const LagrangeEnforcedFlowPreconditioner&)
  {
    BrokenCopy::broken_assign(" LagrangeEnforcedFlowPreconditioner");
  }

  /// \short Setup method for the LagrangeEnforcedFlowPreconditioner.
  void setup();

  /// \short Apply the preconditioner.
  /// r is the residual (rhs), z will contain the solution.
  void preconditioner_solve(const DoubleVector& r, DoubleVector& z);

  /// \short Set the meshes,
  /// the first mesh in the vector must be the bulk mesh.
  void set_meshes(const Vector<Mesh*> &mesh_pt);

  /// \short Set flag to use the infinite norm of the Navier-Stokes F matrix
  /// as the scaling sigma. This is the default behaviour. Note: the norm of
  /// the NS F matrix positive, however, we actually use the negative of 
  /// the norm. This is because the underlying Navier-Stokes Jacobian is 
  /// multiplied by -1. Ask Andrew/Matthias for more detail.
  void use_norm_f_for_scaling_sigma()
  {
    Use_norm_f_for_scaling_sigma = true;
  }

  /// \short Access function to set the scaling sigma. 
  /// Note: this also sets the flag to use the infinite norm of 
  /// the Navier-Stokes F matrix as the scaling sigma to false.
  /// Warning is given if trying to set scaling sigma to be equal to 
  /// or greater than zero.
  void set_scaling_sigma(const double& scaling_sigma)
  {
    // Check if scaling sigma is zero or positive.
#ifdef PARANOID
    if(scaling_sigma == 0.0)
    {
      std::ostringstream warning_stream;
      warning_stream << "WARNING: \n"
        << "Setting scaling_sigma = 0.0 may cause values.\n";
      OomphLibWarning(warning_stream.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
    }
    if(scaling_sigma > 0.0)
    {
      std::ostringstream warning_stream;
      warning_stream << "WARNING: " << std::endl
        << "The scaling (scaling_sigma) is positive: " 
        << Scaling_sigma << "\n"
        << "Performance may be degraded.\n";
      OomphLibWarning(warning_stream.str(),
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
    }
#endif

    Scaling_sigma = scaling_sigma;
    Use_norm_f_for_scaling_sigma = false;
  }

  /// \short Read (const) function to get the scaling sigma.
  double scaling_sigma() const
  {
    return Scaling_sigma;
  }

  /// \short Set a new Navier-Stokes matrix preconditioner 
  /// (inexact solver)
  void set_navier_stokes_preconditioner(
      Preconditioner* new_ns_preconditioner_pt = 0);

  ///\short Set Navier-Stokes matrix preconditioner (inexact
  /// solver) to SuperLU
  void set_superlu_for_navier_stokes_preconditioner()
  {
    if (!Using_superlu_ns_preconditioner)
    {
      delete Navier_stokes_preconditioner_pt;
      Navier_stokes_preconditioner_pt = new SuperLUPreconditioner;
      Using_superlu_ns_preconditioner = true;
    }
  }

  /// \short Clears the memory.
  void clean_up_memory();

private:

  /// \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 Scaling for the augmentation: Scaling_sigma*(LL^T)
  double Scaling_sigma;

  /// \short Flag to indicate if we want to use the infinite norm of the
  /// Navier-Stokes momentum block for the scaling sigma.
  bool Use_norm_f_for_scaling_sigma;

  /// \short Inverse W values
  Vector<Vector<double> > Inv_w_diag_values;

  /// \short Pointer to the 'preconditioner' for the Navier-Stokes block
  Preconditioner* Navier_stokes_preconditioner_pt;

  /// \short Flag to indicate if the preconditioner for the Navier-Stokes 
  /// block is a block preconditioner or not.
  bool Navier_stokes_preconditioner_is_block_preconditioner;

  /// \short Flag to indicate whether the default NS preconditioner is used
  bool Using_superlu_ns_preconditioner;

  /// \short Storage for the meshes. In our implementation, the first mesh 
  /// must always be the Navier-Stokes (bulk) mesh, followed by surface 
  /// meshes.
  Vector<Mesh*> My_mesh_pt;

  /// \short The number of DOF types in each mesh. This is used create 
  /// various lookup lists.
  Vector<unsigned> My_ndof_types_in_mesh;

  /// \short The number of meshes. This is used to create various lookup 
  /// lists.
  unsigned My_nmesh;

  /// \short The number of Lagrange multiplier DOF types.
  unsigned N_lagrange_doftypes;

  /// \short The number of fluid DOF types (including pressure).
  unsigned N_fluid_doftypes;

  /// \short The number of velocity DOF types.
  unsigned N_velocity_doftypes;

}; // end of LagrangeEnforcedFlowPreconditioner class

}
#endif