Public Member Functions | Private Member Functions | Private Attributes | List of all members
oomph::NavierStokesSchurComplementPreconditioner Class Reference

#include <navier_stokes_preconditioners.h>

+ Inheritance diagram for oomph::NavierStokesSchurComplementPreconditioner:

Public Member Functions

 NavierStokesSchurComplementPreconditioner (Problem *problem_pt)
 Constructor - sets defaults for control flags. More...
 
virtual ~NavierStokesSchurComplementPreconditioner ()
 Destructor. More...
 
 NavierStokesSchurComplementPreconditioner (const NavierStokesSchurComplementPreconditioner &)=delete
 Broken copy constructor. More...
 
void set_problem_pt (Problem *problem_pt)
 Broken assignment operator. More...
 
Problem * problem_pt () const
 
void enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements ()
 
void disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements ()
 
void setup ()
 Setup the preconditioner. More...
 
void preconditioner_solve (const DoubleVector &r, DoubleVector &z)
 Apply preconditioner to Vector r. More...
 
void set_navier_stokes_mesh (Mesh *mesh_pt, const bool &allow_multiple_element_type_in_navier_stokes_mesh=false)
 
void set_p_preconditioner (Preconditioner *new_p_preconditioner_pt)
 Function to set a new pressure matrix preconditioner (inexact solver) More...
 
void set_p_superlu_preconditioner ()
 
void set_f_preconditioner (Preconditioner *new_f_preconditioner_pt)
 Function to set a new momentum matrix preconditioner (inexact solver) More...
 
void use_lsc ()
 Use LSC version of the preconditioner. More...
 
void use_fp ()
 Use Fp version of the preconditioner. More...
 
void set_f_superlu_preconditioner ()
 
void enable_doc_time ()
 Enable documentation of time. More...
 
void disable_doc_time ()
 Disable documentation of time. More...
 
void clean_up_memory ()
 Helper function to delete preconditioner data. More...
 
void enable_robin_for_fp ()
 Use Robin BC elements for the Fp preconditioner. More...
 
void disable_robin_for_fp ()
 Don't use Robin BC elements for the Fp preconditioenr. More...
 
void pin_first_pressure_dof_in_press_adv_diff ()
 
void unpin_first_pressure_dof_in_press_adv_diff ()
 
template<class ELEMENT >
void validate (DocInfo &doc_info, Problem *orig_problem_pt)
 
void pin_all_non_pressure_dofs ()
 Pin all non-pressure dofs. More...
 
void get_pressure_advection_diffusion_matrix (CRDoubleMatrix &fp_matrix)
 Get the pressure advection diffusion matrix. More...
 
void reset_pin_status ()
 Reset pin status of all values. More...
 
void setup (DoubleMatrixBase *matrix_pt)
 
void setup (const Problem *problem_pt, DoubleMatrixBase *matrix_pt)
 
virtual void setup ()=0
 
- 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 Member Functions

void assemble_inv_press_and_veloc_mass_matrix_diagonal (CRDoubleMatrix *&inv_p_mass_pt, CRDoubleMatrix *&inv_v_mass_pt, const bool &do_both)
 

Private Attributes

Preconditioner * P_preconditioner_pt
 Pointer to the 'preconditioner' for the pressure matrix. More...
 
Preconditioner * F_preconditioner_pt
 Pointer to the 'preconditioner' for the F matrix. More...
 
bool Using_default_f_preconditioner
 flag indicating whether the default F preconditioner is used More...
 
bool Using_default_p_preconditioner
 flag indicating whether the default P preconditioner is used More...
 
bool Preconditioner_has_been_setup
 
bool Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements
 
bool F_preconditioner_is_block_preconditioner
 
bool Doc_time
 Set Doc_time to true for outputting results of timings. More...
 
MatrixVectorProduct * QBt_mat_vec_pt
 MatrixVectorProduct operator for Qv^{-1} Bt. More...
 
MatrixVectorProduct * Bt_mat_vec_pt
 MatrixVectorProduct operator for Bt. More...
 
MatrixVectorProduct * F_mat_vec_pt
 MatrixVectorProduct operator for F. More...
 
MatrixVectorProduct * E_mat_vec_pt
 MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant) More...
 
Mesh * Navier_stokes_mesh_pt
 
bool Allow_multiple_element_type_in_navier_stokes_mesh
 
bool Use_LSC
 Boolean to indicate use of LSC (true) or Fp (false) variant. More...
 
bool Use_robin_for_fp
 Use Robin BC elements for Fp preconditoner? More...
 
std::map< Data *, std::vector< int > > Eqn_number_backup
 
bool Pin_first_pressure_dof_in_press_adv_diff
 
Problem * Problem_pt
 

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 least-squares commutator (LSC; formerly BFBT) Navier Stokes preconditioner. It uses blocks corresponding to the velocity and pressure unknowns, i.e. there are a total of 2x2 blocks, and all velocity components are treated as a single block of unknowns.

Here are the details: An "ideal" Navier-Stokes preconditioner would solve the system

\[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf D} & {\bf 0} \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right) \]

where \( {\bf F}\) is the momentum block, \( {\bf G} \) the discrete gradient operator, and \( {\bf D}\) the discrete divergence operator. (For unstabilised elements, we have \( {\bf D} = {\bf G}^T \) and in much of the literature the divergence matrix is denoted by \( {\bf B} \) .) The use of this preconditioner would ensure the convergence of any iterative linear solver in a single iteration but its application is, of course, exactly as expensive as a direct solve. The LSC/BFBT preconditioner replaces the exact Jacobian by a block-triangular approximation

\[ \left( \begin{array}{cc} {\bf F} & {\bf G} \\ {\bf 0} & -{\bf M}_s \end{array} \right) \left( \begin{array}{c} {\bf z}_u \\ {\bf z}_p \end{array} \right) = \left( \begin{array}{c} {\bf r}_u \\ {\bf r}_p \end{array} \right), \]

where \({\bf M}_s\) is an approximation to the pressure Schur-complement \( {\bf S} = {\bf D} {\bf F}^{-1}{\bf G}. \) This system can be solved in two steps:

  1. Solve the second row for \( {\bf z}_p\) via

    \[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \]

  2. Given \( {\bf z}_p \) , solve the first row for \( {\bf z}_u\) via

    \[ {\bf z}_u = {\bf F}^{-1} \big( {\bf r}_u - {\bf G} {\bf z}_p \big) \]

In the LSC/BFBT preconditioner, the action of the inverse pressure Schur complement

\[ {\bf z}_p = - {\bf M}_s^{-1} {\bf r}_p \]

is approximated by

\[ {\bf z}_p = - \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} \big({\bf D} \widehat{\bf Q}^{-1}{\bf F} \widehat{\bf Q}^{-1}{\bf G}\big) \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big)^{-1} {\bf r}_p, \]

where \( \widehat{\bf Q} \) is the diagonal of the velocity mass matrix. The evaluation of this expression involves two linear solves involving the matrix

\[ {\bf P} = \big({\bf D} \widehat{\bf Q}^{-1}{\bf G} \big) \]

which has the character of a matrix arising from the discretisation of a Poisson problem on the pressure space. We also have to evaluate matrix-vector products with the matrix

\[ {\bf E}={\bf D}\widehat{\bf Q}^{-1}{\bf F}\widehat{\bf Q}^{-1}{\bf G} \]

Details of the theory can be found in "Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics" by Howard C. Elman, David J. Silvester, and Andrew J. Wathen, published by Oxford University Press, 2006.

In our implementation of the preconditioner, the linear systems can either be solved "exactly", using SuperLU (in its incarnation as an exact preconditioner; this is the default) or by any other Preconditioner (inexact solver) specified via the access functions

NavierStokesSchurComplementPreconditioner::set_f_preconditioner(...)
oomph::NavierStokesSchurComplementPreconditioner::set_f_preconditioner
void set_f_preconditioner(Preconditioner *new_f_preconditioner_pt)
Function to set a new momentum matrix preconditioner (inexact solver)
Definition: navier_stokes_preconditioners.h:769

or

NavierStokesSchurComplementPreconditioner::set_p_preconditioner(...)
oomph::NavierStokesSchurComplementPreconditioner::set_p_preconditioner
void set_p_preconditioner(Preconditioner *new_p_preconditioner_pt)
Function to set a new pressure matrix preconditioner (inexact solver)
Definition: navier_stokes_preconditioners.h:745

Constructor & Destructor Documentation

◆ NavierStokesSchurComplementPreconditioner() [1/2]

oomph::NavierStokesSchurComplementPreconditioner::NavierStokesSchurComplementPreconditioner ( Problem *  problem_pt)
inline

Constructor - sets defaults for control flags.

611  : BlockPreconditioner<CRDoubleMatrix>(), Problem_pt(problem_pt)
612  {
613  // Insist that all elements are of type
614  // NavierStokesElementWithDiagonalMassMatrices and issue a warning
615  // if one is not.
616  Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements = false;
617 
618  // Use Robin BC elements for Fp preconditioner -- yes by default
619  Use_robin_for_fp = true;
620 
621  // Flag to indicate that the preconditioner has been setup
622  // previously -- if setup() is called again, data can
623  // be wiped.
624  Preconditioner_has_been_setup = false;
625 
626  // By default we use the LSC version
627  Use_LSC = true;
628 
629  // By default we use SuperLU for both p and f blocks
630  Using_default_p_preconditioner = true;
631  Using_default_f_preconditioner = true;
632 
633  // Pin pressure dof in press adv diff problem for Fp precond
634  Pin_first_pressure_dof_in_press_adv_diff = true;
635 
636  Navier_stokes_mesh_pt = 0;
637 
638  // Set default preconditioners (inexact solvers) -- they are
639  // members of this class!
640  P_preconditioner_pt = 0;
641  F_preconditioner_pt = 0;
642 
643  // set Doc_time to false
644  Doc_time = false;
645 
646  // null the off diagonal Block matrix pt
647  Bt_mat_vec_pt = 0;
648 
649  // null the F matrix vector product helper
650  F_mat_vec_pt = 0;
651 
652  // null the QBt matrix vector product pt
653  QBt_mat_vec_pt = 0;
654 
655  // null the E matrix vector product helper in Fp
656  E_mat_vec_pt = 0;
657 
658  // Initially assume that there are no multiple element types in the
659  // Navier-Stokes mesh
660  Allow_multiple_element_type_in_navier_stokes_mesh = false;
661  }
oomph::NavierStokesSchurComplementPreconditioner::Doc_time
bool Doc_time
Set Doc_time to true for outputting results of timings.
Definition: navier_stokes_preconditioners.h:1155
oomph::NavierStokesSchurComplementPreconditioner::F_preconditioner_pt
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
Definition: navier_stokes_preconditioners.h:1121
oomph::NavierStokesSchurComplementPreconditioner::Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements
bool Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements
Definition: navier_stokes_preconditioners.h:1137
oomph::NavierStokesSchurComplementPreconditioner::Problem_pt
Problem * Problem_pt
Definition: navier_stokes_preconditioners.h:1196
oomph::NavierStokesSchurComplementPreconditioner::Preconditioner_has_been_setup
bool Preconditioner_has_been_setup
Definition: navier_stokes_preconditioners.h:1132
oomph::NavierStokesSchurComplementPreconditioner::Use_robin_for_fp
bool Use_robin_for_fp
Use Robin BC elements for Fp preconditoner?
Definition: navier_stokes_preconditioners.h:1181
oomph::NavierStokesSchurComplementPreconditioner::Pin_first_pressure_dof_in_press_adv_diff
bool Pin_first_pressure_dof_in_press_adv_diff
Definition: navier_stokes_preconditioners.h:1192
oomph::NavierStokesSchurComplementPreconditioner::Navier_stokes_mesh_pt
Mesh * Navier_stokes_mesh_pt
Definition: navier_stokes_preconditioners.h:1171
oomph::NavierStokesSchurComplementPreconditioner::F_mat_vec_pt
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F.
Definition: navier_stokes_preconditioners.h:1164
oomph::NavierStokesSchurComplementPreconditioner::E_mat_vec_pt
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
Definition: navier_stokes_preconditioners.h:1167
oomph::NavierStokesSchurComplementPreconditioner::Use_LSC
bool Use_LSC
Boolean to indicate use of LSC (true) or Fp (false) variant.
Definition: navier_stokes_preconditioners.h:1178
oomph::NavierStokesSchurComplementPreconditioner::Using_default_p_preconditioner
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
Definition: navier_stokes_preconditioners.h:1127
oomph::NavierStokesSchurComplementPreconditioner::QBt_mat_vec_pt
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for Qv^{-1} Bt.
Definition: navier_stokes_preconditioners.h:1158
oomph::NavierStokesSchurComplementPreconditioner::Using_default_f_preconditioner
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
Definition: navier_stokes_preconditioners.h:1124
oomph::NavierStokesSchurComplementPreconditioner::Allow_multiple_element_type_in_navier_stokes_mesh
bool Allow_multiple_element_type_in_navier_stokes_mesh
Definition: navier_stokes_preconditioners.h:1175
oomph::NavierStokesSchurComplementPreconditioner::P_preconditioner_pt
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
Definition: navier_stokes_preconditioners.h:1118
oomph::NavierStokesSchurComplementPreconditioner::Bt_mat_vec_pt
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt.
Definition: navier_stokes_preconditioners.h:1161
oomph::NavierStokesSchurComplementPreconditioner::problem_pt
Problem * problem_pt() const
Definition: navier_stokes_preconditioners.h:689

References Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements, Allow_multiple_element_type_in_navier_stokes_mesh, Bt_mat_vec_pt, Doc_time, E_mat_vec_pt, F_mat_vec_pt, F_preconditioner_pt, Navier_stokes_mesh_pt, P_preconditioner_pt, Pin_first_pressure_dof_in_press_adv_diff, Preconditioner_has_been_setup, QBt_mat_vec_pt, Use_LSC, Use_robin_for_fp, Using_default_f_preconditioner, and Using_default_p_preconditioner.

◆ ~NavierStokesSchurComplementPreconditioner()

virtual oomph::NavierStokesSchurComplementPreconditioner::~NavierStokesSchurComplementPreconditioner ( )
inlinevirtual

Destructor.

665  {
666  clean_up_memory();
667  }
oomph::NavierStokesSchurComplementPreconditioner::clean_up_memory
void clean_up_memory()
Helper function to delete preconditioner data.
Definition: navier_stokes_preconditioners.cc:1672

References clean_up_memory().

◆ NavierStokesSchurComplementPreconditioner() [2/2]

oomph::NavierStokesSchurComplementPreconditioner::NavierStokesSchurComplementPreconditioner ( const NavierStokesSchurComplementPreconditioner &  )
delete

Broken copy constructor.

Member Function Documentation

◆ assemble_inv_press_and_veloc_mass_matrix_diagonal()

void oomph::NavierStokesSchurComplementPreconditioner::assemble_inv_press_and_veloc_mass_matrix_diagonal ( CRDoubleMatrix *&  inv_p_mass_pt,
CRDoubleMatrix *&  inv_v_mass_pt,
const bool &  do_both 
)
private

Helper function to assemble the inverse diagonals of the pressure and velocity mass matrices from the elemental contributions defined in NavierStokesEquations<DIM>. If do_both=true, both are computed, otherwise only the velocity mass matrix (the LSC version of the preconditioner only needs that one)

Helper function to assemble the inverse diagonals of the pressure and velocity mass matrix from the elemental contributions defined in NavierStokesElementWithDiagonalMassMatrices:: get_pressure_and_velocity_mass_matrix_diagonal(...) If do_both=true, both are computed, otherwise only the velocity mass matrix (the LSC version of the preconditioner only needs that one)

