#include <lagrange_enforced_flow_preconditioner.h>

Inheritance diagram for oomph::LagrangeEnforcedFlowPreconditioner:

Public Types
typedef Preconditioner (	SubsidiaryPreconditionerFctPt) ()

Public Member Functions
	LagrangeEnforcedFlowPreconditioner ()
	Constructor - initialise variables. More...

virtual	~LagrangeEnforcedFlowPreconditioner ()
	Destructor. More...

	LagrangeEnforcedFlowPreconditioner (const LagrangeEnforcedFlowPreconditioner &)=delete
	Broken copy constructor. More...

void	operator= (const LagrangeEnforcedFlowPreconditioner &)=delete
	Broken assignment operator. More...

void	setup ()
	Setup method for the LagrangeEnforcedFlowPreconditioner. More...

void	preconditioner_solve (const DoubleVector &r, DoubleVector &z)

void	set_meshes (const Vector< Mesh * > &mesh_pt)

void	use_norm_f_for_scaling_sigma ()

void	set_scaling_sigma (const double &scaling_sigma)

double	scaling_sigma () const
	Read (const) function to get the scaling sigma. More...

void	set_navier_stokes_preconditioner (Preconditioner *new_ns_preconditioner_pt=0)

void	set_superlu_for_navier_stokes_preconditioner ()

void	clean_up_memory ()
	Clears the memory. More...

Public Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
	BlockPreconditioner ()
	Constructor. More...

	BlockPreconditioner (const BlockPreconditioner &)=delete
	Broken copy constructor. More...

virtual	~BlockPreconditioner ()
	Destructor. More...

void	operator= (const BlockPreconditioner &)=delete
	Broken assignment operator. More...

CRDoubleMatrix *	matrix_pt () const

void	turn_on_recursive_debug_flag ()

void	turn_off_recursive_debug_flag ()

void	turn_on_debug_flag ()
	Toggles on the debug flag. More...

void	turn_off_debug_flag ()
	Toggles off the debug flag. More...

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

virtual void	block_setup ()

void	block_setup (const Vector< unsigned > &dof_to_block_map)

void	get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const

CRDoubleMatrix	get_block (const unsigned &i, const unsigned &j, const bool &ignore_replacement_block=false) const

void	set_master_matrix_pt (CRDoubleMatrix *in_matrix_pt)
	Set the matrix_pt in the upper-most master preconditioner. More...

void	get_block_other_matrix (const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)

void	get_blocks (DenseMatrix< bool > &required_blocks, DenseMatrix< CRDoubleMatrix * > &block_matrix_pt) const

void	get_dof_level_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block=false) const

void	get_dof_level_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block) const

CRDoubleMatrix	get_concatenated_block (const VectorMatrix< BlockSelector > &selected_block)

void	get_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)

void	return_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
	Takes concatenated block ordered vector, b, and copies its. More...

void	get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w)

void	return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	nblock_types () const
	Return the number of block types. More...

unsigned	ndof_types () const
	Return the total number of DOF types. More...

const Mesh *	mesh_pt (const unsigned &i) const

unsigned	nmesh () const

int	block_number (const unsigned &i_dof) const
	Return the block number corresponding to a global index i_dof. More...

int	index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	block_distribution_pt (const unsigned &b) const
	Access function to the block distributions (const version). More...

LinearAlgebraDistribution *	block_distribution_pt (const unsigned b)
	Access function to the block distributions (non-const version). More...

LinearAlgebraDistribution *	dof_block_distribution_pt (const unsigned &b)
	Access function to the dof-level block distributions. More...

const LinearAlgebraDistribution *	master_distribution_pt () const

unsigned	ndof_types_in_mesh (const unsigned &i) const

bool	is_subsidiary_block_preconditioner () const

bool	is_master_block_preconditioner () const

void	set_block_output_to_files (const std::string &basefilename)

void	disable_block_output_to_files ()
	Turn off output of blocks (by clearing the basefilename string). More...

bool	block_output_on () const
	Test if output of blocks is on or not. More...

void	output_blocks_to_files (const std::string &basefilename, const unsigned &precision=8) const

void	post_block_matrix_assembly_partial_clear ()

BlockPreconditioner< CRDoubleMatrix > *	master_block_preconditioner_pt () const
	Access function to the master block preconditioner pt. More...

void	clear_block_preconditioner_base ()

void	document ()

Vector< Vector< unsigned > >	doftype_coarsen_map_fine () const

Vector< unsigned >	get_fine_grain_dof_types_in (const unsigned &i) const

unsigned	nfine_grain_dof_types_in (const unsigned &i) const

MapMatrix< unsigned, CRDoubleMatrix * >	replacement_dof_block_pt () const
	Access function to the replaced dof-level blocks. More...

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const Vector< unsigned > &block_col_indices)

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const unsigned &block_col_index)

void	internal_get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w) const

void	internal_return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	internal_nblock_types () const

unsigned	internal_ndof_types () const

void	internal_return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	internal_get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	internal_get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block) const

void	internal_get_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block) const

int	internal_block_number (const unsigned &i_dof) const

int	internal_index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	internal_block_distribution_pt (const unsigned &b) const
	Access function to the internal block distributions. More...

void	insert_auxiliary_block_distribution (const Vector< unsigned > &block_vec_number, LinearAlgebraDistribution *dist_pt)

void	block_matrix_test (const unsigned &i, const unsigned &j, const CRDoubleMatrix *block_matrix_pt) const

int	get_index_of_value (const Vector< myType > &vec, const myType val, const bool sorted=false) const

Public Member Functions inherited from oomph::Preconditioner
	Preconditioner ()
	Constructor. More...

	Preconditioner (const Preconditioner &)=delete
	Broken copy constructor. More...

void	operator= (const Preconditioner &)=delete
	Broken assignment operator. More...

virtual	~Preconditioner ()
	Destructor (empty) More...

virtual void	preconditioner_solve_transpose (const DoubleVector &r, DoubleVector &z)

void	setup (DoubleMatrixBase *matrix_pt)

void	setup (const Problem problem_pt, DoubleMatrixBase matrix_pt)

void	enable_silent_preconditioner_setup ()
	Set up the block preconditioner quietly! More...

void	disable_silent_preconditioner_setup ()
	Be verbose in the block preconditioner setup. More...

virtual void	set_matrix_pt (DoubleMatrixBase *matrix_pt)
	Set the matrix pointer. More...

virtual const OomphCommunicator *	comm_pt () const
	Get function for comm pointer. More...

virtual void	set_comm_pt (const OomphCommunicator *const comm_pt)
	Set the communicator pointer. More...

double	setup_time () const
	Returns the time to setup the preconditioner. More...

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...

	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...

void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...

virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...

LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...

unsigned	nrow () const
	access function to the number of global rows. More...

unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...

unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...

unsigned	first_row () const
	access function for the first row on this processor More...

unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...

bool	distributed () const
	distribution is serial or distributed More...

bool	distribution_built () const