850  {
851  // determine the velocity rows required by this processor
852  unsigned v_first_row = this->block_distribution_pt(0)->first_row();
853  unsigned v_nrow_local = this->block_distribution_pt(0)->nrow_local();
854  unsigned v_nrow = this->block_distribution_pt(0)->nrow();
855 
856  // create storage for the diagonals
857  double* v_values = new double[v_nrow_local];
858  for (unsigned i = 0; i < v_nrow_local; i++)
859  {
860  v_values[i] = 0.0;
861  }
862 
863  // Equivalent information for pressure mass matrix (only needed for
864  // Fp version)
865  unsigned p_first_row = 0;
866  unsigned p_nrow_local = 0;
867  unsigned p_nrow = 0;
868  double* p_values = 0;
869  if (!Use_LSC)
870  {
871  // determine the pressure rows required by this processor
872  p_first_row = this->block_distribution_pt(1)->first_row();
873  p_nrow_local = this->block_distribution_pt(1)->nrow_local();
874  p_nrow = this->block_distribution_pt(1)->nrow();
875 
876  // create storage for the diagonals
877  p_values = new double[p_nrow_local];
878  for (unsigned i = 0; i < p_nrow_local; i++)
879  {
880  p_values[i] = 0.0;
881  }
882  }
883 
884  // if the problem is distributed
885  bool distributed = false;
886 #ifdef OOMPH_HAS_MPI
887  if (problem_pt()->distributed() ||
888  this->master_distribution_pt()->distributed())
889  {
890  distributed = true;
891  }
892 #endif
893 
894  // next we get the diagonal velocity mass matrix data
895  if (distributed)
896  {
897 #ifdef OOMPH_HAS_MPI
898 
899  // the number of processors
900  unsigned nproc = comm_pt()->nproc();
901 
902  // and my rank
903  unsigned my_rank = comm_pt()->my_rank();
904 
905  // determine the rows for which we have lookup rows
906 
907  // if the problem is NOT distributed then we only classify global
908  // equations on this processor to avoid duplication (as every processor
909  // holds every element)
910  unsigned first_lookup_row = 0;
911  unsigned last_lookup_row = 0;
912  first_lookup_row = this->master_distribution_pt()->first_row();
913  last_lookup_row =
914  first_lookup_row + this->master_distribution_pt()->nrow_local() - 1;
915 
916  // find number of local elements
917  unsigned n_el = Navier_stokes_mesh_pt->nelement();
918 
919  // get the master distribution pt
920  const LinearAlgebraDistribution* master_distribution_pt =
921  this->master_distribution_pt();
922 
923  // Do the two blocks (0: veloc; 1: press)
924  unsigned max_block = 0;
925  if (!Use_LSC) max_block = 1;
926  for (unsigned block_index = 0; block_index <= max_block; block_index++)
927  {
928  // Local working variables: Default to velocity
929  unsigned v_or_p_first_row = v_first_row;
930  double* v_or_p_values = v_values;
931  // Switch to pressure
932  if (block_index == 1)
933  {
934  v_or_p_first_row = p_first_row;
935  v_or_p_values = p_values;
936  }
937 
938 
939  // the diagonal mass matrix contributions that have been
940  // classified and should be sent to another processor
941  Vector<double>* classified_contributions_send =
942  new Vector<double>[nproc];
943 
944  // the corresponding block indices
945  Vector<unsigned>* classified_indices_send = new Vector<unsigned>[nproc];
946 
947  // the matrix contributions that cannot be classified by this processor
948  // and therefore must be sent to another for classification
949  Vector<double>* unclassified_contributions_send =
950  new Vector<double>[nproc];
951 
952  // the corresponding global indices that require classification
953  Vector<unsigned>* unclassified_indices_send =
954  new Vector<unsigned>[nproc];
955 
956  // get the velocity or pressure distribution pt
957  const LinearAlgebraDistribution* velocity_or_press_dist_pt =
958  this->block_distribution_pt(block_index);
959 
960  // get the contribution for each element
961  for (unsigned e = 0; e < n_el; e++)
962  {
963  // Get element
964  GeneralisedElement* el_pt = Navier_stokes_mesh_pt->element_pt(e);
965 
966  // check that the element is not halo
967  if (!el_pt->is_halo())
968  {
969  // find number of degrees of freedom in the element
970  // (this is slightly too big because it includes the
971  // pressure dofs but this doesn't matter)
972  unsigned el_dof = el_pt->ndof();
973 
974  // Allocate local storage for the element's contribution to the
975  // mass matrix diagonal
976  Vector<double> el_vmm_diagonal(el_dof, 0.0);
977  Vector<double> el_pmm_diagonal(el_dof, 0.0);
978 
979  unsigned which_one = 2;
980  if (block_index == 1) which_one = 1;
981 
982  NavierStokesElementWithDiagonalMassMatrices* cast_el_pt = 0;
983  cast_el_pt =
984  dynamic_cast<NavierStokesElementWithDiagonalMassMatrices*>(el_pt);
985  if (cast_el_pt != 0)
986  {
987  cast_el_pt->get_pressure_and_velocity_mass_matrix_diagonal(
988  el_pmm_diagonal, el_vmm_diagonal, which_one);
989  }
990 
991  // get the contribution for each dof
992  for (unsigned i = 0; i < el_dof; i++)
993  {
994  // Get the equation number
995  unsigned eqn_number = el_pt->eqn_number(i);
996 
997  // if I have lookup information on this processor
998  if ((eqn_number >= first_lookup_row) &&
999  (eqn_number <= last_lookup_row))
1000  {
1001  // Only use the dofs that we're dealing with here
1002  if (this->block_number(eqn_number) == int(block_index))
1003  {
1004  // get the index in the block
1005  unsigned index = this->index_in_block(eqn_number);
1006 
1007  // determine which processor requires the block index
1008  for (unsigned p = 0; p < nproc; p++)
1009  {
1010  if ((index >= velocity_or_press_dist_pt->first_row(p)) &&
1011  (index < (velocity_or_press_dist_pt->first_row(p) +
1012  velocity_or_press_dist_pt->nrow_local(p))))
1013  {
1014  // if it is required by this processor then add the
1015  // contribution
1016  if (p == my_rank)
1017  {
1018  if (block_index == 0)
1019  {
1020  v_or_p_values[index - v_or_p_first_row] +=
1021  el_vmm_diagonal[i];
1022  }
1023  else if (block_index == 1)
1024  {
1025  v_or_p_values[index - v_or_p_first_row] +=
1026  el_pmm_diagonal[i];
1027  }
1028  }
1029  // otherwise store it for communication
1030  else
1031  {
1032  if (block_index == 0)
1033  {
1034  classified_contributions_send[p].push_back(
1035  el_vmm_diagonal[i]);
1036  classified_indices_send[p].push_back(index);
1037  }
1038  else if (block_index == 1)
1039  {
1040  classified_contributions_send[p].push_back(
1041  el_pmm_diagonal[i]);
1042  classified_indices_send[p].push_back(index);
1043  }
1044  }
1045  }
1046  }
1047  }
1048  }
1049  // if we do not have the lookup information on this processor
1050  // then we send the mass matrix contribution to a processor
1051  // which we know to have the lookup information
1052  // the assumption: the processor for which the master block
1053  // preconditioner distribution will definitely hold the lookup
1054  // data for eqn_number (although others may)
1055  else if (problem_pt()->distributed())
1056  {
1057  // determine which processor requires the block index
1058  unsigned p = 0;
1059  while (
1060  !(eqn_number >= master_distribution_pt->first_row(p) &&
1061  (eqn_number < (master_distribution_pt->first_row(p) +
1062  master_distribution_pt->nrow_local(p)))))
1063  {
1064  p++;
1065  }
1066 
1067  // store the data
1068  if (block_index == 0)
1069  {
1070  unclassified_contributions_send[p].push_back(
1071  el_vmm_diagonal[i]);
1072  unclassified_indices_send[p].push_back(eqn_number);
1073  }
1074  else if (block_index == 1)
1075  {
1076  unclassified_contributions_send[p].push_back(
1077  el_pmm_diagonal[i]);
1078  unclassified_indices_send[p].push_back(eqn_number);
1079  }
1080  }
1081  }
1082  }
1083  }
1084 
1085  // next the unclassified contributions are communicated to
1086  // processors that can classify them
1087 
1088  // first determine how many unclassified rows are to be sent to
1089  // each processor
1090  unsigned* n_unclassified_send = new unsigned[nproc];
1091  for (unsigned p = 0; p < nproc; p++)
1092  {
1093  if (p == my_rank)
1094  {
1095  n_unclassified_send[p] = 0;
1096  }
1097  else
1098  {
1099  n_unclassified_send[p] = unclassified_contributions_send[p].size();
1100  }
1101  }
1102 
1103  // then all-to-all com number of unclassified to be sent / recv
1104  unsigned* n_unclassified_recv = new unsigned[nproc];
1105  MPI_Alltoall(n_unclassified_send,
1106  1,
1107  MPI_UNSIGNED,
1108  n_unclassified_recv,
1109  1,
1110  MPI_UNSIGNED,
1111  comm_pt()->mpi_comm());
1112 
1113  // the base displacement for the sends
1114  MPI_Aint base_displacement;
1115  MPI_Get_address(v_or_p_values, &base_displacement);
1116 
1117  // allocate storage for the data to be received
1118  // and post the sends and recvs
1119  Vector<double*> unclassified_contributions_recv(nproc);
1120  Vector<unsigned*> unclassified_indices_recv(nproc);
1121  Vector<MPI_Request> unclassified_recv_requests;
1122  Vector<MPI_Request> unclassified_send_requests;
1123  Vector<unsigned> unclassified_recv_proc;
1124  for (unsigned p = 0; p < nproc; p++)
1125  {
1126  if (p != my_rank)
1127  {
1128  // recv
1129  if (n_unclassified_recv[p] > 0)
1130  {
1131  unclassified_contributions_recv[p] =
1132  new double[n_unclassified_recv[p]];
1133  unclassified_indices_recv[p] =
1134  new unsigned[n_unclassified_recv[p]];
1135 
1136  // data for the struct data type
1137  MPI_Datatype recv_types[2];
1138  MPI_Aint recv_displacements[2];
1139  int recv_sz[2];
1140 
1141  // contributions
1142  MPI_Type_contiguous(
1143  n_unclassified_recv[p], MPI_DOUBLE, &recv_types[0]);
1144  MPI_Type_commit(&recv_types[0]);
1145  MPI_Get_address(unclassified_contributions_recv[p],
1146  &recv_displacements[0]);
1147  recv_displacements[0] -= base_displacement;
1148  recv_sz[0] = 1;
1149 
1150  // indices
1151  MPI_Type_contiguous(
1152  n_unclassified_recv[p], MPI_UNSIGNED, &recv_types[1]);
1153  MPI_Type_commit(&recv_types[1]);
1154  MPI_Get_address(unclassified_indices_recv[p],
1155  &recv_displacements[1]);
1156  recv_displacements[1] -= base_displacement;
1157  recv_sz[1] = 1;
1158 
1159  // build the final recv type
1160  MPI_Datatype final_recv_type;
1161  MPI_Type_create_struct(
1162  2, recv_sz, recv_displacements, recv_types, &final_recv_type);
1163  MPI_Type_commit(&final_recv_type);
1164 
1165  // and recv
1166  MPI_Request req;
1167  MPI_Irecv(v_or_p_values,
1168  1,
1169  final_recv_type,
1170  p,
1171  0,
1172  comm_pt()->mpi_comm(),
1173  &req);
1174  unclassified_recv_requests.push_back(req);
1175  unclassified_recv_proc.push_back(p);
1176  MPI_Type_free(&recv_types[0]);
1177  MPI_Type_free(&recv_types[1]);
1178  MPI_Type_free(&final_recv_type);
1179  }
1180 
1181  // send
1182  if (n_unclassified_send[p] > 0)
1183  {
1184  // data for the struct data type
1185  MPI_Datatype send_types[2];
1186  MPI_Aint send_displacements[2];
1187  int send_sz[2];
1188 
1189  // contributions
1190  MPI_Type_contiguous(
1191  n_unclassified_send[p], MPI_DOUBLE, &send_types[0]);
1192  MPI_Type_commit(&send_types[0]);
1193  MPI_Get_address(&unclassified_contributions_send[p][0],
1194  &send_displacements[0]);
1195  send_displacements[0] -= base_displacement;
1196  send_sz[0] = 1;
1197 
1198  // indices
1199  MPI_Type_contiguous(
1200  n_unclassified_send[p], MPI_UNSIGNED, &send_types[1]);
1201  MPI_Type_commit(&send_types[1]);
1202  MPI_Get_address(&unclassified_indices_send[p][0],
1203  &send_displacements[1]);
1204  send_displacements[1] -= base_displacement;
1205  send_sz[1] = 1;
1206 
1207  // build the final send type
1208  MPI_Datatype final_send_type;
1209  MPI_Type_create_struct(
1210  2, send_sz, send_displacements, send_types, &final_send_type);
1211  MPI_Type_commit(&final_send_type);
1212 
1213  // and send
1214  MPI_Request req;
1215  MPI_Isend(v_or_p_values,
1216  1,
1217  final_send_type,
1218  p,
1219  0,
1220  comm_pt()->mpi_comm(),
1221  &req);
1222  unclassified_send_requests.push_back(req);
1223  MPI_Type_free(&send_types[0]);
1224  MPI_Type_free(&send_types[1]);
1225  MPI_Type_free(&final_send_type);
1226  }
1227  }
1228  }
1229 
1230  // next classify the data as it is received
1231  unsigned n_unclassified_recv_req = unclassified_recv_requests.size();
1232  while (n_unclassified_recv_req > 0)
1233  {
1234  // get the processor number and remove the completed request
1235  // for the vector of requests
1236  int req_num;
1237  MPI_Waitany(n_unclassified_recv_req,
1238  &unclassified_recv_requests[0],
1239  &req_num,
1240  MPI_STATUS_IGNORE);
1241  unsigned p = unclassified_recv_proc[req_num];
1242  unclassified_recv_requests.erase(unclassified_recv_requests.begin() +
1243  req_num);
1244  unclassified_recv_proc.erase(unclassified_recv_proc.begin() +
1245  req_num);
1246  n_unclassified_recv_req--;
1247 
1248  // next classify the dofs
1249  // and store them for sending to other processors if required
1250  unsigned n_recv = n_unclassified_recv[p];
1251  for (unsigned i = 0; i < n_recv; i++)
1252  {
1253  unsigned eqn_number = unclassified_indices_recv[p][i];
1254  // Only deal with our block unknowns
1255  if (this->block_number(eqn_number) == int(block_index))
1256  {
1257  // get the index in the block
1258  unsigned index = this->index_in_block(eqn_number);
1259 
1260  // determine which processor requires the block index
1261  for (unsigned pp = 0; pp < nproc; pp++)
1262  {
1263  if ((index >= velocity_or_press_dist_pt->first_row(pp)) &&
1264  (index < (velocity_or_press_dist_pt->first_row(pp) +
1265  velocity_or_press_dist_pt->nrow_local(pp))))
1266  {
1267  // if it is required by this processor then add the
1268  // contribution
1269  if (pp == my_rank)
1270  {
1271  v_or_p_values[index - v_or_p_first_row] +=
1272  unclassified_contributions_recv[p][i];
1273  }
1274  // otherwise store it for communication
1275  else
1276  {
1277  double v = unclassified_contributions_recv[p][i];
1278  classified_contributions_send[pp].push_back(v);
1279  classified_indices_send[pp].push_back(index);
1280  }
1281  }
1282  }
1283  }
1284  }
1285 
1286  // clean up
1287  delete[] unclassified_contributions_recv[p];
1288  delete[] unclassified_indices_recv[p];
1289  }
1290  delete[] n_unclassified_recv;
1291 
1292  // now all indices have been classified
1293 
1294  // next the classified contributions are communicated to
1295  // processors that require them
1296 
1297  // first determine how many classified rows are to be sent to
1298  // each processor
1299  unsigned* n_classified_send = new unsigned[nproc];
1300  for (unsigned p = 0; p < nproc; p++)
1301  {
1302  if (p == my_rank)
1303  {
1304  n_classified_send[p] = 0;
1305  }
1306  else
1307  {
1308  n_classified_send[p] = classified_contributions_send[p].size();
1309  }
1310  }
1311 
1312  // then all-to-all number of classified to be sent / recv
1313  unsigned* n_classified_recv = new unsigned[nproc];
1314  MPI_Alltoall(n_classified_send,
1315  1,
1316  MPI_UNSIGNED,
1317  n_classified_recv,
1318  1,
1319  MPI_UNSIGNED,
1320  comm_pt()->mpi_comm());
1321 
1322  // allocate storage for the data to be received
1323  // and post the sends and recvs
1324  Vector<double*> classified_contributions_recv(nproc);
1325  Vector<unsigned*> classified_indices_recv(nproc);
1326  Vector<MPI_Request> classified_recv_requests;
1327  Vector<MPI_Request> classified_send_requests;
1328  Vector<unsigned> classified_recv_proc;
1329  for (unsigned p = 0; p < nproc; p++)
1330  {
1331  if (p != my_rank)
1332  {
1333  // recv
1334  if (n_classified_recv[p] > 0)
1335  {
1336  classified_contributions_recv[p] =
1337  new double[n_classified_recv[p]];
1338  classified_indices_recv[p] = new unsigned[n_classified_recv[p]];
1339 
1340  // data for the struct data type
1341  MPI_Datatype recv_types[2];
1342  MPI_Aint recv_displacements[2];
1343  int recv_sz[2];
1344 
1345  // contributions
1346  MPI_Type_contiguous(
1347  n_classified_recv[p], MPI_DOUBLE, &recv_types[0]);
1348  MPI_Type_commit(&recv_types[0]);
1349  MPI_Get_address(classified_contributions_recv[p],
1350  &recv_displacements[0]);
1351  recv_displacements[0] -= base_displacement;
1352  recv_sz[0] = 1;
1353 
1354  // indices
1355  MPI_Type_contiguous(
1356  n_classified_recv[p], MPI_UNSIGNED, &recv_types[1]);
1357  MPI_Type_commit(&recv_types[1]);
1358  MPI_Get_address(classified_indices_recv[p],
1359  &recv_displacements[1]);
1360  recv_displacements[1] -= base_displacement;
1361  recv_sz[1] = 1;
1362 
1363  // build the final recv type
1364  MPI_Datatype final_recv_type;
1365  MPI_Type_create_struct(
1366  2, recv_sz, recv_displacements, recv_types, &final_recv_type);
1367  MPI_Type_commit(&final_recv_type);
1368 
1369  // and recv
1370  MPI_Request req;
1371  MPI_Irecv(v_or_p_values,
1372  1,
1373  final_recv_type,
1374  p,
1375  0,
1376  comm_pt()->mpi_comm(),
1377  &req);
1378  classified_recv_requests.push_back(req);
1379  classified_recv_proc.push_back(p);
1380  MPI_Type_free(&recv_types[0]);
1381  MPI_Type_free(&recv_types[1]);
1382  MPI_Type_free(&final_recv_type);
1383  }
1384 
1385  // send
1386  if (n_classified_send[p] > 0)
1387  {
1388  // data for the struct data type
1389  MPI_Datatype send_types[2];
1390  MPI_Aint send_displacements[2];
1391  int send_sz[2];
1392 
1393  // contributions
1394  MPI_Type_contiguous(
1395  n_classified_send[p], MPI_DOUBLE, &send_types[0]);
1396  MPI_Type_commit(&send_types[0]);
1397  MPI_Get_address(&classified_contributions_send[p][0],
1398  &send_displacements[0]);
1399  send_displacements[0] -= base_displacement;
1400  send_sz[0] = 1;
1401 
1402  // indices
1403  MPI_Type_contiguous(
1404  n_classified_send[p], MPI_UNSIGNED, &send_types[1]);
1405  MPI_Type_commit(&send_types[1]);
1406  MPI_Get_address(&classified_indices_send[p][0],
1407  &send_displacements[1]);
1408  send_displacements[1] -= base_displacement;
1409  send_sz[1] = 1;
1410 
1411  // build the final send type
1412  MPI_Datatype final_send_type;
1413  MPI_Type_create_struct(
1414  2, send_sz, send_displacements, send_types, &final_send_type);
1415  MPI_Type_commit(&final_send_type);
1416 
1417  // and send
1418  MPI_Request req;
1419  MPI_Isend(v_or_p_values,
1420  1,
1421  final_send_type,
1422  p,
1423  0,
1424  comm_pt()->mpi_comm(),
1425  &req);
1426  classified_send_requests.push_back(req);
1427  MPI_Type_free(&send_types[0]);
1428  MPI_Type_free(&send_types[1]);
1429  MPI_Type_free(&final_send_type);
1430  }
1431  }
1432  }
1433 
1434  // next classify the data as it is received
1435  unsigned n_classified_recv_req = classified_recv_requests.size();
1436  while (n_classified_recv_req > 0)
1437  {
1438  // get the processor number and remove the completed request
1439  // for the vector of requests
1440  int req_num;
1441  MPI_Waitany(n_classified_recv_req,
1442  &classified_recv_requests[0],
1443  &req_num,
1444  MPI_STATUS_IGNORE);
1445  unsigned p = classified_recv_proc[req_num];
1446  classified_recv_requests.erase(classified_recv_requests.begin() +
1447  req_num);
1448  classified_recv_proc.erase(classified_recv_proc.begin() + req_num);
1449  n_classified_recv_req--;
1450 
1451  // next classify the dofs
1452  // and store them for sending to other processors if required
1453  unsigned n_recv = n_classified_recv[p];
1454  for (unsigned i = 0; i < n_recv; i++)
1455  {
1456  v_or_p_values[classified_indices_recv[p][i] - v_or_p_first_row] +=
1457  classified_contributions_recv[p][i];
1458  }
1459 
1460  // clean up
1461  delete[] classified_contributions_recv[p];
1462  delete[] classified_indices_recv[p];
1463  }
1464 
1465  // wait for the unclassified sends to complete
1466  unsigned n_unclassified_send_req = unclassified_send_requests.size();
1467  if (n_unclassified_send_req > 0)
1468  {
1469  MPI_Waitall(n_unclassified_send_req,
1470  &unclassified_send_requests[0],
1471  MPI_STATUS_IGNORE);
1472  }
1473  delete[] unclassified_contributions_send;
1474  delete[] unclassified_indices_send;
1475  delete[] n_unclassified_send;
1476 
1477  // wait for the classified sends to complete
1478  unsigned n_classified_send_req = classified_send_requests.size();
1479  if (n_classified_send_req > 0)
1480  {
1481  MPI_Waitall(n_classified_send_req,
1482  &classified_send_requests[0],
1483  MPI_STATUS_IGNORE);
1484  }
1485  delete[] classified_indices_send;
1486  delete[] classified_contributions_send;
1487  delete[] n_classified_recv;
1488  delete[] n_classified_send;
1489 
1490  // Copy the values back where they belong
1491  if (block_index == 0)
1492  {
1493  v_values = v_or_p_values;
1494  }
1495  else if (block_index == 1)
1496  {
1497  p_values = v_or_p_values;
1498  }
1499  }
1500 
1501 #endif
1502  }
1503  // or if the problem is not distributed
1504  else
1505  {
1506  // find number of elements
1507  unsigned n_el = Navier_stokes_mesh_pt->nelement();
1508 
1509  // Fp needs pressure and velocity mass matrices
1510  unsigned which_one = 0;
1511  if (Use_LSC) which_one = 2;
1512 
1513  // get the contribution for each element
1514  for (unsigned e = 0; e < n_el; e++)
1515  {
1516  // Get element
1517  GeneralisedElement* el_pt = Navier_stokes_mesh_pt->element_pt(e);
1518 
1519  // find number of degrees of freedom in the element
1520  // (this is slightly too big because it includes the
1521  // pressure dofs but this doesn't matter)
1522  unsigned el_dof = el_pt->ndof();
1523 
1524  // allocate local storage for the element's contribution to the
1525  // pressure and velocity mass matrix diagonal
1526  Vector<double> el_vmm_diagonal(el_dof, 0.0);
1527  Vector<double> el_pmm_diagonal(el_dof, 0.0);
1528 
1529  NavierStokesElementWithDiagonalMassMatrices* cast_el_pt = 0;
1530  cast_el_pt =
1531  dynamic_cast<NavierStokesElementWithDiagonalMassMatrices*>(el_pt);
1532  if (cast_el_pt != 0)
1533  {
1534  cast_el_pt->get_pressure_and_velocity_mass_matrix_diagonal(
1535  el_pmm_diagonal, el_vmm_diagonal, which_one);
1536  }
1537  else
1538  {
1539 #ifdef PARANOID
1540  std::ostringstream error_message;
1541  error_message
1542  << "Navier-Stokes mesh contains element that is not of type \n"
1543  << "NavierStokesElementWithDiagonalMassMatrices. \n"
1544  << "The element is in fact of type " << typeid(*el_pt).name()
1545  << "\nWe'll assume that it does not make a used contribution \n"
1546  << "to the inverse diagonal mass matrix used in the "
1547  "preconditioner\n"
1548  << "and (to avoid divisions by zero) fill in dummy unit entries.\n"
1549  << "[This case currently arises with flux control problems\n"
1550  << "where for simplicity the NetFluxControlElement has been added "
1551  "\n"
1552  << "to the Navier Stokes mesh -- this should probably be changed "
1553  "at\n"
1554  << "some point -- if you get this warning in any other context\n"
1555  << "you should check your code very carefully]\n";
1556  OomphLibWarning(error_message.str(),
1557  "NavierStokesSchurComplementPreconditioner::assemble_"
1558  "inv_press_and_veloc_mass_matrix_diagonal()",
1559  OOMPH_EXCEPTION_LOCATION);
1560 #endif
1561 
1562  // Fill in dummy entries to stop division by zero below
1563  for (unsigned j = 0; j < el_dof; j++)
1564  {
1565  el_vmm_diagonal[j] = 1.0;
1566  el_pmm_diagonal[j] = 1.0;
1567  }
1568  }
1569 
1570  // Get the contribution for each dof
1571  for (unsigned i = 0; i < el_dof; i++)
1572  {
1573  // Get the equation number
1574  unsigned eqn_number = el_pt->eqn_number(i);
1575 
1576  // Get the velocity dofs
1577  if (this->block_number(eqn_number) == 0)
1578  {
1579  // get the index in the block
1580  unsigned index = this->index_in_block(eqn_number);
1581 
1582  // if it is required on this processor
1583  if ((index >= v_first_row) &&
1584  (index < (v_first_row + v_nrow_local)))
1585  {
1586  v_values[index - v_first_row] += el_vmm_diagonal[i];
1587  }
1588  }
1589  // Get the pressure dofs
1590  else if (this->block_number(eqn_number) == 1)
1591  {
1592  if (!Use_LSC)
1593  {
1594  // get the index in the block
1595  unsigned index = this->index_in_block(eqn_number);
1596 
1597  // if it is required on this processor
1598  if ((index >= p_first_row) &&
1599  (index < (p_first_row + p_nrow_local)))
1600  {
1601  p_values[index - p_first_row] += el_pmm_diagonal[i];
1602  }
1603  }
1604  }
1605  }
1606  }
1607  }
1608 
1609  // Create column index and row start for velocity mass matrix
1610  int* v_column_index = new int[v_nrow_local];
1611  int* v_row_start = new int[v_nrow_local + 1];
1612  for (unsigned i = 0; i < v_nrow_local; i++)
1613  {
1614 #ifdef PARANOID
1615  if (v_values[i] == 0.0)
1616  {
1617  std::ostringstream error_message;
1618  error_message << "Zero entry in diagonal of velocity mass matrix\n"
1619  << "Index: " << i << std::endl;
1620  throw OomphLibError(error_message.str(),
1621  OOMPH_CURRENT_FUNCTION,
1622  OOMPH_EXCEPTION_LOCATION);
1623  }
1624 #endif
1625  v_values[i] = 1.0 / v_values[i];
1626  v_column_index[i] = v_first_row + i;
1627  v_row_start[i] = i;
1628  }
1629  v_row_start[v_nrow_local] = v_nrow_local;
1630 
1631  // Build the velocity mass matrix
1632  inv_v_mass_pt = new CRDoubleMatrix(this->block_distribution_pt(0));
1633  inv_v_mass_pt->build_without_copy(
1634  v_nrow, v_nrow_local, v_values, v_column_index, v_row_start);
1635 
1636  // Create pressure mass matrix
1637  if (!Use_LSC)
1638  {
1639  // Create column index and row start for pressure mass matrix
1640  int* p_column_index = new int[p_nrow_local];
1641  int* p_row_start = new int[p_nrow_local + 1];
1642  for (unsigned i = 0; i < p_nrow_local; i++)
1643  {
1644 #ifdef PARANOID
1645  if (p_values[i] == 0.0)
1646  {
1647  std::ostringstream error_message;
1648  error_message << "Zero entry in diagonal of pressure mass matrix\n"
1649  << "Index: " << i << std::endl;
1650  throw OomphLibError(error_message.str(),
1651  OOMPH_CURRENT_FUNCTION,
1652  OOMPH_EXCEPTION_LOCATION);
1653  }
1654 #endif
1655  p_values[i] = 1.0 / p_values[i];
1656 
1657  p_column_index[i] = p_first_row + i;
1658  p_row_start[i] = i;
1659  }
1660  p_row_start[p_nrow_local] = p_nrow_local;
1661 
1662  // Build the pressure mass matrix
1663  inv_p_mass_pt = new CRDoubleMatrix(this->block_distribution_pt(1));
1664  inv_p_mass_pt->build_without_copy(
1665  p_nrow, p_nrow_local, p_values, p_column_index, p_row_start);
1666  }
1667  }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
p
float * p
Definition: Tutorial_Map_using.cpp:9
oomph::BlockPreconditioner< CRDoubleMatrix >::master_distribution_pt
const LinearAlgebraDistribution * master_distribution_pt() const
Definition: block_preconditioner.h:2017
oomph::BlockPreconditioner< CRDoubleMatrix >::index_in_block
int index_in_block(const unsigned &i_dof) const
Definition: block_preconditioner.h:1850
oomph::BlockPreconditioner< CRDoubleMatrix >::block_distribution_pt
const LinearAlgebraDistribution * block_distribution_pt(const unsigned &b) const
Access function to the block distributions (const version).
Definition: block_preconditioner.h:1903
oomph::BlockPreconditioner< CRDoubleMatrix >::block_number
int block_number(const unsigned &i_dof) const
Return the block number corresponding to a global index i_dof.
Definition: block_preconditioner.h:1822
oomph::DistributableLinearAlgebraObject::distributed
bool distributed() const
distribution is serial or distributed
Definition: linear_algebra_distribution.h:493
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::LinearAlgebraDistribution::first_row
unsigned first_row() const
Definition: linear_algebra_distribution.h:261
oomph::LinearAlgebraDistribution::nrow
unsigned nrow() const
access function to the number of global rows.
Definition: linear_algebra_distribution.h:186
oomph::LinearAlgebraDistribution::nrow_local
unsigned nrow_local() const
Definition: linear_algebra_distribution.h:193
oomph::Mesh::element_pt
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
oomph::Mesh::nelement
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
oomph::OomphCommunicator::my_rank
int my_rank() const
my rank
Definition: communicator.h:176
oomph::OomphCommunicator::nproc
int nproc() const
number of processors
Definition: communicator.h:157
oomph::Preconditioner::comm_pt
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
Definition: preconditioner.h:171
int
return int(ret)+1
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References oomph::BlockPreconditioner< CRDoubleMatrix >::block_distribution_pt(), oomph::BlockPreconditioner< CRDoubleMatrix >::block_number(), oomph::CRDoubleMatrix::build_without_copy(), oomph::Preconditioner::comm_pt(), oomph::DistributableLinearAlgebraObject::distributed(), e(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::eqn_number(), oomph::LinearAlgebraDistribution::first_row(), oomph::NavierStokesElementWithDiagonalMassMatrices::get_pressure_and_velocity_mass_matrix_diagonal(), i, oomph::BlockPreconditioner< CRDoubleMatrix >::index_in_block(), int(), j, oomph::BlockPreconditioner< CRDoubleMatrix >::master_distribution_pt(), oomph::OomphCommunicator::my_rank(), Navier_stokes_mesh_pt, oomph::GeneralisedElement::ndof(), oomph::Mesh::nelement(), oomph::OomphCommunicator::nproc(), oomph::LinearAlgebraDistribution::nrow(), oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, problem_pt(), Use_LSC, and v.

◆ clean_up_memory()

void oomph::NavierStokesSchurComplementPreconditioner::clean_up_memory ( )
virtual

Helper function to delete preconditioner data.

Reimplemented from oomph::Preconditioner.

1673  {
1674  if (Preconditioner_has_been_setup)
1675  {
1676  // delete matvecs
1677  delete Bt_mat_vec_pt;
1678  Bt_mat_vec_pt = 0;
1679 
1680  delete F_mat_vec_pt;
1681  F_mat_vec_pt = 0;
1682 
1683  delete QBt_mat_vec_pt;
1684  QBt_mat_vec_pt = 0;
1685 
1686  delete E_mat_vec_pt;
1687  E_mat_vec_pt = 0;
1688 
1689  // delete stuff from velocity solve
1690  if (Using_default_f_preconditioner)
1691  {
1692  delete F_preconditioner_pt;
1693  F_preconditioner_pt = 0;
1694  }
1695 
1696  // delete stuff from Schur complement approx
1697  if (Using_default_p_preconditioner)
1698  {
1699  delete P_preconditioner_pt;
1700  P_preconditioner_pt = 0;
1701  }
1702  }
1703  }

References Bt_mat_vec_pt, E_mat_vec_pt, F_mat_vec_pt, F_preconditioner_pt, P_preconditioner_pt, Preconditioner_has_been_setup, QBt_mat_vec_pt, Using_default_f_preconditioner, and Using_default_p_preconditioner.

Referenced by oomph::PseudoElasticFSIPreconditioner::clean_up_memory(), and ~NavierStokesSchurComplementPreconditioner().

◆ disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()

void oomph::NavierStokesSchurComplementPreconditioner::disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements ( )
inline

Do not accept presence of elements that are not of type NavierStokesElementWithDiagonalMassMatrices without issuing a warning.

714  {
715  Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements = false;
716  }

References Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements.

◆ disable_doc_time()

void oomph::NavierStokesSchurComplementPreconditioner::disable_doc_time ( )
inline

Disable documentation of time.

812  {
813  Doc_time = false;
814  }

References Doc_time.

◆ disable_robin_for_fp()

void oomph::NavierStokesSchurComplementPreconditioner::disable_robin_for_fp ( )
inline

Don't use Robin BC elements for the Fp preconditioenr.

827  {
828  Use_robin_for_fp = false;
829  }

References Use_robin_for_fp.

Referenced by FpTestProblem::FpTestProblem().

◆ enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements()

void oomph::NavierStokesSchurComplementPreconditioner::enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements ( )
inline

Accept presence of elements that are not of type NavierStokesElementWithDiagonalMassMatrices without issuing a warning.

706  {
707  Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements = true;
708  }

References Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements.

◆ enable_doc_time()

void oomph::NavierStokesSchurComplementPreconditioner::enable_doc_time ( )
inline

Enable documentation of time.

806  {
807  Doc_time = true;
808  }

References Doc_time.

Referenced by main().

◆ enable_robin_for_fp()

void oomph::NavierStokesSchurComplementPreconditioner::enable_robin_for_fp ( )
inline

Use Robin BC elements for the Fp preconditioner.

821  {
822  Use_robin_for_fp = true;
823  }

References Use_robin_for_fp.

Referenced by FpTestProblem::FpTestProblem().

◆ get_pressure_advection_diffusion_matrix()

void oomph::NavierStokesSchurComplementPreconditioner::get_pressure_advection_diffusion_matrix ( CRDoubleMatrix &  fp_matrix)
inline

Get the pressure advection diffusion matrix.

881  {
882  // Pin all non-pressure dofs
883  pin_all_non_pressure_dofs();
884 
885  // Get the spatial dimension
886  unsigned ndim = Navier_stokes_mesh_pt->finite_element_pt(0)->dim();
887 
888  // Backup assembly handler and set new one
889  AssemblyHandler* backed_up_assembly_handler_pt =
890  problem_pt()->assembly_handler_pt();
891  problem_pt()->assembly_handler_pt() =
892  new FpPreconditionerAssemblyHandler(ndim);
893 
894  // Backup submeshes (if any)
895  unsigned n_sub_mesh = problem_pt()->nsub_mesh();
896  Vector<Mesh*> backed_up_sub_mesh_pt(n_sub_mesh);
897  for (unsigned i = 0; i < n_sub_mesh; i++)
898  {
899  backed_up_sub_mesh_pt[i] = problem_pt()->mesh_pt(i);
900  }
901  // Flush sub-meshes but don't call rebuild_global_mesh()
902  // so we can simply re-instate the sub-meshes below.
903  problem_pt()->flush_sub_meshes();
904 
905  // Backup the problem's own mesh pointer
906  Mesh* backed_up_mesh_pt = problem_pt()->mesh_pt();
907  problem_pt()->mesh_pt() = Navier_stokes_mesh_pt;
908 
909 #ifdef OOMPH_HAS_MPI
910 
911  // Backup the start and end elements for the distributed
912  // assembly process
913  Vector<unsigned> backed_up_first_el_for_assembly;
914  Vector<unsigned> backed_up_last_el_for_assembly;
915  problem_pt()->get_first_and_last_element_for_assembly(
916  backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
917 
918  // Now re-assign
919  problem_pt()->set_default_first_and_last_element_for_assembly();
920 
921 #endif
922 
923  // Identify pinned pressure dof
924  int pinned_pressure_eqn = -2;
925  if (Pin_first_pressure_dof_in_press_adv_diff)
926  {
927  // Loop over all Navier Stokes elements to find first non-pinned
928  // pressure dof
929  unsigned nel = Navier_stokes_mesh_pt->nelement();
930  for (unsigned e = 0; e < nel; e++)
931  {
932  TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
933  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
934  Navier_stokes_mesh_pt->element_pt(e));
935  int local_eqn = bulk_elem_pt->p_local_eqn(0);
936  if (local_eqn >= 0)
937  {
938  pinned_pressure_eqn = bulk_elem_pt->eqn_number(local_eqn);
939  break;
940  }
941  }
942 
943 
944 #ifdef OOMPH_HAS_MPI
945  if (problem_pt()->distributed())
946  {
947  int pinned_pressure_eqn_local = pinned_pressure_eqn;
948  int pinned_pressure_eqn_global = pinned_pressure_eqn;
949  MPI_Allreduce(&pinned_pressure_eqn_local,
950  &pinned_pressure_eqn_global,
951  1,
952  MPI_INT,
953  MPI_MIN,
954  this->comm_pt()->mpi_comm());
955  pinned_pressure_eqn = pinned_pressure_eqn_global;
956  }
957 
958 #endif
959  }
960 
961 
962  // Loop over all Navier Stokes elements
963  unsigned nel = Navier_stokes_mesh_pt->nelement();
964  for (unsigned e = 0; e < nel; e++)
965  {
966  TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
967  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
968  Navier_stokes_mesh_pt->element_pt(e));
969 
970  // Set pinned pressure equation
971  bulk_elem_pt->pinned_fp_pressure_eqn() = pinned_pressure_eqn;
972  }
973 
974 
975  // Attach Robin BC elements
976  unsigned nbound = Navier_stokes_mesh_pt->nboundary();
977  if (Use_robin_for_fp)
978  {
979  // Loop over all boundaries of Navier Stokes mesh
980  for (unsigned b = 0; b < nbound; b++)
981  {
982  // How many bulk elements are adjacent to boundary b?
983  unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
984 
985  // Loop over the bulk elements adjacent to boundary b?
986  for (unsigned e = 0; e < n_element; e++)
987  {
988  TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
989  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
990  Navier_stokes_mesh_pt->boundary_element_pt(b, e));
991 
992  // What is the index of the face of the bulk element e on bondary b
993  int face_index =
994  Navier_stokes_mesh_pt->face_index_at_boundary(b, e);
995 
996  // Build face element
997  bulk_elem_pt->build_fp_press_adv_diff_robin_bc_element(face_index);
998 
999  } // end of loop over bulk elements adjacent to boundary b
1000  }
1001  }
1002 
1003  // Get "Jacobian" of the modified system
1004  DoubleVector dummy_residuals;
1005  problem_pt()->get_jacobian(dummy_residuals, fp_matrix);
1006 
1007  // Kill Robin BC elements
1008  if (Use_robin_for_fp)
1009  {
1010  // Loop over all boundaries of Navier Stokes mesh
1011  for (unsigned b = 0; b < nbound; b++)
1012  {
1013  // How many bulk elements are adjacent to boundary b?
1014  unsigned n_element = Navier_stokes_mesh_pt->nboundary_element(b);
1015 
1016  // Loop over the bulk elements adjacent to boundary b?
1017  for (unsigned e = 0; e < n_element; e++)
1018  {
1019  TemplateFreeNavierStokesEquationsBase* bulk_elem_pt =
1020  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
1021  Navier_stokes_mesh_pt->boundary_element_pt(b, e));
1022 
1023  // Kill associated Robin elements
1024  bulk_elem_pt->delete_pressure_advection_diffusion_robin_elements();
1025 
1026  } // end of loop over bulk elements adjacent to boundary b
1027  }
1028  }
1029 
1030  // Reset pin status
1031  reset_pin_status();
1032 
1033 #ifdef OOMPH_HAS_MPI
1034 
1035  // Reset start and end elements for the distributed
1036  // assembly process
1037  problem_pt()->set_first_and_last_element_for_assembly(
1038  backed_up_first_el_for_assembly, backed_up_last_el_for_assembly);
1039 
1040 #endif
1041 
1042  // Cleanup and reset assembly handler
1043  delete problem_pt()->assembly_handler_pt();
1044  problem_pt()->assembly_handler_pt() = backed_up_assembly_handler_pt;
1045 
1046  // Re-instate submeshes. (No need to call rebuild_global_mesh()
1047  // as it was never unbuilt).
1048  for (unsigned i = 0; i < n_sub_mesh; i++)
1049  {
1050  problem_pt()->add_sub_mesh(backed_up_sub_mesh_pt[i]);
1051  }
1052 
1053 
1054  // Reset the problem's mesh pointer
1055  problem_pt()->mesh_pt() = backed_up_mesh_pt;
1056  }
DoubleVector
std::vector< double > DoubleVector
loads clump configuration
Definition: ClumpInput.h:26
b
Scalar * b
Definition: benchVecAdd.cpp:17
oomph::FiniteElement::dim
unsigned dim() const
Definition: elements.h:2611
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::face_index_at_boundary
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
Definition: mesh.h:896
oomph::Mesh::nboundary_element
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
oomph::Mesh::nboundary
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
oomph::Mesh::boundary_element_pt
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
oomph::NavierStokesSchurComplementPreconditioner::pin_all_non_pressure_dofs
void pin_all_non_pressure_dofs()
Pin all non-pressure dofs.
Definition: navier_stokes_preconditioners.h:862
oomph::NavierStokesSchurComplementPreconditioner::reset_pin_status
void reset_pin_status()
Reset pin status of all values.
Definition: navier_stokes_preconditioners.h:1060
oomph::Problem::add_sub_mesh
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Definition: problem.h:1330
oomph::Problem::assembly_handler_pt
AssemblyHandler *& assembly_handler_pt()
Return a pointer to the assembly handler object.
Definition: problem.h:1570
oomph::Problem::flush_sub_meshes
void flush_sub_meshes()
Definition: problem.h:1339
oomph::Problem::mesh_pt
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
oomph::Problem::nsub_mesh
unsigned nsub_mesh() const
Return number of submeshes.
Definition: problem.h:1323
oomph::Problem::distributed
bool distributed() const
Definition: problem.h:808
oomph::Problem::get_jacobian
virtual void get_jacobian(DoubleVector &residuals, DenseDoubleMatrix &jacobian)
Definition: problem.cc:3890