void	build_distribution (const LinearAlgebraDistribution *const dist_pt)

void	build_distribution (const LinearAlgebraDistribution &dist)

Private Attributes
bool	Preconditioner_has_been_setup

double	Scaling_sigma
	Scaling for the augmentation: Scaling_sigma*(LL^T) More...

bool	Use_norm_f_for_scaling_sigma

Vector< Vector< double > >	Inv_w_diag_values
	Inverse W values. More...

Preconditioner *	Navier_stokes_preconditioner_pt
	Pointer to the 'preconditioner' for the Navier-Stokes block. More...

bool	Navier_stokes_preconditioner_is_block_preconditioner

bool	Using_superlu_ns_preconditioner
	Flag to indicate whether the default NS preconditioner is used. More...

Vector< Mesh * >	My_mesh_pt

Vector< unsigned >	My_ndof_types_in_mesh

unsigned	My_nmesh

unsigned	N_lagrange_doftypes
	The number of Lagrange multiplier DOF types. More...

unsigned	N_fluid_doftypes
	The number of fluid DOF types (including pressure). More...

unsigned	N_velocity_doftypes
	The number of velocity DOF types. More...

Additional Inherited Members
Protected Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
void	set_nmesh (const unsigned &n)

void	set_mesh (const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)

void	set_replacement_dof_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)

bool	any_mesh_distributed () const

int	internal_dof_number (const unsigned &i_dof) const

unsigned	internal_index_in_dof (const unsigned &i_dof) const

unsigned	internal_block_dimension (const unsigned &b) const

unsigned	internal_dof_block_dimension (const unsigned &i) const

unsigned	master_nrow () const

unsigned	internal_master_dof_number (const unsigned &b) const

const LinearAlgebraDistribution *	internal_preconditioner_matrix_distribution_pt () const

const LinearAlgebraDistribution *	preconditioner_matrix_distribution_pt () const

Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()

Protected Attributes inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
MapMatrix< unsigned, CRDoubleMatrix * >	Replacement_dof_block_pt
	The replacement dof-level blocks. More...

Vector< LinearAlgebraDistribution * >	Block_distribution_pt
	The distribution for the blocks. More...

Vector< Vector< unsigned > >	Block_to_dof_map_coarse

Vector< Vector< unsigned > >	Block_to_dof_map_fine
	Mapping for the block types to the most fine grain dof types. More...

Vector< Vector< unsigned > >	Doftype_coarsen_map_coarse

Vector< Vector< unsigned > >	Doftype_coarsen_map_fine

Vector< LinearAlgebraDistribution * >	Internal_block_distribution_pt
	Storage for the default distribution for each internal block. More...

Vector< LinearAlgebraDistribution * >	Dof_block_distribution_pt

Vector< unsigned >	Allow_multiple_element_type_in_mesh

Vector< const Mesh * >	Mesh_pt

Vector< unsigned >	Ndof_types_in_mesh

unsigned	Internal_nblock_types

unsigned	Internal_ndof_types

Protected Attributes inherited from oomph::Preconditioner
bool	Silent_preconditioner_setup
	Boolean to indicate whether or not the build should be done silently. More...

std::ostream *	Stream_pt
	Pointer to the output stream – defaults to std::cout. More...

Detailed Description

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]

Member Typedef Documentation

◆ SubsidiaryPreconditionerFctPt

typedef Preconditioner*(* oomph::LagrangeEnforcedFlowPreconditioner::SubsidiaryPreconditionerFctPt) ()

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.

Constructor & Destructor Documentation

◆ LagrangeEnforcedFlowPreconditioner() [1/2]

oomph::LagrangeEnforcedFlowPreconditioner::LagrangeEnforcedFlowPreconditioner ( )

inline

Constructor - initialise variables.

                                          : 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

References My_mesh_pt, My_ndof_types_in_mesh, My_nmesh, N_fluid_doftypes, N_lagrange_doftypes, N_velocity_doftypes, Navier_stokes_preconditioner_is_block_preconditioner, Navier_stokes_preconditioner_pt, Preconditioner_has_been_setup, Scaling_sigma, Use_norm_f_for_scaling_sigma, and Using_superlu_ns_preconditioner.

◆ ~LagrangeEnforcedFlowPreconditioner()

virtual oomph::LagrangeEnforcedFlowPreconditioner::~LagrangeEnforcedFlowPreconditioner ( )

inlinevirtual

Destructor.

     {
       this->clean_up_memory();
     }

References clean_up_memory().

◆ LagrangeEnforcedFlowPreconditioner() [2/2]

oomph::LagrangeEnforcedFlowPreconditioner::LagrangeEnforcedFlowPreconditioner ( const LagrangeEnforcedFlowPreconditioner & )

delete

Broken copy constructor.

Member Function Documentation

◆ clean_up_memory()

void oomph::LagrangeEnforcedFlowPreconditioner::clean_up_memory ( )

virtual

Clears the memory.

Reimplemented from oomph::Preconditioner.

   {
     // clean the block preconditioner base class memory
     this->clear_block_preconditioner_base();
  
     // Delete the Navier-Stokes preconditioner pointer if we have created it.
     if (Using_superlu_ns_preconditioner && Navier_stokes_preconditioner_pt != 0)
     {
       delete Navier_stokes_preconditioner_pt;
       Navier_stokes_preconditioner_pt = 0;
     }
  
     Preconditioner_has_been_setup = false;
   } // func LagrangeEnforcedFlowPreconditioner::clean_up_memory

References oomph::BlockPreconditioner< CRDoubleMatrix >::clear_block_preconditioner_base(), Navier_stokes_preconditioner_pt, Preconditioner_has_been_setup, and Using_superlu_ns_preconditioner.

Referenced by setup(), and ~LagrangeEnforcedFlowPreconditioner().

◆ operator=()

void oomph::LagrangeEnforcedFlowPreconditioner::operator= ( const LagrangeEnforcedFlowPreconditioner & )

delete

Broken assignment operator.

◆ preconditioner_solve()

void oomph::LagrangeEnforcedFlowPreconditioner::preconditioner_solve	(	const DoubleVector &	r,
		DoubleVector &	z
	)

virtual

Apply the preconditioner. r is the residual (rhs), z will contain the solution.

///////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////// Apply the preconditioner. r is the residual (rhs), z will contain the solution.

Implements oomph::Preconditioner.

   {
 #ifdef PARANOID
     if (Preconditioner_has_been_setup == false)
     {
       std::ostringstream error_message;
       error_message << "setup() must be called before using "
                     << "preconditioner_solve()";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
     if (z.built())
     {
       if (z.nrow() != r.nrow())
       {
         std::ostringstream error_message;
         error_message << "The vectors z and r must have the same number of "
                       << "of global rows";
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // if z is not setup then give it the same distribution
     if (!z.distribution_pt()->built())
     {
       z.build(r.distribution_pt(), 0.0);
     }
  
     // Working vectors.
     DoubleVector temp_vec;
     DoubleVector another_temp_vec;
     DoubleVector yet_another_temp_vec;
  
  
     // -----------------------------------------------------------------------
     // Step 1 - apply approximate W block inverse to Lagrange multiplier
     // unknowns
     // -----------------------------------------------------------------------
     // For each subsystem associated with each Lagrange multiplier, we loop
     // through and:
     // 1) extract the block entries from r
     // 2) apply the inverse
     // 3) return the entries to z.
     for (unsigned l_i = 0; l_i < N_lagrange_doftypes; l_i++)
     {
       // The Lagrange multiplier block type.
       const unsigned l_ii = N_fluid_doftypes + l_i;
  
       // Extract the block
       this->get_block_vector(l_ii, r, temp_vec);
  
       // Apply the inverse.
       const unsigned vec_nrow_local = temp_vec.nrow_local();
       double* vec_values_pt = temp_vec.values_pt();
       for (unsigned i = 0; i < vec_nrow_local; i++)
       {
         vec_values_pt[i] = vec_values_pt[i] * Inv_w_diag_values[l_i][i];
       } // for
  
       // Return the unknowns
       this->return_block_vector(l_ii, temp_vec, z);
  
       // Clear vectors.
       temp_vec.clear();
     } // for
  
     // -----------------------------------------------------------------------
     // Step 2 - apply the augmented Navier-Stokes matrix inverse to the
     // velocity and pressure unknowns
     // -----------------------------------------------------------------------
  
     // At this point, all vectors are cleared.
     if (Using_superlu_ns_preconditioner)
     {
       // Which block types corresponds to the fluid block types.
       Vector<unsigned> fluid_block_indices(N_fluid_doftypes, 0);
       for (unsigned b = 0; b < N_fluid_doftypes; b++)
       {
         fluid_block_indices[b] = b;
       }
  
       this->get_concatenated_block_vector(fluid_block_indices, r, temp_vec);
  
       // temp_vec contains the (concatenated) fluid rhs.
       Navier_stokes_preconditioner_pt->preconditioner_solve(temp_vec,
                                                             another_temp_vec);
  
       temp_vec.clear();
  
       // Return it to the unknowns.
       this->return_concatenated_block_vector(
         fluid_block_indices, another_temp_vec, z);
  
       another_temp_vec.clear();
     }
     else
     {
       // The Navier-Stokes preconditioner is a block preconditioner.
       // Thus is handles all of the block vector extraction and returns.
       Navier_stokes_preconditioner_pt->preconditioner_solve(r, z);
     }
   } // end of preconditioner_solve

References b, oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::built(), oomph::DoubleVector::clear(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vector(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_concatenated_block_vector(), i, Inv_w_diag_values, N_fluid_doftypes, N_lagrange_doftypes, Navier_stokes_preconditioner_pt, oomph::DistributableLinearAlgebraObject::nrow(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Preconditioner_has_been_setup, oomph::Preconditioner::preconditioner_solve(), UniformPSDSelfTest::r, oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vector(), oomph::BlockPreconditioner< CRDoubleMatrix >::return_concatenated_block_vector(), Using_superlu_ns_preconditioner, and oomph::DoubleVector::values_pt().

◆ scaling_sigma()

double oomph::LagrangeEnforcedFlowPreconditioner::scaling_sigma ( ) const

inline

Read (const) function to get the scaling sigma.

     {
       return Scaling_sigma;
     }

References Scaling_sigma.

Referenced by set_scaling_sigma().

◆ set_meshes()

void oomph::LagrangeEnforcedFlowPreconditioner::set_meshes ( const Vector< Mesh * > & mesh_pt )

Set the meshes, the first mesh in the vector must be the bulk mesh.

   {
     // There should be at least two meshes passed to this preconditioner.
     const unsigned nmesh = mesh_pt.size();
  
 #ifdef PARANOID
     if (nmesh < 2)
     {
       std::ostringstream err_msg;
       err_msg << "There should be at least two meshes.\n";
       throw OomphLibError(
         err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // Check that all pointers are not null
     for (unsigned mesh_i = 0; mesh_i < nmesh; mesh_i++)
     {
       if (mesh_pt[mesh_i] == 0)
       {
         std::ostringstream err_msg;
         err_msg << "The pointer mesh_pt[" << mesh_i << "] is null.\n";
         throw OomphLibError(
           err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
     }
  
     // We assume that the first mesh is the Navier-Stokes "bulk" mesh.
     // To check this, the elemental dimension must be the same as the
     // nodal (spatial) dimension.
     //
     // We store the elemental dimension i.e. the number of local coordinates
     // required to parametrise its geometry.
     const unsigned elemental_dim = mesh_pt[0]->elemental_dimension();
  
     // The dimension of the nodes in the first element in the (supposedly)
     // bulk mesh.
     const unsigned nodal_dim = mesh_pt[0]->nodal_dimension();
  
     // Check if the first mesh is the "bulk" mesh.
     // Here we assume only one mesh contains "bulk" elements.
     // All subsequent meshes contain block preconditionable elements which
     // re-classify the bulk velocity DOFs to constrained velocity DOFs.
     if (elemental_dim != nodal_dim)
     {
       std::ostringstream err_msg;
       err_msg << "In the first mesh, the elements have elemental dimension "
               << "of " << elemental_dim << ",\n"
               << "with a nodal dimension of " << nodal_dim << ".\n"
               << "The first mesh does not contain 'bulk' elements.\n"
               << "Please re-order your mesh_pt vector.\n";
  
       throw OomphLibError(
         err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Set the number of meshes
     this->set_nmesh(nmesh);
  
     // Set the meshes
     for (unsigned mesh_i = 0; mesh_i < nmesh; mesh_i++)
     {
       this->set_mesh(mesh_i, mesh_pt[mesh_i]);
     }
  
     // We also store the meshes and number of meshes locally in this class.
     // This may seem slightly redundant, since we always have all the meshes
     // stored in the upper most master block preconditioner.
     // But at the moment there is no mapping/look up scheme between
     // master and subsidiary block preconditioners for meshes.
     //
     // If this is a subsidiary block preconditioner, we don't know which of
     // the master's meshes belong to us. We need this information to set up
     // look up lists in the function setup(...).
     // Thus we store them local to this class.
     My_mesh_pt = mesh_pt;
     My_nmesh = nmesh;
   } // function set_meshes

References oomph::Mesh::elemental_dimension(), oomph::BlockPreconditioner< CRDoubleMatrix >::mesh_pt(), My_mesh_pt, My_nmesh, oomph::BlockPreconditioner< CRDoubleMatrix >::nmesh(), oomph::Mesh::nodal_dimension(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::BlockPreconditioner< CRDoubleMatrix >::set_mesh(), and oomph::BlockPreconditioner< CRDoubleMatrix >::set_nmesh().

Referenced by TiltedCavityProblem< ELEMENT >::TiltedCavityProblem().

◆ set_navier_stokes_preconditioner()

void oomph::LagrangeEnforcedFlowPreconditioner::set_navier_stokes_preconditioner ( Preconditioner * new_ns_preconditioner_pt = 0 )

Set a new Navier-Stokes matrix preconditioner (inexact solver)

Function to set a new momentum matrix preconditioner (inexact solver)

   {
     // Check if pointer is non-zero.
     if (new_ns_preconditioner_pt == 0)
     {
       std::ostringstream warning_stream;
       warning_stream << "WARNING: \n"
                      << "The LSC preconditioner point is null.\n"
                      << "Using default (SuperLU) preconditioner.\n"
                      << std::endl;
       OomphLibWarning(
         warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
  
       Navier_stokes_preconditioner_pt = 0;
       Using_superlu_ns_preconditioner = true;
     }
     else
     {
       // If the default SuperLU preconditioner has been used
       // clean it up now...
       if (Using_superlu_ns_preconditioner &&
           Navier_stokes_preconditioner_pt != 0)
       {
         delete Navier_stokes_preconditioner_pt;
         Navier_stokes_preconditioner_pt = 0;
       }
  
       Navier_stokes_preconditioner_pt = new_ns_preconditioner_pt;
       Using_superlu_ns_preconditioner = false;
     }
   } // func set_navier_stokes_preconditioner

References Navier_stokes_preconditioner_pt, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Using_superlu_ns_preconditioner.

Referenced by TiltedCavityProblem< ELEMENT >::TiltedCavityProblem().

◆ set_scaling_sigma()

void oomph::LagrangeEnforcedFlowPreconditioner::set_scaling_sigma ( const double & scaling_sigma )

inline

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.

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

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, scaling_sigma(), Scaling_sigma, and Use_norm_f_for_scaling_sigma.

◆ set_superlu_for_navier_stokes_preconditioner()

void oomph::LagrangeEnforcedFlowPreconditioner::set_superlu_for_navier_stokes_preconditioner ( )

inline

Set Navier-Stokes matrix preconditioner (inexact solver) to SuperLU

     {
       if (!Using_superlu_ns_preconditioner)
       {
         delete Navier_stokes_preconditioner_pt;
         Navier_stokes_preconditioner_pt = new SuperLUPreconditioner;
         Using_superlu_ns_preconditioner = true;
       }
     }

References Navier_stokes_preconditioner_pt, and Using_superlu_ns_preconditioner.

Referenced by TiltedCavityProblem< ELEMENT >::TiltedCavityProblem().

◆ setup()

void oomph::LagrangeEnforcedFlowPreconditioner::setup ( )

virtual

Setup method for the LagrangeEnforcedFlowPreconditioner.

Setup the Lagrange enforced flow preconditioner. This extracts blocks corresponding to the velocity and Lagrange multiplier unknowns, creates the matrices actually needed in the application of the preconditioner and deletes what can be deleted... Note that this preconditioner needs a CRDoubleMatrix.

Implements oomph::Preconditioner.

   {
     // clean
     this->clean_up_memory();
  
 #ifdef PARANOID
     // Paranoid check that meshes have been set. In this preconditioner, we
     // always have to set meshes.
     if (My_nmesh == 0)
     {
       std::ostringstream err_msg;
       err_msg << "There are no meshes set. Please call set_meshes(...)\n"
               << "with at least two mesh.";
       throw OomphLibError(
         err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // -----------------------------------------------------------------------
     // Step 1 - Construct the dof_to_block_map vector.
     // -----------------------------------------------------------------------
     // Assumption: The first mesh is always the "bulk" mesh
     // (Navier-Stokes mesh), which contains block preconditionable
     // Navier-Stokes elements. Thus the bulk elements classify the velocity
     // and pressure degrees of freedom (ndof_types = 3(4) in 2(3)D).
     // All subsequent meshes contain the constrained velocity DOF types,
     // then the Lagrange multiplier DOF types.
     //
     // Thus, a general ordering of DOF types (in 3D) follows the ordering:
     //
     //  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 ...] ...
     //
     // where the square brackets [] represent the DOF types in each mesh.
     //
     // Example:
     // Consider the case of imposing parallel outflow (3 constrained velocity
     // DOF types and 2 Lagrange multiplier DOF types) and tangential flow (3
     // constrained velocity DOF types and 1 Lagrange multiplier DOF type)
     // along two different boundaries in 3D. The resulting natural ordering of
     // the DOF types is:
     //
     // [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]
     //
     //
     // In our implementation, the desired block structure is:
     // | u v w | up vp wp | ut vt wt | p | Lp1 Lp2 Lt1 |
     //
     // The dof_to_block_map should have the following construction:
     //
     //    dof_to_block_map[$dof_number] = $block_number
     //
     // Thus, the dof_to_block_map for the example above should be:
     //
     // To achieve this, we use the dof_to_block_map:
     // dof_to_block_map = [0 1 2 9 3  4  5  10  11  6  7  8  12]
     //
     // To generalise the construction of the dof_to_block_map vector
     // (to work for any number of meshes), we first require some auxiliary
     // variables to aid us in this endeavour.
     // -----------------------------------------------------------------------
  
     // Set up the My_ndof_types_in_mesh vector.
     // If this is already constructed, we reuse it instead.
     if (My_ndof_types_in_mesh.size() == 0)
     {
       for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
       {
         My_ndof_types_in_mesh.push_back(My_mesh_pt[mesh_i]->ndof_types());
       }
     }
  
     // Get the spatial dimension of the problem.
     unsigned spatial_dim = My_mesh_pt[0]->nodal_dimension();
  
     // Get the number of DOF types.
     unsigned n_dof_types = ndof_types();
  
 #ifdef PARANOID
     // We cannot check which DOF types are "correct" in the sense that there
     // is a distinction between bulk and constrained velocity DOF tyoes.
     // But we can at least check if the ndof_types matches the total number
     // of DOF types from the meshes.
     //
     // This check is not necessary for the upper most master block
     // preconditioner since the ndof_types() is calculated by looping
     // through the meshes!
     //
     // This check is useful if this is a subsidiary block preconditioner and
     // incorrect DOF type merging has taken place.
     if (is_subsidiary_block_preconditioner())
     {
       unsigned tmp_ndof_types = 0;
       for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
       {
         tmp_ndof_types += My_ndof_types_in_mesh[mesh_i];
       }
  
       if (tmp_ndof_types != n_dof_types)
       {
         std::ostringstream err_msg;
         err_msg << "The number of DOF types are incorrect.\n"
                 << "The total DOF types from the meshes is: " << tmp_ndof_types
                 << ".\n"
                 << "The number of DOF types from "
                 << "BlockPreconditioner::ndof_types() is " << n_dof_types
                 << ".\n";
         throw OomphLibError(
           err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // The number of velocity DOF types: We assume that all meshes classify
     // at least some of the velocity DOF types (bulk/constrained), thus the
     // total number of velocity DOF types is the spatial dimension multiplied
     // by the number of meshes.
     N_velocity_doftypes = My_nmesh * spatial_dim;
  
     // Fluid has +1 for the pressure.
     N_fluid_doftypes = N_velocity_doftypes + 1;
  
     // The rest are Lagrange multiplier DOF types.
     N_lagrange_doftypes = n_dof_types - N_fluid_doftypes;
  
     //   Now construct the DOF to block map for block_setup()   //
  
     // Now that we have
     //
     // (1) spatial dimension of the problem and
     // (2) the number of DOF types in each of the meshes.
     //
     // We observe that the problem dimension is 3 and
     // My_ndof_type_in_mesh = [4, 5, 4].
     //
     // With these information we can construct the desired block structure:
     // | u v w | up vp wp | ut vt wt | p | Lp1 Lp2 Lt1 |
     //
     // The block structure is determined by the vector dof_to_block_map we
     // give to the function block_setup(...).
  
     // This preconditioner permutes the DOF numbers to get the block numbers.
     // I.e. nblock_types = ndof_types, but they are re-ordered
     // so that we have (in order):
     // 1) bulk velocity DOF types
     // 2) constrained velocity DOF types
     // 3) pressure DOF type
     // 4) Lagrange multiplier DOF types
     //
     // Recall that the natural ordering of the DOF types are ordered first by
     // their meshes, then the elemental DOF type as described by the
     // function get_dof_numbers_for_unknowns().
     //
     // Also recall that (for this problem), we assume/require that every mesh
     // (re-)classify at least some of the velocity DOF types, furthermore,
     // the velocity DOF type classification comes before the
     // pressure / Lagrange multiplier DOF types.
     //
     // Consider the same example as shown previously, repeated here for your
     // convenience:
     //
     // Natural DOF type ordering:
     //  0 1 2 3   4  5  6   7   8    9  10 11 12   <- DOF number.
     // [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1] <- DOF type.
     //
     // Desired block structure:
     //  0 1 2   3  4  5    6  7  8     9   10  11  12   <- block number
     // [u v w | up vp wp | ut vt wt ] [p | Lp1 Lp2 Lt1] <- DOF type
     //
     // dof_to_block_map[$dof_number] = $block_number
     //
     // Required dof_to_block_map:
     //  0 1 2 9 3  4  5  10  11  6  7  8  12   <- block number
     // [u v w p up vp wp Lp1 Lp2 ut vt wt Lt1] <- DOF type
     //
     // Consider the 4th entry of the dof_to_block_map, this represents the
     // pressure DOF type, from the desired block structure we see this takes
     // the block number 9.
     //
     // One way to generalise the construction of the dof_to_block_map  is to
     // lump the first $spatial_dimension number of DOF types
     // from each mesh together, then lump the remaining DOF types together.
     //
     // Notice that the values in the velocity DOF types of the
     // dof_to_block_map vector increases sequentially, from 0 to 8. The values
     // of the Lagrange multiplier DOF types follow the same pattern, from 9 to
     // 12. We follow this construction.
  
     // Storage for the dof_to_block_map vector.
     Vector<unsigned> dof_to_block_map(n_dof_types, 0);
  
     // Index for the dof_to_block_map vector.
     unsigned temp_index = 0;
  
     // Value for the velocity DOF type.
     unsigned velocity_number = 0;
  
     // Value for the pressure/Lagrange multiplier DOF type.
     unsigned pres_lgr_number = N_velocity_doftypes;
  
     // Loop through the number of meshes.
     for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
     {
       // Fill in the velocity DOF types.
       for (unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
       {
         dof_to_block_map[temp_index++] = velocity_number++;
       } // for
  
       // Loop through the pressure/Lagrange multiplier DOF types.
       unsigned ndof_type_in_mesh_i = My_ndof_types_in_mesh[mesh_i];
       for (unsigned doftype_i = spatial_dim; doftype_i < ndof_type_in_mesh_i;
            doftype_i++)
       {
         dof_to_block_map[temp_index++] = pres_lgr_number++;
       } // for
     } // for
  
     // Call the block setup
     this->block_setup(dof_to_block_map);
  
  
     // -----------------------------------------------------------------------
     // Step 2 - Get the infinite norm of Navier-Stokes F block.
     // -----------------------------------------------------------------------
  
     // Extract the velocity blocks.
     DenseMatrix<CRDoubleMatrix*> v_aug_pt(
       N_velocity_doftypes, N_velocity_doftypes, 0);
  
     for (unsigned row_i = 0; row_i < N_velocity_doftypes; row_i++)
     {
       for (unsigned col_i = 0; col_i < N_velocity_doftypes; col_i++)
       {
         v_aug_pt(row_i, col_i) = new CRDoubleMatrix;
         this->get_block(row_i, col_i, *v_aug_pt(row_i, col_i));
       } // for
     } // for
  
     // Now get the infinite norm.
     if (Use_norm_f_for_scaling_sigma)
     {
       Scaling_sigma = -CRDoubleMatrixHelpers::inf_norm(v_aug_pt);
     } // if
  
 #ifdef PARANOID
     // Warning for division by zero.
     if (Scaling_sigma == 0.0)
     {
       std::ostringstream warning_stream;
       warning_stream << "WARNING: " << std::endl
                      << "The scaling (Scaling_sigma) is " << Scaling_sigma
                      << ".\n"
                      << "Division by 0 will occur."
                      << "\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 << std::endl
                      << "Performance may be degraded." << std::endl;
       OomphLibWarning(
         warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // -----------------------------------------------------------------------
     // Step 3 - Augment the constrained fluid blocks.
     // -----------------------------------------------------------------------
     // Loop through the Lagrange multipliers and do three things:
     // For each Lagrange block:
     //   3.1) Extract the mass matrices and store the location of non-zero mass
     //        matrices.
     //
     //   3.2) a) Create inv_w (for the augmentation)
     //        b) Store the diagonal values of inv_w (for preconditioner solve)
     //
     //   3.3) Perform the augmentation: v_aug + m_i * inv(w_i) * m_j
     //
     //   3.4) Clean up memory.
  
     // Storage for the inv w diag values.
     Inv_w_diag_values.clear();
     for (unsigned l_i = 0; l_i < N_lagrange_doftypes; l_i++)
     {
       // Step 3.1: Location and extraction of non-zero mass matrices.
  
       // Storage for the current Lagrange block mass matrices.
       Vector<CRDoubleMatrix*> mm_pt;
       Vector<CRDoubleMatrix*> mmt_pt;
  
       // Block type for the Lagrange multiplier.
       const unsigned lgr_block_num = N_fluid_doftypes + l_i;
  
       // Store the mass matrix locations for the current Lagrange block.
       Vector<unsigned> mm_locations;
  
       // Store the number of mass matrices.
       unsigned n_mm = 0;
  
       // Go along the block columns for the current Lagrange block ROW.
       for (unsigned col_i = 0; col_i < N_velocity_doftypes; col_i++)
       {
         // Get the block matrix for this block column.
         CRDoubleMatrix* mm_temp_pt = new CRDoubleMatrix;
         this->get_block(lgr_block_num, col_i, *mm_temp_pt);
  
         // Check if this is non-zero
         if (mm_temp_pt->nnz() > 0)
         {
           mm_locations.push_back(col_i);
           mm_pt.push_back(mm_temp_pt);
           n_mm++;
         }
         else
         {
           // This is just an empty matrix. No need to keep this.
           delete mm_temp_pt;
         }
       } // loop through the columns of the Lagrange row.
  
 #ifdef PARANOID
       if (n_mm == 0)
       {
         std::ostringstream warning_stream;
         warning_stream << "WARNING:\n"
                        << "There are no mass matrices on Lagrange block row "
                        << l_i << ".\n"
                        << "Perhaps the problem setup is incorrect."
                        << std::endl;
         OomphLibWarning(warning_stream.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Get the transpose of the mass matrices.
       for (unsigned mm_i = 0; mm_i < n_mm; mm_i++)
       {
         // Get the block matrix for this block column.
         CRDoubleMatrix* mm_temp_pt = new CRDoubleMatrix;
         this->get_block(mm_locations[mm_i], lgr_block_num, *mm_temp_pt);
  
         if (mm_temp_pt->nnz() > 0)
         {
           mmt_pt.push_back(mm_temp_pt);
         }
         else
         {
           // There should be a non-zero mass matrix here, since L=(L^T)^T
 #ifdef PARANOID
           {
             std::ostringstream warning_stream;
             warning_stream << "WARNING:\n"
                            << "The mass matrix block " << mm_locations[mm_i]
                            << " in L^T block " << l_i << " is zero.\n"
                            << "Perhaps the problem setup is incorrect."
                            << std::endl;
             OomphLibWarning(warning_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
           }
 #endif
         }
       } // loop through the ROW of the Lagrange COLUMN.
  
       //  Step 3.2: Create inv_w and store its diagonal entries.
  
       // Storage for inv_w
       CRDoubleMatrix* inv_w_pt =
         new CRDoubleMatrix(this->Block_distribution_pt[lgr_block_num]);
  
       // Get the number of local rows for this Lagrange block.
       unsigned long l_i_nrow_local =
         this->Block_distribution_pt[lgr_block_num]->nrow_local();
  
       // The first row, for the column offset (required in parallel).
       unsigned l_i_first_row =
         this->Block_distribution_pt[lgr_block_num]->first_row();
  
       // A vector to contain the results of mass matrices squared.
       Vector<double> w_i_diag_values(l_i_nrow_local, 0);
  
       // Create mm*mm^T, and component-wise add the mass matrices
       for (unsigned m_i = 0; m_i < n_mm; m_i++)
       {
         // Create mm*mm^T
         CRDoubleMatrix temp_mm_sqrd;
         temp_mm_sqrd.build(mm_pt[m_i]->distribution_pt());
         mm_pt[m_i]->multiply(*mmt_pt[m_i], temp_mm_sqrd);
  
         // Extract diagonal entries
         Vector<double> m_diag = temp_mm_sqrd.diagonal_entries();
  
         // Loop through the entries, add them.
         for (unsigned long row_i = 0; row_i < l_i_nrow_local; row_i++)
         {
           w_i_diag_values[row_i] += m_diag[row_i];
         }
       }
  
       // Storage for inv_w matrix vectors
       Vector<double> invw_i_diag_values(l_i_nrow_local, 0);
       Vector<int> w_i_column_indices(l_i_nrow_local);
       Vector<int> w_i_row_start(l_i_nrow_local + 1);
  
       // Divide by Scaling_sigma and create the inverse of w.
       for (unsigned long row_i = 0; row_i < l_i_nrow_local; row_i++)
       {
         invw_i_diag_values[row_i] = Scaling_sigma / w_i_diag_values[row_i];
  
         w_i_column_indices[row_i] = row_i + l_i_first_row;
         w_i_row_start[row_i] = row_i;
       }
  
       w_i_row_start[l_i_nrow_local] = l_i_nrow_local;
  
       // Theses are square matrices. So we can use the l_i_nrow_global for the
       // number of columns.
       unsigned long l_i_nrow_global =
         this->Block_distribution_pt[lgr_block_num]->nrow();
       inv_w_pt->build(
         l_i_nrow_global, invw_i_diag_values, w_i_column_indices, w_i_row_start);
  
       Inv_w_diag_values.push_back(invw_i_diag_values);
  
  
       // Step 3.3: Perform the augmentation: v_aug + m_i * inv(w_i) * m_j
  
       // Now we create the augmented matrix in v_aug_pt.
       // v_aug_pt is already re-ordered
       // Loop through the mm_locations
       for (unsigned ii = 0; ii < n_mm; ii++)
       {
         // Location of the mass matrix
         unsigned aug_i = mm_locations[ii];
  
         for (unsigned jj = 0; jj < n_mm; jj++)
         {
           // Location of the mass matrix
           unsigned aug_j = mm_locations[jj];
  
           // Storage for intermediate results.
           CRDoubleMatrix aug_block;
           CRDoubleMatrix another_aug_block;
  
           // aug_block = mmt*inv_w
           mmt_pt[ii]->multiply((*inv_w_pt), (aug_block));
  
           // another_aug_block = aug_block*mm = mmt*inv_w*mm
           aug_block.multiply(*mm_pt[jj], another_aug_block);
  
           // Add the augmentation.
           v_aug_pt(aug_i, aug_j)
             ->add(another_aug_block, *v_aug_pt(aug_i, aug_j));
         } // loop jj
       } // loop ii
  
       // Step 3.4: Clean up memory.
       delete inv_w_pt;
       inv_w_pt = 0;
       for (unsigned m_i = 0; m_i < n_mm; m_i++)
       {
         delete mm_pt[m_i];
         mm_pt[m_i] = 0;
         delete mmt_pt[m_i];
         mmt_pt[m_i] = 0;
       }
     } // loop through Lagrange multipliers.
  
     // -----------------------------------------------------------------------
     // Step 4 - Replace all the velocity blocks
     // -----------------------------------------------------------------------
     // When we replace (DOF) blocks, the indices have to respect the DOF type
     // ordering. This is because only DOF type information is passed between
     // preconditioners. So we need to create the inverse of the
     // dof_to_block_map... a block_to_dof_map!
     //
     // Note: Once the DOF blocks have been replaced, further calls to
     // get_block at this preconditioning hierarchy level and/or lower will use
     // the (nearest) replaced blocks.
     Vector<unsigned> block_to_dof_map(n_dof_types, 0);
     for (unsigned dof_i = 0; dof_i < n_dof_types; dof_i++)
     {
       block_to_dof_map[dof_to_block_map[dof_i]] = dof_i;
     }
  
     // Now do the replacement of all blocks in v_aug_pt
     for (unsigned block_row_i = 0; block_row_i < N_velocity_doftypes;
          block_row_i++)
     {
       unsigned temp_dof_row_i = block_to_dof_map[block_row_i];
       for (unsigned block_col_i = 0; block_col_i < N_velocity_doftypes;
            block_col_i++)
       {
         unsigned temp_dof_col_i = block_to_dof_map[block_col_i];
         this->set_replacement_dof_block(
           temp_dof_row_i, temp_dof_col_i, v_aug_pt(block_row_i, block_col_i));
       }
     }
  
     // -----------------------------------------------------------------------
     // Step 5 - Set up Navier-Stokes preconditioner
     // If the subsidiary fluid preconditioner is a block preconditioner:
     //   5.1) Set up the dof_number_in_master_map
     //   5.2) Set up the doftype_coarsen_map
     // otherwise:
     //   5.3) Concatenate the fluid matrices into one big matrix and pass it
     //        to the solver.
     // -----------------------------------------------------------------------
  
     // First we determine if we're using a block preconditioner or not.
     BlockPreconditioner<CRDoubleMatrix>* navier_stokes_block_preconditioner_pt =
       dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
         Navier_stokes_preconditioner_pt);
     Navier_stokes_preconditioner_is_block_preconditioner = true;
     if (navier_stokes_block_preconditioner_pt == 0)
     {
       Navier_stokes_preconditioner_is_block_preconditioner = false;
     }
  
     // If the Navier-Stokes preconditioner is a block preconditioner, then we
     // need to turn it into a subsidiary block preconditioner.
     if (Navier_stokes_preconditioner_is_block_preconditioner)
     {
       // Assumption: All Navier-Stokes block preconditioners take $dim
       // number of velocity DOF types and 1 pressure DOF type.
  
       // Step 5.1: Set up the dof_number_in_master_map
  
       // First we create the dof_number_in_master_map vector to pass to the
       // subsidiary preconditioner's
       // turn_into_subsidiary_block_preconditioner(...) function.
       //
       // This vector maps the subsidiary DOF numbers and the master DOF
       // numbers and has the construction:
       //
       //   dof_number_in_master_map[subsidiary DOF number] = master DOF number
       //
       // Example: Using the example above, our problem has the natural
       // DOF type ordering:
       //
       //  0 1 2 3   4  5  6   7   8    9  10 11 12   <- DOF number in master
       // [u v w p] [up vp wp Lp1 Lp2] [ut vt wt Lt1] <- DOF type
       //
       // For now, we ignore the fact that this preconditioner's number of
       // DOF types may be more fine grain than what is assumed by the
       // subsidiary block preconditioner. Normally, the order of the values in
       // dof_number_in_master_map matters, since the indices must match the
       // DOF type in the subsidiary block preconditioner (see the assumption
       // above). For example, if the (subsidiary) LSC block preconditioner
       // requires the DOF type ordering:
       //
       // [0 1 2 3]
       //  u v w p
       //
       // Then the dof_number_in_master_map vector must match the u velocity
       // DOF type in the subsidiary preconditioner with the u velocity in the
       // master preconditioner, etc...
       //
       // However, we shall see (later) that it does not matter in this
       // instance because the DOF type coarsening feature overrides the
       // ordering provided here.
       //
       // For now, we only need to ensure that the subsidiary preconditioner
       // knows about 10 master DOF types (9 velocity and 1 pressure), the
       // ordering does not matter.
       //
       // We pass to the subsidiary block preconditioner the following
       // dof_number_in_master_map:
       // [0 1 2 4  5  6  9  10 11 3] <- DOF number in master
       //  u v w up vp wp ut vt wt p  <- corresponding DOF type
  
       // Variable to keep track of the DOF number.
       unsigned temp_dof_number = 0;
  
       // Storage for the dof_number_in_master_map
       Vector<unsigned> dof_number_in_master_map;
       for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
       {
         // Store the velocity dof types.
         for (unsigned dim_i = 0; dim_i < spatial_dim; dim_i++)
         {
           dof_number_in_master_map.push_back(temp_dof_number + dim_i);
         } // for spatial_dim
  
         // Update the DOF index
         temp_dof_number += My_ndof_types_in_mesh[mesh_i];
       } // for My_nmesh
  
       // Push back the pressure DOF type
       dof_number_in_master_map.push_back(spatial_dim);
  
  
       // Step 5.2 DOF type coarsening.
  
       // Since this preconditioner works with more fine grained DOF types than
       // the subsidiary block preconditioner, we need to tell the subsidiary
       // preconditioner which DOF types to coarsen together.
       //
       // The LSC preconditioner expects 2(3) velocity dof types, and 1
       // pressure DOF types. Thus, we give it this list:
       // u [0, 3, 6]
       // v [1, 4, 7]
       // w [2, 5, 8]
       // p [9]
       //
       // See how the ordering of dof_number_in_master_map (constructed above)
       // does not matter as long as we construct the doftype_coarsen_map
       // correctly.
  
       // Storage for which subsidiary DOF types to coarsen
       Vector<Vector<unsigned>> doftype_coarsen_map;
       for (unsigned direction = 0; direction < spatial_dim; direction++)
       {
         Vector<unsigned> dir_doftypes_vec(My_nmesh, 0);
         for (unsigned mesh_i = 0; mesh_i < My_nmesh; mesh_i++)
         {
           dir_doftypes_vec[mesh_i] = spatial_dim * mesh_i + direction;
         }
         doftype_coarsen_map.push_back(dir_doftypes_vec);
       }
  
       Vector<unsigned> ns_p_vec(1, 0);
  
       // This is simply the number of velocity dof types,
       ns_p_vec[0] = My_nmesh * spatial_dim;
  
       doftype_coarsen_map.push_back(ns_p_vec);
  
       // Turn the Navier-Stokes block preconditioner into a subsidiary block
       // preconditioner.
       navier_stokes_block_preconditioner_pt
         ->turn_into_subsidiary_block_preconditioner(
           this, dof_number_in_master_map, doftype_coarsen_map);
  
       // Call the setup function
       navier_stokes_block_preconditioner_pt->setup(matrix_pt());
     }
     else
     {
       // Step 5.3: This is not a block preconditioner, thus we need to
       // concatenate all the fluid matrices and pass them to the solver.
  
       // Select all the fluid blocks (velocity and pressure)
       VectorMatrix<BlockSelector> f_aug_blocks(N_fluid_doftypes,
                                                N_fluid_doftypes);
       for (unsigned block_i = 0; block_i < N_fluid_doftypes; block_i++)
       {
         for (unsigned block_j = 0; block_j < N_fluid_doftypes; block_j++)
         {
           f_aug_blocks[block_i][block_j].select_block(block_i, block_j, true);
         }
       }
  
       // Note: This will use the replaced blocks.
       CRDoubleMatrix f_aug_block = this->get_concatenated_block(f_aug_blocks);
  
       if (Using_superlu_ns_preconditioner)
       {
         if (Navier_stokes_preconditioner_pt == 0)
         {
           Navier_stokes_preconditioner_pt = new SuperLUPreconditioner;
         }
       }
       else
       {
         if (Navier_stokes_preconditioner_pt == 0)
         {
           std::ostringstream err_msg;
           err_msg << "Not using SuperLUPreconditioner for NS block,\n"
                   << "but the Navier_stokes_preconditioner_pt is null.\n";
           throw OomphLibError(
             err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
         }
       }
  
       // Setup the solver.
       Navier_stokes_preconditioner_pt->setup(&f_aug_block);
  
     } // else
  
     // Clean up memory
     const unsigned v_aug_nrow = v_aug_pt.nrow();
     const unsigned v_aug_ncol = v_aug_pt.ncol();
     for (unsigned v_row = 0; v_row < v_aug_nrow; v_row++)
     {
       for (unsigned v_col = 0; v_col < v_aug_ncol; v_col++)
       {
         delete v_aug_pt(v_row, v_col);
         v_aug_pt(v_row, v_col) = 0;
       }
     }
  
     Preconditioner_has_been_setup = true;
   } // func setup

◆ use_norm_f_for_scaling_sigma()

void oomph::LagrangeEnforcedFlowPreconditioner::use_norm_f_for_scaling_sigma ( )

inline

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.

     {
       Use_norm_f_for_scaling_sigma = true;
     }

References Use_norm_f_for_scaling_sigma.

Member Data Documentation

◆ Inv_w_diag_values

Vector<Vector<double> > oomph::LagrangeEnforcedFlowPreconditioner::Inv_w_diag_values

private

Inverse W values.

Referenced by preconditioner_solve(), and setup().

◆ My_mesh_pt

Vector<Mesh*> oomph::LagrangeEnforcedFlowPreconditioner::My_mesh_pt

private

Storage for the meshes. In our implementation, the first mesh must always be the Navier-Stokes (bulk) mesh, followed by surface meshes.

Referenced by LagrangeEnforcedFlowPreconditioner(), set_meshes(), and setup().

◆ My_ndof_types_in_mesh

Vector<unsigned> oomph::LagrangeEnforcedFlowPreconditioner::My_ndof_types_in_mesh

private

The number of DOF types in each mesh. This is used create various lookup lists.

Referenced by LagrangeEnforcedFlowPreconditioner(), and setup().

◆ My_nmesh

unsigned oomph::LagrangeEnforcedFlowPreconditioner::My_nmesh

private

The number of meshes. This is used to create various lookup lists.

Referenced by LagrangeEnforcedFlowPreconditioner(), set_meshes(), and setup().

◆ N_fluid_doftypes

unsigned oomph::LagrangeEnforcedFlowPreconditioner::N_fluid_doftypes

private

The number of fluid DOF types (including pressure).

Referenced by LagrangeEnforcedFlowPreconditioner(), preconditioner_solve(), and setup().

◆ N_lagrange_doftypes

unsigned oomph::LagrangeEnforcedFlowPreconditioner::N_lagrange_doftypes

private

The number of Lagrange multiplier DOF types.

Referenced by LagrangeEnforcedFlowPreconditioner(), preconditioner_solve(), and setup().

◆ N_velocity_doftypes

unsigned oomph::LagrangeEnforcedFlowPreconditioner::N_velocity_doftypes

private

The number of velocity DOF types.

Referenced by LagrangeEnforcedFlowPreconditioner(), and setup().

◆ Navier_stokes_preconditioner_is_block_preconditioner

bool oomph::LagrangeEnforcedFlowPreconditioner::Navier_stokes_preconditioner_is_block_preconditioner

private

Flag to indicate if the preconditioner for the Navier-Stokes block is a block preconditioner or not.

Referenced by LagrangeEnforcedFlowPreconditioner(), and setup().

◆ Navier_stokes_preconditioner_pt

Preconditioner* oomph::LagrangeEnforcedFlowPreconditioner::Navier_stokes_preconditioner_pt

private

Pointer to the 'preconditioner' for the Navier-Stokes block.

Referenced by clean_up_memory(), LagrangeEnforcedFlowPreconditioner(), preconditioner_solve(), set_navier_stokes_preconditioner(), set_superlu_for_navier_stokes_preconditioner(), and setup().

◆ Preconditioner_has_been_setup

bool oomph::LagrangeEnforcedFlowPreconditioner::Preconditioner_has_been_setup

private

Control flag is true if the preconditioner has been setup (used so we can wipe the data when the preconditioner is called again)

Referenced by clean_up_memory(), LagrangeEnforcedFlowPreconditioner(), preconditioner_solve(), and setup().

◆ Scaling_sigma

double oomph::LagrangeEnforcedFlowPreconditioner::Scaling_sigma

private

Scaling for the augmentation: Scaling_sigma*(LL^T)

Referenced by LagrangeEnforcedFlowPreconditioner(), scaling_sigma(), set_scaling_sigma(), and setup().

◆ Use_norm_f_for_scaling_sigma

bool oomph::LagrangeEnforcedFlowPreconditioner::Use_norm_f_for_scaling_sigma

private

Flag to indicate if we want to use the infinite norm of the Navier-Stokes momentum block for the scaling sigma.

Referenced by LagrangeEnforcedFlowPreconditioner(), set_scaling_sigma(), setup(), and use_norm_f_for_scaling_sigma().

◆ Using_superlu_ns_preconditioner

bool oomph::LagrangeEnforcedFlowPreconditioner::Using_superlu_ns_preconditioner

private

Flag to indicate whether the default NS preconditioner is used.

Referenced by clean_up_memory(), LagrangeEnforcedFlowPreconditioner(), preconditioner_solve(), set_navier_stokes_preconditioner(), set_superlu_for_navier_stokes_preconditioner(), and setup().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ SubsidiaryPreconditionerFctPt

Constructor & Destructor Documentation

◆ LagrangeEnforcedFlowPreconditioner() [1/2]

◆ ~LagrangeEnforcedFlowPreconditioner()

◆ LagrangeEnforcedFlowPreconditioner() [2/2]

Member Function Documentation

◆ clean_up_memory()

◆ operator=()

◆ preconditioner_solve()

◆ scaling_sigma()

◆ set_meshes()

◆ set_navier_stokes_preconditioner()

◆ set_scaling_sigma()

◆ set_superlu_for_navier_stokes_preconditioner()

◆ setup()

◆ use_norm_f_for_scaling_sigma()

Member Data Documentation

◆ Inv_w_diag_values

◆ My_mesh_pt

◆ My_ndof_types_in_mesh

◆ My_nmesh

◆ N_fluid_doftypes

◆ N_lagrange_doftypes

◆ N_velocity_doftypes

◆ Navier_stokes_preconditioner_is_block_preconditioner

◆ Navier_stokes_preconditioner_pt

◆ Preconditioner_has_been_setup

◆ Scaling_sigma

◆ Use_norm_f_for_scaling_sigma

◆ Using_superlu_ns_preconditioner