References oomph::Problem::add_sub_mesh(), oomph::Problem::assembly_handler_pt(), b, oomph::Mesh::boundary_element_pt(), oomph::TemplateFreeNavierStokesEquationsBase::build_fp_press_adv_diff_robin_bc_element(), oomph::Preconditioner::comm_pt(), oomph::TemplateFreeNavierStokesEquationsBase::delete_pressure_advection_diffusion_robin_elements(), oomph::FiniteElement::dim(), oomph::Problem::distributed(), e(), oomph::Mesh::element_pt(), oomph::GeneralisedElement::eqn_number(), oomph::Mesh::face_index_at_boundary(), oomph::Mesh::finite_element_pt(), oomph::Problem::flush_sub_meshes(), oomph::Problem::get_jacobian(), i, oomph::Problem::mesh_pt(), Navier_stokes_mesh_pt, oomph::Mesh::nboundary(), oomph::Mesh::nboundary_element(), oomph::Mesh::nelement(), oomph::Problem::nsub_mesh(), oomph::TemplateFreeNavierStokesEquationsBase::p_local_eqn(), pin_all_non_pressure_dofs(), Pin_first_pressure_dof_in_press_adv_diff, oomph::TemplateFreeNavierStokesEquationsBase::pinned_fp_pressure_eqn(), problem_pt(), reset_pin_status(), and Use_robin_for_fp.

◆ pin_all_non_pressure_dofs()

void oomph::NavierStokesSchurComplementPreconditioner::pin_all_non_pressure_dofs ( )
inline

Pin all non-pressure dofs.

863  {
864  // Backup pin status and then pin all non-pressure degrees of freedom
865  unsigned nelem = Navier_stokes_mesh_pt->nelement();
866  for (unsigned e = 0; e < nelem; e++)
867  {
868  // Get pointer to the bulk element that is adjacent to boundary b
869  TemplateFreeNavierStokesEquationsBase* el_pt =
870  dynamic_cast<TemplateFreeNavierStokesEquationsBase*>(
871  Navier_stokes_mesh_pt->element_pt(e));
872 
873  // Pin
874  el_pt->pin_all_non_pressure_dofs(Eqn_number_backup);
875  }
876  }
oomph::NavierStokesSchurComplementPreconditioner::Eqn_number_backup
std::map< Data *, std::vector< int > > Eqn_number_backup
Definition: navier_stokes_preconditioners.h:1185

References e(), oomph::Mesh::element_pt(), Eqn_number_backup, Navier_stokes_mesh_pt, oomph::Mesh::nelement(), and oomph::TemplateFreeNavierStokesEquationsBase::pin_all_non_pressure_dofs().

Referenced by get_pressure_advection_diffusion_matrix().

◆ pin_first_pressure_dof_in_press_adv_diff()

void oomph::NavierStokesSchurComplementPreconditioner::pin_first_pressure_dof_in_press_adv_diff ( )
inline

Set boolean indicating that we want to pin first pressure dof in Navier Stokes mesh when assembling the pressure advection diffusion system for Fp preconditoner – needed at zero Reynolds number for non-enclosed flows but seems harmless in any case

837  {
838  Pin_first_pressure_dof_in_press_adv_diff = true;
839  }

References Pin_first_pressure_dof_in_press_adv_diff.

◆ preconditioner_solve()

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

Apply preconditioner to Vector r.

Apply preconditioner to r.

Implements oomph::Preconditioner.

385  {
386 
387 #ifdef PARANOID
388  if (Preconditioner_has_been_setup==false)
389  {
390  std::ostringstream error_message;
391  error_message << "setup must be called before using preconditioner_solve";
392  throw OomphLibError(
393  error_message.str(),
394  OOMPH_CURRENT_FUNCTION,
395  OOMPH_EXCEPTION_LOCATION);
396  }
397  if (z.built())
398  {
399  if (z.nrow() != r.nrow())
400  {
401  std::ostringstream error_message;
402  error_message << "The vectors z and r must have the same number of "
403  << "of global rows";
404  throw OomphLibError(
405  error_message.str(),
406  OOMPH_CURRENT_FUNCTION,
407  OOMPH_EXCEPTION_LOCATION);
408  }
409  }
410 #endif
411 
412  // if z is not setup then give it the same distribution
413  if (!z.distribution_pt()->built())
414  {
415  z.build(r.distribution_pt(),0.0);
416  }
417 
418  // Step 1 - apply approximate Schur inverse to pressure unknowns (block 1)
419  // -----------------------------------------------------------------------
420 
421  // Working vectors
422  DoubleVector temp_vec;
423  DoubleVector another_temp_vec;
424  DoubleVector yet_another_temp_vec;
425 
426  // Copy pressure values from residual vector to temp_vec:
427  // Loop over all entries in the global vector (this one
428  // includes velocity and pressure dofs in some random fashion)
429  this->get_block_vector(1,r,temp_vec);
430 
431  // NOTE: The vector temp_vec now contains the vector r_p.
432 
433 
434  // LSC version
435  if (Use_LSC)
436  {
437 
438  // Solve first pressure Poisson system
439 #ifdef PARANOID
440  // check a solver has been set
441  if (P_preconditioner_pt==0)
442  {
443  std::ostringstream error_message;
444  error_message << "P_preconditioner_pt has not been set.";
445  throw OomphLibError(
446  error_message.str(),
447  OOMPH_CURRENT_FUNCTION,
448  OOMPH_EXCEPTION_LOCATION);
449  }
450 #endif
451 
452  // use some Preconditioner's preconditioner_solve function
453  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
454 
455  // NOTE: The vector another_temp_vec now contains the vector P^{-1} r_p
456 
457  // Multiply another_temp_vec by matrix E and stick the result into temp_vec
458  temp_vec.clear();
459  QBt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
460  another_temp_vec.clear();
461  F_mat_vec_pt->multiply(temp_vec,another_temp_vec);
462  temp_vec.clear();
463  QBt_mat_vec_pt->multiply_transpose(another_temp_vec, temp_vec);
464 
465 
466  // NOTE: The vector temp_vec now contains E P^{-1} r_p
467 
468  // Solve second pressure Poisson system using preconditioner_solve
469  another_temp_vec.clear();
470  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
471 
472  // NOTE: The vector another_temp_vec now contains z_p = P^{-1} E P^{-1} r_p
473  // as required (apart from the sign which we'll fix in the
474  // next step.
475 
476  }
477  // Fp version
478  else
479  {
480 
481  // Multiply temp_vec by matrix E and stick the result into
482  // yet_another_temp_vec
483  E_mat_vec_pt->multiply(temp_vec,yet_another_temp_vec);
484 
485  // NOTE: The vector yet_another_temp_vec now contains Fp Qp^{-1} r_p
486 
487  // Solve pressure Poisson system
488 #ifdef PARANOID
489  // check a solver has been set
490  if (P_preconditioner_pt==0)
491  {
492  std::ostringstream error_message;
493  error_message << "P_preconditioner_pt has not been set.";
494  throw OomphLibError(
495  error_message.str(),
496  OOMPH_CURRENT_FUNCTION,
497  OOMPH_EXCEPTION_LOCATION);
498  }
499 #endif
500 
501  // Solve second pressure Poisson system using preconditioner_solve
502  another_temp_vec.clear();
503  P_preconditioner_pt->preconditioner_solve(yet_another_temp_vec,
504  another_temp_vec);
505 
506  // NOTE: The vector another_temp_vec now contains
507  // z_p = P^{-1} Fp Qp^{-1} r_p
508  // as required (apart from the sign which we'll fix in the
509  // next step.
510 
511  }
512 
513  // Now copy another_temp_vec (i.e. z_p) back into the global vector z.
514  // Loop over all entries in the global results vector z:
515  temp_vec.build(another_temp_vec.distribution_pt(),0.0);
516  temp_vec -= another_temp_vec;
517  return_block_vector(1,temp_vec,z);
518 
519 
520 
521  // Step 2 - apply preconditioner to velocity unknowns (block 0)
522  // ------------------------------------------------------------
523 
524  // Recall that another_temp_vec (computed above) contains the
525  // negative of the solution of the Schur complement systen, -z_p.
526  // Multiply by G (stored in Block_matrix_pt(0,1) and store
527  // result in temp_vec (vector resizes itself).
528  temp_vec.clear();
529  Bt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
530 
531  // NOTE: temp_vec now contains -G z_p
532 
533  // The vector another_temp_vec is no longer needed -- re-use it to store
534  // velocity quantities:
535  another_temp_vec.clear();
536 
537  // Loop over all enries in the global vector and find the
538  // entries associated with the velocities:
539  get_block_vector(0,r,another_temp_vec);
540  another_temp_vec += temp_vec;
541 
542  // NOTE: The vector another_temp_vec now contains r_u - G z_p
543 
544  // Solve momentum system
545 #ifdef PARANOID
546  // check a solver has been set
547  if (F_preconditioner_pt==0)
548  {
549  std::ostringstream error_message;
550  error_message << "F_preconditioner_pt has not been set.";
551  throw OomphLibError(
552  error_message.str(),
553  OOMPH_CURRENT_FUNCTION,
554  OOMPH_EXCEPTION_LOCATION);
555  }
556 #endif
557 
558  // use some Preconditioner's preconditioner solve
559  // and return
560  if (F_preconditioner_is_block_preconditioner)
561  {
562  return_block_vector(0,another_temp_vec,z);
563  F_preconditioner_pt->preconditioner_solve(z,z);
564  }
565  else
566  {
567  F_preconditioner_pt->preconditioner_solve(another_temp_vec, temp_vec);
568  return_block_vector(0,temp_vec,z);
569  }
570  }
oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vector
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Definition: block_preconditioner.cc:4489
oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vector
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Definition: block_preconditioner.cc:4186
oomph::DistributableLinearAlgebraObject::distribution_pt
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Definition: linear_algebra_distribution.h:457
oomph::DistributableLinearAlgebraObject::nrow
unsigned nrow() const
access function to the number of global rows.
Definition: linear_algebra_distribution.h:463
oomph::DoubleVector
Definition: double_vector.h:58
oomph::DoubleVector::build
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
Definition: double_vector.cc:35
oomph::DoubleVector::built
bool built() const
Definition: double_vector.h:154
oomph::DoubleVector::clear
void clear()
wipes the DoubleVector
Definition: double_vector.h:142
oomph::LinearAlgebraDistribution::built
bool built() const
Definition: linear_algebra_distribution.h:342
oomph::MatrixVectorProduct::multiply_transpose
void multiply_transpose(const DoubleVector &x, DoubleVector &y) const
Definition: matrix_vector_product.cc:177
oomph::MatrixVectorProduct::multiply
void multiply(const DoubleVector &x, DoubleVector &y) const
Definition: matrix_vector_product.cc:108
oomph::NavierStokesSchurComplementPreconditioner::F_preconditioner_is_block_preconditioner
bool F_preconditioner_is_block_preconditioner
Definition: navier_stokes_preconditioners.h:1152
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::Preconditioner::preconditioner_solve
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20

References oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::built(), oomph::DoubleVector::clear(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and UniformPSDSelfTest::r.

Referenced by oomph::FSIPreconditioner::preconditioner_solve(), and oomph::PseudoElasticFSIPreconditioner::preconditioner_solve().

◆ problem_pt()

Problem* oomph::NavierStokesSchurComplementPreconditioner::problem_pt ( ) const
inline
690  {
691 #ifdef PARANOID
692  if (Problem_pt == 0)
693  {
694  std::ostringstream error_msg;
695  error_msg << "Problem pointer is null, did you set it yet?";
696  throw OomphLibError(
697  error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
698  }
699 #endif
700  return Problem_pt;
701  }

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and Problem_pt.

Referenced by assemble_inv_press_and_veloc_mass_matrix_diagonal(), get_pressure_advection_diffusion_matrix(), and set_problem_pt().

◆ reset_pin_status()

void oomph::NavierStokesSchurComplementPreconditioner::reset_pin_status ( )
inline

Reset pin status of all values.

1061  {
1062  // Reset pin status for nodes
1063  unsigned nnod = Navier_stokes_mesh_pt->nnode();
1064  for (unsigned j = 0; j < nnod; j++)
1065  {
1066  Node* nod_pt = Navier_stokes_mesh_pt->node_pt(j);
1067  unsigned nval = nod_pt->nvalue();
1068  for (unsigned i = 0; i < nval; i++)
1069  {
1070  nod_pt->eqn_number(i) = Eqn_number_backup[nod_pt][i];
1071  }
1072 
1073  // If it's a solid node deal with its positional data too
1074  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1075  if (solid_nod_pt != 0)
1076  {
1077  Data* solid_posn_data_pt = solid_nod_pt->variable_position_pt();
1078  unsigned nval = solid_posn_data_pt->nvalue();
1079  for (unsigned i = 0; i < nval; i++)
1080  {
1081  solid_posn_data_pt->eqn_number(i) =
1082  Eqn_number_backup[solid_posn_data_pt][i];
1083  }
1084  }
1085  }
1086 
1087  // ... and internal data
1088  unsigned nelem = Navier_stokes_mesh_pt->nelement();
1089  for (unsigned e = 0; e < nelem; e++)
1090  {
1091  // Pointer to element
1092  GeneralisedElement* el_pt = Navier_stokes_mesh_pt->element_pt(e);
1093 
1094  // Pin/unpin internal data
1095  unsigned nint = el_pt->ninternal_data();
1096  for (unsigned j = 0; j < nint; j++)
1097  {
1098  Data* data_pt = el_pt->internal_data_pt(j);
1099  unsigned nvalue = data_pt->nvalue();
1100  for (unsigned i = 0; i < nvalue; i++)
1101  {
1102  data_pt->eqn_number(i) = Eqn_number_backup[data_pt][i];
1103  }
1104  }
1105  }
1106 
1107  // Free up storage
1108  Eqn_number_backup.clear();
1109  }
oomph::Data::nvalue
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
oomph::GeneralisedElement::ninternal_data
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
oomph::Mesh::nnode
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436

References e(), oomph::Mesh::element_pt(), oomph::Data::eqn_number(), Eqn_number_backup, i, oomph::GeneralisedElement::internal_data_pt(), j, Navier_stokes_mesh_pt, oomph::Mesh::nelement(), oomph::GeneralisedElement::ninternal_data(), oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), and oomph::SolidNode::variable_position_pt().

Referenced by get_pressure_advection_diffusion_matrix().

◆ set_f_preconditioner()

void oomph::NavierStokesSchurComplementPreconditioner::set_f_preconditioner ( Preconditioner *  new_f_preconditioner_pt)
inline

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

770  {
771  // If the default preconditioner has been used
772  // clean it up now...
773  if (Using_default_f_preconditioner)
774  {
775  delete F_preconditioner_pt;
776  }
777  F_preconditioner_pt = new_f_preconditioner_pt;
778  Using_default_f_preconditioner = false;
779  }

References F_preconditioner_pt, and Using_default_f_preconditioner.

Referenced by FpTestProblem::FpTestProblem(), main(), FSIChannelWithLeafletProblem< ELEMENT >::set_iterative_solver(), PseudoElasticCollapsibleChannelProblem< FLUID_ELEMENT, SOLID_ELEMENT >::set_pseudo_elastic_fsi_solver(), set_pseudo_elastic_fsi_solver(), TiltedCavityProblem< ELEMENT >::TiltedCavityProblem(), and TorusProblem< ELEMENT >::TorusProblem().

◆ set_f_superlu_preconditioner()

void oomph::NavierStokesSchurComplementPreconditioner::set_f_superlu_preconditioner ( )
inline

Function to (re-)set momentum matrix preconditioner (inexact solver) to SuperLU

796  {
797  if (!Using_default_f_preconditioner)
798  {
799  F_preconditioner_pt = new SuperLUPreconditioner;
800  Using_default_f_preconditioner = true;
801  }
802  }

References F_preconditioner_pt, and Using_default_f_preconditioner.

◆ set_navier_stokes_mesh()

void oomph::NavierStokesSchurComplementPreconditioner::set_navier_stokes_mesh ( Mesh *  mesh_pt,
const bool &  allow_multiple_element_type_in_navier_stokes_mesh = false 
)
inline

Specify the mesh containing the block-preconditionable Navier-Stokes elements. The optional argument indicates if there are multiple types of elements in the same mesh.

735  {
736  // Store the mesh pointer.
737  Navier_stokes_mesh_pt = mesh_pt;
738 
739  // Are there multiple element types in this mesh?
740  Allow_multiple_element_type_in_navier_stokes_mesh =
741  allow_multiple_element_type_in_navier_stokes_mesh;
742  }
oomph::BlockPreconditioner< CRDoubleMatrix >::mesh_pt
const Mesh * mesh_pt(const unsigned &i) const
Definition: block_preconditioner.h:1782

References Allow_multiple_element_type_in_navier_stokes_mesh, oomph::BlockPreconditioner< CRDoubleMatrix >::mesh_pt(), and Navier_stokes_mesh_pt.

Referenced by FpTestProblem::FpTestProblem(), main(), oomph::FSIPreconditioner::setup(), oomph::PseudoElasticFSIPreconditioner::setup(), TiltedCavityProblem< ELEMENT >::TiltedCavityProblem(), and TorusProblem< ELEMENT >::TorusProblem().

◆ set_p_preconditioner()

void oomph::NavierStokesSchurComplementPreconditioner::set_p_preconditioner ( Preconditioner *  new_p_preconditioner_pt)
inline

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

746  {
747  // If the default preconditioner has been used
748  // clean it up now...
749  if (Using_default_p_preconditioner)
750  {
751  delete P_preconditioner_pt;
752  }
753  P_preconditioner_pt = new_p_preconditioner_pt;
754  Using_default_p_preconditioner = false;
755  }

References P_preconditioner_pt, and Using_default_p_preconditioner.

Referenced by FpTestProblem::FpTestProblem(), main(), FSIChannelWithLeafletProblem< ELEMENT >::set_iterative_solver(), PseudoElasticCollapsibleChannelProblem< FLUID_ELEMENT, SOLID_ELEMENT >::set_pseudo_elastic_fsi_solver(), set_pseudo_elastic_fsi_solver(), TiltedCavityProblem< ELEMENT >::TiltedCavityProblem(), and TorusProblem< ELEMENT >::TorusProblem().

◆ set_p_superlu_preconditioner()

void oomph::NavierStokesSchurComplementPreconditioner::set_p_superlu_preconditioner ( )
inline

Function to (re-)set pressure matrix preconditioner (inexact solver) to SuperLU

760  {
761  if (!Using_default_p_preconditioner)
762  {
763  P_preconditioner_pt = new SuperLUPreconditioner;
764  Using_default_p_preconditioner = true;
765  }
766  }

References P_preconditioner_pt, and Using_default_p_preconditioner.

◆ set_problem_pt()

void oomph::NavierStokesSchurComplementPreconditioner::set_problem_pt ( Problem *  problem_pt)
inline

Broken assignment operator.

Set the problem pointer (non-const pointer, the problem WILL be changed) for use in get_pressure_advection_diffusion_matrix().

685  {
686  Problem_pt = problem_pt;
687  }

References problem_pt(), and Problem_pt.

◆ setup() [1/4]

void oomph::NavierStokesSchurComplementPreconditioner::setup ( )
virtual

Setup the preconditioner.

Setup the least-squares commutator Navier Stokes preconditioner. This extracts blocks corresponding to the velocity and pressure 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.

///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////// Setup the least-squares commutator Navier Stokes preconditioner. This extracts blocks corresponding to the velocity and pressure 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.

168  {
169 
170  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
171  // NOTE: In the interest of minimising memory usage, several containers
172  // are recycled, therefore their content/meaning changes
173  // throughout this function. The code is carefully annotated
174  // but you'll have to read it line by line!
175  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
176 
177  // make sure any old data is deleted
178  clean_up_memory(); // hierher
179 
180 #ifdef PARANOID
181  // paranoid check that the navier stokes mesh pt has been set
182  if (Navier_stokes_mesh_pt == 0)
183  {
184  std::ostringstream error_message;
185  error_message << "The navier stokes elements mesh pointer must be set.\n"
186  << "Use method set_navier_stokes_mesh(...)";
187  throw OomphLibError(error_message.str(),
188  OOMPH_CURRENT_FUNCTION,
189  OOMPH_EXCEPTION_LOCATION);
190  }
191 #endif
192 
193  // set the mesh
194  this->set_nmesh(1);
195  this->set_mesh(0,Navier_stokes_mesh_pt);
196 
197  // Get blocks
198  // ----------
199 
200  // In comes the current Jacobian. Recast it to a CR double matrix;
201  // shout if that can't be done.
202  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
203 
204 
205 #ifdef PARANOID
206  if (cr_matrix_pt==0)
207  {
208  std::ostringstream error_message;
209  error_message
210  << "NavierStokesSchurComplementPreconditioner only works with "
211  << "CRDoubleMatrix matrices" << std::endl;
212  throw OomphLibError(error_message.str(),
213  OOMPH_CURRENT_FUNCTION,
214  OOMPH_EXCEPTION_LOCATION);
215  }
216 #endif
217 
218  // Set up block look up schemes
219  //-----------------------------
220 
221  // this preconditioner has two types of block:
222  // type 0: velocity - corresponding to DOFs 0 to n-2
223  // type 1: pressure - corresponding to DOF n-1
224  unsigned ndof_types = 0;
225  if (this->is_subsidiary_block_preconditioner())
226  {
227  ndof_types = this->ndof_types();
228  }
229  else
230  {
231  ndof_types = this->ndof_types_in_mesh(0);
232  }
233 
234  // hierher
235  Vector<unsigned> dof_to_block_map(ndof_types);
236  dof_to_block_map[ndof_types-1]=1;
237  this->block_setup(dof_to_block_map);
238 
239  // Get B (the divergence block)
240  CRDoubleMatrix* b_pt = 0;
241  this->get_block(1,0,b_pt);
242 
243  // // get the inverse velocity and pressure mass matrices
244  // CRDoubleMatrix* inv_v_mass_pt = 0;
245  // CRDoubleMatrix* inv_p_mass_pt = 0;
246  // double ivmm_assembly_start_t = TimingHelpers::timer();
247  // if (Use_LSC)
248  // {
249  // // We only need the velocity mass matrix
250  // bool do_both=false;
251  // assemble_inv_press_and_veloc_mass_matrix_diagonal(inv_p_mass_pt,
252  // inv_v_mass_pt,
253  // do_both);
254  // }
255  // else
256  // {
257  // // We only need both mass matrices
258  // bool do_both=true;
259  // assemble_inv_press_and_veloc_mass_matrix_diagonal(inv_p_mass_pt,
260  // inv_v_mass_pt,
261  // do_both);
262  // }
263 
264 
265  // Get gradient matrix Bt
266  CRDoubleMatrix* bt_pt = 0;
267  this->get_block(0,1,bt_pt);
268 
269  // Multiply inverse velocity mass matrix by gradient matrix B^T
270  double t_QBt_matrix_start = TimingHelpers::timer();
271  CRDoubleMatrix* qbt_pt = new CRDoubleMatrix;
272  inv_v_mass_pt->multiply(*bt_pt, *qbt_pt);
273  delete bt_pt;
274 
275  // Store product in bt_pt
276  bt_pt = qbt_pt;
277  delete inv_v_mass_pt;
278 
279  // Multiply B from left by divergence matrix B and store result in
280  // pressure Poisson matrix.
281  b_pt->multiply(*bt_pt, *p_matrix_pt);
282 
283  // Kill divergence matrix because we don't need it any more
284  delete b_pt;
285 
286  // Build the matvec operator for QBt
287  QBt_mat_vec_pt = new MatrixVectorProduct;
288  //QBt_mat_vec_pt->setup(bt_pt);
289  this->setup_matrix_vector_product(QBt_mat_vec_pt,
290  bt_pt, 1);
291 
292  // Kill gradient matrix B^T (it's been overwritten anyway and
293  // needs to be recomputed afresh below)
294  delete bt_pt;
295 
296  // Get momentum block F
297  CRDoubleMatrix* f_pt = 0;
298  this->get_block(0,0,f_pt);
299 
300  // form the matrix vector product helper
301  F_mat_vec_pt = new MatrixVectorProduct;
302  //F_mat_vec_pt->setup(f_pt);
303  this->setup_matrix_vector_product(F_mat_vec_pt,f_pt,0);
304 
305  // if F is a block preconditioner then we can delete the F matrix
306  if (F_preconditioner_is_block_preconditioner)
307  {
308  delete f_pt;
309  }
310 
311  // Rebuild Bt (remember that we temporarily overwrote
312  // it by its product with the inverse velocity mass matrix)
313  bt_pt = 0;
314  this->get_block(0,1,bt_pt);
315 
316  // Form the matrix vector operator for Bt
317  Bt_mat_vec_pt = new MatrixVectorProduct;
318  //Bt_mat_vec_pt->setup(bt_pt);
319  this->setup_matrix_vector_product(Bt_mat_vec_pt,bt_pt,1);
320  delete bt_pt;
321 
322  // If the P preconditioner has not been setup
323  P_preconditioner_pt = new SuperLUPreconditioner;
324  P_preconditioner_pt->setup(p_matrix_pt,comm_pt());
325  delete p_matrix_pt;
326  p_matrix_pt=0;
327 
328  // Set up solver for solution of system with momentum matrix
329  // ----------------------------------------------------------
330 
331  F_preconditioner_pt = new SuperLUPreconditioner;
332  Using_default_f_preconditioner = true;
333 
334 
335  // if F is a block preconditioner
336  double t_f_prec_start = TimingHelpers::timer();
337  if (F_preconditioner_is_block_preconditioner)
338  {
339  unsigned ndof_types = this->ndof_types();
340  ndof_types--;
341  Vector<unsigned> dof_map(ndof_types);
342  for (unsigned i = 0; i < ndof_types; i++)
343  {
344  dof_map[i] = i;
345  }
346  F_block_preconditioner_pt->
347  turn_into_subsidiary_block_preconditioner(this,dof_map);
348 
349  // Set the mesh in the subsidiary preconditioner.
350  // RAYRAY This is incorrect. We never set the mesh for a subsidiary
351  // block preconditioner. The self test still passes since this part of
352  // the if-else condition is never executed.
353  F_block_preconditioner_pt->set_nmesh(1);
354  F_block_preconditioner_pt->set_mesh(0, Navier_stokes_mesh_pt);
355 
356  F_block_preconditioner_pt->setup(matrix_pt(),comm_pt());
357  }
358  // otherwise F is not a block preconditioner
359  else
360  {
361  F_preconditioner_pt->setup(f_pt,comm_pt());
362  delete f_pt;
363  }
364  double t_f_prec_finish = TimingHelpers::timer();
365  if(Doc_time)
366  {
367  double t_f_prec_time = t_f_prec_finish - t_f_prec_start;
368  oomph_info << "F sub-preconditioner setup time [sec]: "
369  << t_f_prec_time << "\n";
370  }
371 
372  // Remember that the preconditioner has been setup so
373  // the stored information can be wiped when we
374  // come here next...
375  Preconditioner_has_been_setup = true;
376  }
oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types_in_mesh
unsigned ndof_types_in_mesh(const unsigned &i) const
Definition: block_preconditioner.h:2037
oomph::BlockPreconditioner< CRDoubleMatrix >::get_block
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Definition: block_preconditioner.h:673
oomph::BlockPreconditioner< CRDoubleMatrix >::matrix_pt
CRDoubleMatrix * matrix_pt() const
Definition: block_preconditioner.h:520
oomph::BlockPreconditioner< CRDoubleMatrix >::is_subsidiary_block_preconditioner
bool is_subsidiary_block_preconditioner() const
Definition: block_preconditioner.h:2061
oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types
unsigned ndof_types() const
Return the total number of DOF types.
Definition: block_preconditioner.h:1694
oomph::BlockPreconditioner< CRDoubleMatrix >::turn_into_subsidiary_block_preconditioner
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Definition: block_preconditioner.cc:2376
oomph::BlockPreconditioner< CRDoubleMatrix >::setup_matrix_vector_product
void setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
Definition: block_preconditioner.h:2397
oomph::BlockPreconditioner< CRDoubleMatrix >::set_nmesh
void set_nmesh(const unsigned &n)
Definition: block_preconditioner.h:2851
oomph::BlockPreconditioner< CRDoubleMatrix >::block_setup
virtual void block_setup()
Definition: block_preconditioner.cc:2483
oomph::BlockPreconditioner< CRDoubleMatrix >::set_mesh
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Definition: block_preconditioner.h:2866
oomph::CRDoubleMatrix
Definition: matrices.h:888
oomph::CRDoubleMatrix::multiply
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1782
oomph::MatrixVectorProduct
Definition: matrix_vector_product.h:51
oomph::Preconditioner::setup
void setup(DoubleMatrixBase *matrix_pt)
Definition: preconditioner.h:94
oomph::SuperLUPreconditioner
An interface to allow SuperLU to be used as an (exact) Preconditioner.
Definition: SuperLU_preconditioner.h:40
oomph::TimingHelpers::timer
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319

References oomph::TerminateHelper::clean_up_memory(), i, oomph::CRDoubleMatrix::multiply(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, and oomph::SuperLUPreconditioner::setup().

Referenced by oomph::FSIPreconditioner::setup().

◆ setup() [2/4]

virtual void oomph::Preconditioner::setup
virtual

for some reason we have to remind the compiler that there is a setup() function in Preconditioner base class.

Implements oomph::Preconditioner.

◆ setup() [3/4]

void oomph::Preconditioner::setup
inlinevirtual

for some reason we have to remind the compiler that there is a setup() function in Preconditioner base class.

Implements oomph::Preconditioner.

121  {
122  ObsoleteCode::obsolete();
123  setup(matrix_pt);
124  }
oomph::NavierStokesSchurComplementPreconditioner::setup
void setup()
Setup the preconditioner.
Definition: driven_cavity_with_simple_lsc_preconditioner.cc:167
oomph::ObsoleteCode::obsolete
void obsolete()
Output warning message.
Definition: oomph_utilities.cc:1102

◆ setup() [4/4]

void oomph::Preconditioner::setup
inlinevirtual

for some reason we have to remind the compiler that there is a setup() function in Preconditioner base class.

Implements oomph::Preconditioner.

95  {
96  // Store matrix pointer
97  set_matrix_pt(matrix_pt);
98 
99  // Extract and store communicator pointer
100  DistributableLinearAlgebraObject* dist_obj_pt =
101  dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
102  if (dist_obj_pt != 0)
103  {
104  set_comm_pt(dist_obj_pt->distribution_pt()->communicator_pt());
105  }
106  else
107  {
108  set_comm_pt(0);
109  }
110 
111  double setup_time_start = TimingHelpers::timer();
112  setup();
113  double setup_time_finish = TimingHelpers::timer();
114  Setup_time = setup_time_finish - setup_time_start;
115  }
oomph::DistributableLinearAlgebraObject::DistributableLinearAlgebraObject
DistributableLinearAlgebraObject()
Default constructor - create a distribution.
Definition: linear_algebra_distribution.h:438
oomph::Preconditioner::Setup_time
double Setup_time
The time it takes to set up this preconditioner.
Definition: preconditioner.h:243
oomph::Preconditioner::set_comm_pt
virtual void set_comm_pt(const OomphCommunicator *const comm_pt)
Set the communicator pointer.
Definition: preconditioner.h:193
oomph::Preconditioner::set_matrix_pt
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
Definition: preconditioner.h:165

◆ unpin_first_pressure_dof_in_press_adv_diff()

void oomph::NavierStokesSchurComplementPreconditioner::unpin_first_pressure_dof_in_press_adv_diff ( )
inline

Set boolean indicating that we do not want to pin first pressure dof in Navier Stokes mesh when assembling the pressure advection diffusion system for Fp preconditoner – needed at zero Reynolds number for non-enclosed flows but seems harmless in any case

847  {
848  Pin_first_pressure_dof_in_press_adv_diff = false;
849  }

References Pin_first_pressure_dof_in_press_adv_diff.

◆ use_fp()

void oomph::NavierStokesSchurComplementPreconditioner::use_fp ( )
inline

Use Fp version of the preconditioner.

789  {
790  Use_LSC = false;
791  }

References Use_LSC.

Referenced by FpTestProblem::FpTestProblem().

◆ use_lsc()

void oomph::NavierStokesSchurComplementPreconditioner::use_lsc ( )
inline

Use LSC version of the preconditioner.

783  {
784  Use_LSC = true;
785  }

References Use_LSC.

Referenced by FpTestProblem::FpTestProblem(), and TiltedCavityProblem< ELEMENT >::TiltedCavityProblem().

◆ validate()

template<class ELEMENT >
void oomph::NavierStokesSchurComplementPreconditioner::validate ( DocInfo &  doc_info,
Problem *  orig_problem_pt 
)
inline

Validate auxiliary pressure advection diffusion problem in 2D

855  {
856  FpPressureAdvectionDiffusionProblem<ELEMENT> aux_problem(
857  Navier_stokes_mesh_pt, orig_problem_pt);
858  aux_problem.validate(doc_info);
859  }

References Navier_stokes_mesh_pt, and oomph::FpPressureAdvectionDiffusionProblem< ELEMENT >::validate().

Member Data Documentation

◆ Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements

bool oomph::NavierStokesSchurComplementPreconditioner::Accept_non_NavierStokesElementWithDiagonalMassMatrices_elements
private

Boolean to indicate that presence of elements that are not of type NavierStokesElementWithDiagonalMassMatrices is acceptable (suppresses warning that issued otherwise).

Referenced by disable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements(), enable_accept_non_NavierStokesElementWithDiagonalMassMatrices_elements(), and NavierStokesSchurComplementPreconditioner().

◆ Allow_multiple_element_type_in_navier_stokes_mesh

bool oomph::NavierStokesSchurComplementPreconditioner::Allow_multiple_element_type_in_navier_stokes_mesh
private

Flag to indicate if there are multiple element types in the Navier-Stokes mesh.

Referenced by NavierStokesSchurComplementPreconditioner(), and set_navier_stokes_mesh().

◆ Bt_mat_vec_pt

MatrixVectorProduct* oomph::NavierStokesSchurComplementPreconditioner::Bt_mat_vec_pt
private

MatrixVectorProduct operator for Bt.

Referenced by clean_up_memory(), and NavierStokesSchurComplementPreconditioner().

◆ Doc_time

bool oomph::NavierStokesSchurComplementPreconditioner::Doc_time
private

Set Doc_time to true for outputting results of timings.

Referenced by disable_doc_time(), enable_doc_time(), and NavierStokesSchurComplementPreconditioner().

◆ E_mat_vec_pt

MatrixVectorProduct* oomph::NavierStokesSchurComplementPreconditioner::E_mat_vec_pt
private

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

Referenced by clean_up_memory(), and NavierStokesSchurComplementPreconditioner().

◆ Eqn_number_backup

std::map<Data*, std::vector<int> > oomph::NavierStokesSchurComplementPreconditioner::Eqn_number_backup
private

Map to store original eqn numbers of various Data values when assembling pressure advection diffusion matrix

Referenced by pin_all_non_pressure_dofs(), and reset_pin_status().

◆ F_mat_vec_pt

MatrixVectorProduct* oomph::NavierStokesSchurComplementPreconditioner::F_mat_vec_pt
private

MatrixVectorProduct operator for F.

Referenced by clean_up_memory(), and NavierStokesSchurComplementPreconditioner().

◆ F_preconditioner_is_block_preconditioner

bool oomph::NavierStokesSchurComplementPreconditioner::F_preconditioner_is_block_preconditioner
private

Boolean indicating whether the momentum system preconditioner is a block preconditioner

◆ F_preconditioner_pt

Preconditioner* oomph::NavierStokesSchurComplementPreconditioner::F_preconditioner_pt
private

Pointer to the 'preconditioner' for the F matrix.

Referenced by clean_up_memory(), NavierStokesSchurComplementPreconditioner(), set_f_preconditioner(), and set_f_superlu_preconditioner().

◆ Navier_stokes_mesh_pt

Mesh* oomph::NavierStokesSchurComplementPreconditioner::Navier_stokes_mesh_pt
private

the pointer to the mesh of block preconditionable Navier Stokes elements.

Referenced by assemble_inv_press_and_veloc_mass_matrix_diagonal(), get_pressure_advection_diffusion_matrix(), NavierStokesSchurComplementPreconditioner(), pin_all_non_pressure_dofs(), reset_pin_status(), set_navier_stokes_mesh(), and validate().

◆ P_preconditioner_pt

Preconditioner* oomph::NavierStokesSchurComplementPreconditioner::P_preconditioner_pt
private

Pointer to the 'preconditioner' for the pressure matrix.

Referenced by clean_up_memory(), NavierStokesSchurComplementPreconditioner(), set_p_preconditioner(), and set_p_superlu_preconditioner().

◆ Pin_first_pressure_dof_in_press_adv_diff

bool oomph::NavierStokesSchurComplementPreconditioner::Pin_first_pressure_dof_in_press_adv_diff
private

Boolean indicating if we want to pin first pressure dof in Navier Stokes mesh when assembling the pressure advection diffusion system for Fp preconditoner – needed at zero Reynolds number for non-enclosed flows but seems harmless in any case

Referenced by get_pressure_advection_diffusion_matrix(), NavierStokesSchurComplementPreconditioner(), pin_first_pressure_dof_in_press_adv_diff(), and unpin_first_pressure_dof_in_press_adv_diff().

◆ Preconditioner_has_been_setup

bool oomph::NavierStokesSchurComplementPreconditioner::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(), and NavierStokesSchurComplementPreconditioner().

◆ Problem_pt

Problem* oomph::NavierStokesSchurComplementPreconditioner::Problem_pt
private

Storage for the (non-const!) problem pointer for use in get_pressure_advection_diffusion_matrix().

Referenced by problem_pt(), and set_problem_pt().

◆ QBt_mat_vec_pt

MatrixVectorProduct* oomph::NavierStokesSchurComplementPreconditioner::QBt_mat_vec_pt
private

MatrixVectorProduct operator for Qv^{-1} Bt.

Referenced by clean_up_memory(), and NavierStokesSchurComplementPreconditioner().

◆ Use_LSC

bool oomph::NavierStokesSchurComplementPreconditioner::Use_LSC
private

Boolean to indicate use of LSC (true) or Fp (false) variant.

Referenced by assemble_inv_press_and_veloc_mass_matrix_diagonal(), NavierStokesSchurComplementPreconditioner(), use_fp(), and use_lsc().

◆ Use_robin_for_fp

bool oomph::NavierStokesSchurComplementPreconditioner::Use_robin_for_fp
private

Use Robin BC elements for Fp preconditoner?

Referenced by disable_robin_for_fp(), enable_robin_for_fp(), get_pressure_advection_diffusion_matrix(), and NavierStokesSchurComplementPreconditioner().

◆ Using_default_f_preconditioner

bool oomph::NavierStokesSchurComplementPreconditioner::Using_default_f_preconditioner
private

flag indicating whether the default F preconditioner is used

Referenced by clean_up_memory(), NavierStokesSchurComplementPreconditioner(), set_f_preconditioner(), and set_f_superlu_preconditioner().

◆ Using_default_p_preconditioner

bool oomph::NavierStokesSchurComplementPreconditioner::Using_default_p_preconditioner
private

flag indicating whether the default P preconditioner is used

Referenced by clean_up_memory(), NavierStokesSchurComplementPreconditioner(), set_p_preconditioner(), and set_p_superlu_preconditioner().

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