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

#include <solid_preconditioners.h>

+ Inheritance diagram for oomph::PressureBasedSolidLSCPreconditioner:

Public Member Functions

 PressureBasedSolidLSCPreconditioner ()
 Constructor - sets defaults for control flags. More...
 
 ~PressureBasedSolidLSCPreconditioner ()
 Destructor. More...
 
 PressureBasedSolidLSCPreconditioner (const PressureBasedSolidLSCPreconditioner &)=delete
 Broken copy constructor. More...
 
void setup ()
 Broken assignment operator. More...
 
void preconditioner_solve (const DoubleVector &r, DoubleVector &z)
 Apply preconditioner to Vector r. More...
 
void set_solid_mesh (Mesh *mesh_pt)
 
void enable_p_matrix_scaling ()
 
void disable_p_matrix_scaling ()
 
bool is_p_matrix_using_scaling () const
 
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 set_f_superlu_preconditioner ()
 
void enable_doc_time ()
 Enable documentation of time. More...
 
void disable_doc_time ()
 Disable documentation of time. More...
 
void enable_form_BFBt_product ()
 
void disable_form_BFBt_product ()
 
void clean_up_memory ()
 Helper function to delete preconditioner data. More...
 
- Public Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
 BlockPreconditioner ()
 Constructor. More...
 
 BlockPreconditioner (const BlockPreconditioner &)=delete
 Broken copy constructor. More...
 
virtual ~BlockPreconditioner ()
 Destructor. More...
 
void operator= (const BlockPreconditioner &)=delete
 Broken assignment operator. More...
 
CRDoubleMatrix * matrix_pt () const
 
void turn_on_recursive_debug_flag ()
 
void turn_off_recursive_debug_flag ()
 
void turn_on_debug_flag ()
 Toggles on the debug flag. More...
 
void turn_off_debug_flag ()
 Toggles off the debug flag. More...
 
void turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
 
void turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)
 
virtual void block_setup ()
 
void block_setup (const Vector< unsigned > &dof_to_block_map)
 
void get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
 
CRDoubleMatrix get_block (const unsigned &i, const unsigned &j, const bool &ignore_replacement_block=false) const
 
void set_master_matrix_pt (CRDoubleMatrix *in_matrix_pt)
 Set the matrix_pt in the upper-most master preconditioner. More...
 
void get_block_other_matrix (const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)
 
void get_blocks (DenseMatrix< bool > &required_blocks, DenseMatrix< CRDoubleMatrix * > &block_matrix_pt) const
 
void get_dof_level_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block=false) const
 
void get_dof_level_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block) const
 
CRDoubleMatrix get_concatenated_block (const VectorMatrix< BlockSelector > &selected_block)
 
void get_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)
 
void return_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
 Takes concatenated block ordered vector, b, and copies its. More...
 
void get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
 
void get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const
 
void return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
 
void return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const
 
void get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const
 
void return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const
 
void get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w)
 
void return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const
 
unsigned nblock_types () const
 Return the number of block types. More...
 
unsigned ndof_types () const
 Return the total number of DOF types. More...
 
const Mesh * mesh_pt (const unsigned &i) const
 
unsigned nmesh () const
 
int block_number (const unsigned &i_dof) const
 Return the block number corresponding to a global index i_dof. More...
 
int index_in_block (const unsigned &i_dof) const
 
const LinearAlgebraDistribution * block_distribution_pt (const unsigned &b) const
 Access function to the block distributions (const version). More...
 
LinearAlgebraDistribution * block_distribution_pt (const unsigned b)
 Access function to the block distributions (non-const version). More...
 
LinearAlgebraDistribution * dof_block_distribution_pt (const unsigned &b)
 Access function to the dof-level block distributions. More...
 
const LinearAlgebraDistribution * master_distribution_pt () const
 
unsigned ndof_types_in_mesh (const unsigned &i) const
 
bool is_subsidiary_block_preconditioner () const
 
bool is_master_block_preconditioner () const
 
void set_block_output_to_files (const std::string &basefilename)
 
void disable_block_output_to_files ()
 Turn off output of blocks (by clearing the basefilename string). More...
 
bool block_output_on () const
 Test if output of blocks is on or not. More...
 
void output_blocks_to_files (const std::string &basefilename, const unsigned &precision=8) const
 
void post_block_matrix_assembly_partial_clear ()
 
BlockPreconditioner< CRDoubleMatrix > * master_block_preconditioner_pt () const
 Access function to the master block preconditioner pt. More...
 
void clear_block_preconditioner_base ()
 
void document ()
 
Vector< Vector< unsigned > > doftype_coarsen_map_fine () const
 
Vector< unsigned > get_fine_grain_dof_types_in (const unsigned &i) const
 
unsigned nfine_grain_dof_types_in (const unsigned &i) const
 
MapMatrix< unsigned, CRDoubleMatrix * > replacement_dof_block_pt () const
 Access function to the replaced dof-level blocks. More...
 
void setup_matrix_vector_product (MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
 
void setup_matrix_vector_product (MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const unsigned &block_col_index)
 
void internal_get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w) const
 
void internal_return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const
 
unsigned internal_nblock_types () const
 
unsigned internal_ndof_types () const
 
void internal_return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const
 
void internal_get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const
 
void internal_get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
 
void internal_get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const
 
void internal_return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
 
void internal_return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const
 
void internal_get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block) const
 
void internal_get_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block) const
 
int internal_block_number (const unsigned &i_dof) const
 
int internal_index_in_block (const unsigned &i_dof) const
 
const LinearAlgebraDistribution * internal_block_distribution_pt (const unsigned &b) const
 Access function to the internal block distributions. More...
 
void insert_auxiliary_block_distribution (const Vector< unsigned > &block_vec_number, LinearAlgebraDistribution *dist_pt)
 
void block_matrix_test (const unsigned &i, const unsigned &j, const CRDoubleMatrix *block_matrix_pt) const
 
int get_index_of_value (const Vector< myType > &vec, const myType val, const bool sorted=false) const
 
- Public Member Functions inherited from oomph::Preconditioner
 Preconditioner ()
 Constructor. More...
 
 Preconditioner (const Preconditioner &)=delete
 Broken copy constructor. More...
 
void operator= (const Preconditioner &)=delete
 Broken assignment operator. More...
 
virtual ~Preconditioner ()
 Destructor (empty) More...
 
virtual void preconditioner_solve_transpose (const DoubleVector &r, DoubleVector &z)
 
void setup (DoubleMatrixBase *matrix_pt)
 
void setup (const Problem *problem_pt, DoubleMatrixBase *matrix_pt)
 
void enable_silent_preconditioner_setup ()
 Set up the block preconditioner quietly! More...
 
void disable_silent_preconditioner_setup ()
 Be verbose in the block preconditioner setup. More...
 
virtual void set_matrix_pt (DoubleMatrixBase *matrix_pt)
 Set the matrix pointer. More...
 
virtual const OomphCommunicator * comm_pt () const
 Get function for comm pointer. More...
 
virtual void set_comm_pt (const OomphCommunicator *const comm_pt)
 Set the communicator pointer. More...
 
double setup_time () const
 Returns the time to setup the preconditioner. More...
 
virtual void turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
 
virtual void turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)
 
- Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
 DistributableLinearAlgebraObject ()
 Default constructor - create a distribution. More...
 
 DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
 Broken copy constructor. More...
 
void operator= (const DistributableLinearAlgebraObject &)=delete
 Broken assignment operator. More...
 
virtual ~DistributableLinearAlgebraObject ()
 Destructor. More...
 
LinearAlgebraDistribution * distribution_pt () const
 access to the LinearAlgebraDistribution More...
 
unsigned nrow () const
 access function to the number of global rows. More...
 
unsigned nrow_local () const
 access function for the num of local rows on this processor. More...
 
unsigned nrow_local (const unsigned &p) const
 access function for the num of local rows on this processor. More...
 
unsigned first_row () const
 access function for the first row on this processor More...
 
unsigned first_row (const unsigned &p) const
 access function for the first row on this processor More...
 
bool distributed () const
 distribution is serial or distributed More...
 
bool distribution_built () const
 
void build_distribution (const LinearAlgebraDistribution *const dist_pt)
 
void build_distribution (const LinearAlgebraDistribution &dist)
 

Private Member Functions

CRDoubleMatrix * assemble_mass_matrix_diagonal ()
 

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 P_matrix_using_scaling
 
bool F_preconditioner_is_block_preconditioner
 
bool Doc_time
 Set Doc_time to true for outputting results of timings. More...
 
MatrixVectorProduct * F_mat_vec_pt
 MatrixVectorProduct operator for F if BFBt is not to be formed. More...
 
MatrixVectorProduct * QBt_mat_vec_pt
 MatrixVectorProduct operator for QBt if BFBt is not to be formed. More...
 
MatrixVectorProduct * Bt_mat_vec_pt
 MatrixVectorProduct operator for Bt;. More...
 
MatrixVectorProduct * E_mat_vec_pt
 MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed. More...
 
bool Form_BFBt_product
 
Mesh * Solid_mesh_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) preconditioner. It uses blocks corresponding to the displacement/position and pressure unknowns, i.e. there are a total of 2x2 blocks, and all displacement/position components are treated as a single block of unknowns.

Here are the details: An "ideal" preconditioner would solve the saddle point 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}\), \( {\bf G} \), and \( {\bf D}\) are the blocks that arise in the Jacobian of the pressure-based equations of linear and nonlinear elasticity (with dofs in order of displacement/position and pressure). 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 displacement/position 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

PressureBasedSolidLSCPreconditioner::set_f_preconditioner(...)
oomph::PressureBasedSolidLSCPreconditioner::set_f_preconditioner
void set_f_preconditioner(Preconditioner *new_f_preconditioner_pt)
Function to set a new momentum matrix preconditioner (inexact solver)
Definition: solid_preconditioners.h:233

or

PressureBasedSolidLSCPreconditioner::set_p_preconditioner(...)
oomph::PressureBasedSolidLSCPreconditioner::set_p_preconditioner
void set_p_preconditioner(Preconditioner *new_p_preconditioner_pt)
Function to set a new pressure matrix preconditioner (inexact solver)
Definition: solid_preconditioners.h:209

Constructor & Destructor Documentation

◆ PressureBasedSolidLSCPreconditioner() [1/2]

oomph::PressureBasedSolidLSCPreconditioner::PressureBasedSolidLSCPreconditioner ( )
inline

Constructor - sets defaults for control flags.

111  : BlockPreconditioner<CRDoubleMatrix>()
112  {
113  // Flag to indicate that the preconditioner has been setup
114  // previously -- if setup() is called again, data can
115  // be wiped.
116  Preconditioner_has_been_setup = false;
117 
118  // By default we use SuperLU for both p and f blocks
119  Using_default_p_preconditioner = true;
120  Using_default_f_preconditioner = true;
121 
122  // resize the mesh pt
123  // note: meaningless if subsidiary preconditioner
124  this->set_nmesh(1);
125  Solid_mesh_pt = 0;
126 
127  // Set default preconditioners (inexact solvers) -- they are
128  // members of this class!
129  P_preconditioner_pt = 0;
130  F_preconditioner_pt = 0;
131 
132  // Flag to determine if mass matrix diagonal Q^{-1}
133  // is used for scaling.
134  P_matrix_using_scaling = true;
135 
136  // set Doc_time to false
137  Doc_time = false;
138 
139  // null the off diagonal Block matrix pt
140  Bt_mat_vec_pt = 0;
141 
142  // null the F matrix vector product helper
143  F_mat_vec_pt = 0;
144 
145  // null the QBt matrix vector product pt
146  QBt_mat_vec_pt = 0;
147 
148  // null the E matrix vector product helper
149  E_mat_vec_pt = 0;
150 
151  // by default we do not form the E matrix (BQFQBt)
152  Form_BFBt_product = false;
153  }
oomph::BlockPreconditioner< CRDoubleMatrix >::set_nmesh
void set_nmesh(const unsigned &n)
Definition: block_preconditioner.h:2851
oomph::PressureBasedSolidLSCPreconditioner::QBt_mat_vec_pt
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for QBt if BFBt is not to be formed.
Definition: solid_preconditioners.h:336
oomph::PressureBasedSolidLSCPreconditioner::Form_BFBt_product
bool Form_BFBt_product
Definition: solid_preconditioners.h:356
oomph::PressureBasedSolidLSCPreconditioner::Preconditioner_has_been_setup
bool Preconditioner_has_been_setup
Definition: solid_preconditioners.h:314
oomph::PressureBasedSolidLSCPreconditioner::F_preconditioner_pt
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
Definition: solid_preconditioners.h:303
oomph::PressureBasedSolidLSCPreconditioner::E_mat_vec_pt
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.
Definition: solid_preconditioners.h:342
oomph::PressureBasedSolidLSCPreconditioner::Doc_time
bool Doc_time
Set Doc_time to true for outputting results of timings.
Definition: solid_preconditioners.h:330
oomph::PressureBasedSolidLSCPreconditioner::Solid_mesh_pt
Mesh * Solid_mesh_pt
Definition: solid_preconditioners.h:360
oomph::PressureBasedSolidLSCPreconditioner::Bt_mat_vec_pt
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt;.
Definition: solid_preconditioners.h:339
oomph::PressureBasedSolidLSCPreconditioner::Using_default_p_preconditioner
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
Definition: solid_preconditioners.h:309
oomph::PressureBasedSolidLSCPreconditioner::F_mat_vec_pt
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F if BFBt is not to be formed.
Definition: solid_preconditioners.h:333
oomph::PressureBasedSolidLSCPreconditioner::Using_default_f_preconditioner
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
Definition: solid_preconditioners.h:306
oomph::PressureBasedSolidLSCPreconditioner::P_preconditioner_pt
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
Definition: solid_preconditioners.h:300
oomph::PressureBasedSolidLSCPreconditioner::P_matrix_using_scaling
bool P_matrix_using_scaling
Definition: solid_preconditioners.h:318

References Bt_mat_vec_pt, Doc_time, E_mat_vec_pt, F_mat_vec_pt, F_preconditioner_pt, Form_BFBt_product, P_matrix_using_scaling, P_preconditioner_pt, Preconditioner_has_been_setup, QBt_mat_vec_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::set_nmesh(), Solid_mesh_pt, Using_default_f_preconditioner, and Using_default_p_preconditioner.

◆ ~PressureBasedSolidLSCPreconditioner()

oomph::PressureBasedSolidLSCPreconditioner::~PressureBasedSolidLSCPreconditioner ( )
inline

Destructor.

157  {
158  clean_up_memory();
159  }
oomph::PressureBasedSolidLSCPreconditioner::clean_up_memory
void clean_up_memory()
Helper function to delete preconditioner data.
Definition: solid_preconditioners.cc:1357

References clean_up_memory().

◆ PressureBasedSolidLSCPreconditioner() [2/2]

oomph::PressureBasedSolidLSCPreconditioner::PressureBasedSolidLSCPreconditioner ( const PressureBasedSolidLSCPreconditioner &  )
delete

Broken copy constructor.

Member Function Documentation

◆ assemble_mass_matrix_diagonal()

CRDoubleMatrix * oomph::PressureBasedSolidLSCPreconditioner::assemble_mass_matrix_diagonal ( )
private

Helper function to assemble the diagonal of the mass matrix from the elemental contributions defined in PressureBasedSolidEquations<DIM>::get_mass_matrix_diagonal(...).

Helper function to assemble the diagonal of the mass matrix from the elemental contributions defined in SolidElementWithDiagonalMassMatrix::get_mass_matrix_diagonal(...).

706  {
707  // determine the rows required by this processor
708  unsigned first_row = this->block_distribution_pt(0)->first_row();
709  unsigned nrow_local = this->block_distribution_pt(0)->nrow_local();
710  unsigned nrow = this->block_distribution_pt(0)->nrow();
711 
712  // create storage for the diagonals
713  double* m_values = new double[nrow_local];
714  for (unsigned i = 0; i < nrow_local; i++)
715  {
716  m_values[i] = 0;
717  }
718 
719  // if the problem is distributed
720  bool distributed = false;
721 #ifdef OOMPH_HAS_MPI
722  if (any_mesh_distributed() || this->master_distribution_pt()->distributed())
723  {
724  distributed = true;
725  }
726 #endif
727 
728  // next we get the diagonal mass matrix data
729  if (distributed)
730  {
731 #ifdef OOMPH_HAS_MPI
732  // the number of processors
733  unsigned nproc = this->comm_pt()->nproc();
734 
735  // and my rank
736  unsigned my_rank = this->comm_pt()->my_rank();
737 
738  // determine the rows for which we have lookup rows
739  // if the problem is NOT distributed then we only classify global equation
740  // on this processor to avoid duplication (as every processor holds
741  // every element)
742  unsigned first_lookup_row = 0;
743  unsigned last_lookup_row = 0;
744  first_lookup_row = this->master_distribution_pt()->first_row();
745  last_lookup_row =
746  first_lookup_row + this->master_distribution_pt()->nrow_local() - 1;
747 
748  // find number of local elements
749  unsigned n_el = Solid_mesh_pt->nelement();
750 
751  // the diagonal mass matrix contributions that have been
752  // classified and should be sent to another processor
753  Vector<double>* classified_contributions_send = new Vector<double>[nproc];
754 
755  // the corresponding block indices
756  Vector<unsigned>* classified_indices_send = new Vector<unsigned>[nproc];
757 
758  // the maitrix contributions that cannot be classified by this processor
759  // and therefore must be sent to another for classification
760  Vector<double>* unclassified_contributions_send =
761  new Vector<double>[nproc];
762 
763  // the corresponding global indices that require classification
764  Vector<unsigned>* unclassified_indices_send = new Vector<unsigned>[nproc];
765 
766  // get the master distribution pt
767  const LinearAlgebraDistribution* master_distribution_pt =
768  this->master_distribution_pt();
769 
770  // get the displacement/position distribution pt
771  const LinearAlgebraDistribution* displ_dist_pt =
772  this->block_distribution_pt(0);
773 
774  // get the contribution for each element
775  for (unsigned e = 0; e < n_el; e++)
776  {
777  // check that the element is not halo d
778  if (!Solid_mesh_pt->element_pt(e)->is_halo())
779  {
780  // find number of degrees of freedom in the element
781  // (this is slightly too big because it includes the
782  // pressure dofs but this doesn't matter)
783  unsigned el_dof = Solid_mesh_pt->element_pt(e)->ndof();
784 
785  // allocate local storage for the element's contribution to the
786  // mass matrix diagonal
787  Vector<double> el_vmm_diagonal(el_dof);
788 
789  SolidElementWithDiagonalMassMatrix* cast_el_pt =
790  dynamic_cast<SolidElementWithDiagonalMassMatrix*>(
791  Solid_mesh_pt->element_pt(e));
792 
793 
794  if (cast_el_pt == 0)
795  {
796 #ifdef PARANOID
797  std::ostringstream error_message;
798  error_message << "Failed cast to "
799  << "SolidElementWithDiagonalMassMatrix*\n"
800  << "Element is of type: "
801  << typeid(*(Solid_mesh_pt->element_pt(e))).name()
802  << "\n"
803  << typeid(Solid_mesh_pt->element_pt(e)).name()
804  << std::endl;
805  OomphLibWarning(error_message.str(),
806  "PressureBasedSolidLSCPreconditioner::assemble_"
807  "mass_matrix_diagonal()",
808  OOMPH_EXCEPTION_LOCATION);
809 #endif
810  }
811  else
812  {
813  cast_el_pt->get_mass_matrix_diagonal(el_vmm_diagonal);
814  }
815 
816  // get the contribution for each dof
817  for (unsigned i = 0; i < el_dof; i++)
818  {
819  // Get the equation number
820  unsigned eqn_number = Solid_mesh_pt->element_pt(e)->eqn_number(i);
821 
822  // if I have lookup information on this processor
823  if (eqn_number >= first_lookup_row && eqn_number <= last_lookup_row)
824  {
825  // bypass non displacement/position DOFs
826  if (this->block_number(eqn_number) == 0)
827  {
828  // get the index in the block
829  unsigned index = this->index_in_block(eqn_number);
830 
831  // determine which processor requires the block index
832  for (unsigned p = 0; p < nproc; p++)
833  {
834  if (index >= displ_dist_pt->first_row(p) &&
835  (index < (displ_dist_pt->first_row(p) +
836  displ_dist_pt->nrow_local(p))))
837  {
838  // if it is required by this processor then add the
839  // contribution
840  if (p == my_rank)
841  {
842  m_values[index - first_row] += el_vmm_diagonal[i];
843  }
844  // other wise store it for communication
845  else
846  {
847  classified_contributions_send[p].push_back(
848  el_vmm_diagonal[i]);
849  classified_indices_send[p].push_back(index);
850  }
851  }
852  }
853  }
854  }
855 
856  // if we do not have the lookup information on this processor
857  // then we send the mass matrix contribution to a processor
858  // which we know does have the lookup information
859  // the assumption: the processor for which the master block
860  // preconditioner distribution will definitely hold the lookup
861  // data for eqn_number (although others may)
862  else if (any_mesh_distributed())
863  {
864  // determine which processor requires the block index
865  unsigned p = 0;
866  while (!(eqn_number >= master_distribution_pt->first_row(p) &&
867  (eqn_number < (master_distribution_pt->first_row(p) +
868  master_distribution_pt->nrow_local(p)))))
869  {
870  p++;
871  }
872 
873  // store the data
874  unclassified_contributions_send[p].push_back(el_vmm_diagonal[i]);
875  unclassified_indices_send[p].push_back(eqn_number);
876  }
877  }
878  }
879  }
880 
881  // next the unclassified contributions are communicated to
882  // processors that can classify them
883 
884  // first determine how many unclassified rows are to be sent to
885  // each processor
886  unsigned* n_unclassified_send = new unsigned[nproc];
887  for (unsigned p = 0; p < nproc; p++)
888  {
889  if (p == my_rank)
890  {
891  n_unclassified_send[p] = 0;
892  }
893  else
894  {
895  n_unclassified_send[p] = unclassified_contributions_send[p].size();
896  }
897  }
898 
899  // then all-to-all number of unclassified to be sent / recv
900  unsigned* n_unclassified_recv = new unsigned[nproc];
901  MPI_Alltoall(n_unclassified_send,
902  1,
903  MPI_UNSIGNED,
904  n_unclassified_recv,
905  1,
906  MPI_UNSIGNED,
907  this->comm_pt()->mpi_comm());
908 
909  // the base displacement for the sends
910  MPI_Aint base_displacement;
911  MPI_Get_address(m_values, &base_displacement);
912 
913  // allocate storage for the data to be recieved
914  // and post the sends and recvs
915  Vector<double*> unclassified_contributions_recv(nproc);
916  Vector<unsigned*> unclassified_indices_recv(nproc);
917  Vector<MPI_Request> unclassified_recv_requests;
918  Vector<MPI_Request> unclassified_send_requests;
919  Vector<unsigned> unclassified_recv_proc;
920  for (unsigned p = 0; p < nproc; p++)
921  {
922  if (p != my_rank)
923  {
924  // recv
925  if (n_unclassified_recv[p] > 0)
926  {
927  unclassified_contributions_recv[p] =
928  new double[n_unclassified_recv[p]];
929  unclassified_indices_recv[p] = new unsigned[n_unclassified_recv[p]];
930 
931  // data for the struct data type
932  MPI_Datatype recv_types[2];
933  MPI_Aint recv_displacements[2];
934  int recv_sz[2];
935 
936  // contributions
937  MPI_Type_contiguous(
938  n_unclassified_recv[p], MPI_DOUBLE, &recv_types[0]);
939  MPI_Type_commit(&recv_types[0]);
940  MPI_Get_address(unclassified_contributions_recv[p],
941  &recv_displacements[0]);
942  recv_displacements[0] -= base_displacement;
943  recv_sz[0] = 1;
944 
945  // indices
946  MPI_Type_contiguous(
947  n_unclassified_recv[p], MPI_UNSIGNED, &recv_types[1]);
948  MPI_Type_commit(&recv_types[1]);
949  MPI_Get_address(unclassified_indices_recv[p],
950  &recv_displacements[1]);
951  recv_displacements[1] -= base_displacement;
952  recv_sz[1] = 1;
953 
954  // build the final recv type
955  MPI_Datatype final_recv_type;
956  MPI_Type_create_struct(
957  2, recv_sz, recv_displacements, recv_types, &final_recv_type);
958  MPI_Type_commit(&final_recv_type);
959 
960  // and recv
961  MPI_Request req;
962  MPI_Irecv(
963  m_values, 1, final_recv_type, p, 0, comm_pt()->mpi_comm(), &req);
964  unclassified_recv_requests.push_back(req);
965  unclassified_recv_proc.push_back(p);
966  MPI_Type_free(&recv_types[0]);
967  MPI_Type_free(&recv_types[1]);
968  MPI_Type_free(&final_recv_type);
969  }
970 
971  // send
972  if (n_unclassified_send[p] > 0)
973  {
974  // data for the struct data type
975  MPI_Datatype send_types[2];
976  MPI_Aint send_displacements[2];
977  int send_sz[2];
978 
979  // contributions
980  MPI_Type_contiguous(
981  n_unclassified_send[p], MPI_DOUBLE, &send_types[0]);
982  MPI_Type_commit(&send_types[0]);
983  MPI_Get_address(&unclassified_contributions_send[p][0],
984  &send_displacements[0]);
985  send_displacements[0] -= base_displacement;
986  send_sz[0] = 1;
987 
988  // indices
989  MPI_Type_contiguous(
990  n_unclassified_send[p], MPI_UNSIGNED, &send_types[1]);
991  MPI_Type_commit(&send_types[1]);
992  MPI_Get_address(&unclassified_indices_send[p][0],
993  &send_displacements[1]);
994  send_displacements[1] -= base_displacement;
995  send_sz[1] = 1;
996 
997  // build the final send type
998  MPI_Datatype final_send_type;
999  MPI_Type_create_struct(
1000  2, send_sz, send_displacements, send_types, &final_send_type);
1001  MPI_Type_commit(&final_send_type);
1002 
1003  // and send
1004  MPI_Request req;
1005  MPI_Isend(
1006  m_values, 1, final_send_type, p, 0, comm_pt()->mpi_comm(), &req);
1007  unclassified_send_requests.push_back(req);
1008  MPI_Type_free(&send_types[0]);
1009  MPI_Type_free(&send_types[1]);
1010  MPI_Type_free(&final_send_type);
1011  }
1012  }
1013  }
1014 
1015  // next classify the data as it is received
1016  unsigned n_unclassified_recv_req = unclassified_recv_requests.size();
1017  while (n_unclassified_recv_req > 0)
1018  {
1019  // get the processor number and remove the completed request
1020  // for the vector of requests
1021  int req_num;
1022  MPI_Waitany(n_unclassified_recv_req,
1023  &unclassified_recv_requests[0],
1024  &req_num,
1025  MPI_STATUS_IGNORE);
1026  unsigned p = unclassified_recv_proc[req_num];
1027  unclassified_recv_requests.erase(unclassified_recv_requests.begin() +
1028  req_num);
1029  unclassified_recv_proc.erase(unclassified_recv_proc.begin() + req_num);
1030  n_unclassified_recv_req--;
1031 
1032  // next classify the dofs
1033  // and store them for sending to other processors if required
1034  unsigned n_recv = n_unclassified_recv[p];
1035  for (unsigned i = 0; i < n_recv; i++)
1036  {
1037  unsigned eqn_number = unclassified_indices_recv[p][i];
1038  // bypass non displacement/position DOFs
1039  if (this->block_number(eqn_number) == 0)
1040  {
1041  // get the index in the block
1042  unsigned index = this->index_in_block(eqn_number);
1043 
1044  // determine which processor requires the block index
1045  for (unsigned pp = 0; pp < nproc; pp++)
1046  {
1047  if (index >= displ_dist_pt->first_row(pp) &&
1048  (index < (displ_dist_pt->first_row(pp) +
1049  displ_dist_pt->nrow_local(pp))))
1050  {
1051  // if it is required by this processor then add the
1052  // contribution
1053  if (pp == my_rank)
1054  {
1055  m_values[index - first_row] +=
1056  unclassified_contributions_recv[p][i];
1057  }
1058  // other wise store it for communication
1059  else
1060  {
1061  double v = unclassified_contributions_recv[p][i];
1062  classified_contributions_send[pp].push_back(v);
1063  classified_indices_send[pp].push_back(index);
1064  }
1065  }
1066  }
1067  }
1068  }
1069 
1070  // clean up
1071  delete[] unclassified_contributions_recv[p];
1072  delete[] unclassified_indices_recv[p];
1073  }
1074  delete[] n_unclassified_recv;
1075 
1076  // now all indices have been classified
1077 
1078  // next the classified contributions are communicated to
1079  // processors that require them
1080 
1081  // first determine how many classified rows are to be sent to
1082  // each processor
1083  unsigned* n_classified_send = new unsigned[nproc];
1084  for (unsigned p = 0; p < nproc; p++)
1085  {
1086  if (p == my_rank)
1087  {
1088  n_classified_send[p] = 0;
1089  }
1090  else
1091  {
1092  n_classified_send[p] = classified_contributions_send[p].size();
1093  }
1094  }
1095 
1096  // then all-to-all com number of classified to be sent / recv
1097  unsigned* n_classified_recv = new unsigned[nproc];
1098  MPI_Alltoall(n_classified_send,
1099  1,
1100  MPI_UNSIGNED,
1101  n_classified_recv,
1102  1,
1103  MPI_UNSIGNED,
1104  this->comm_pt()->mpi_comm());
1105 
1106  // allocate storage for the data to be recieved
1107  // and post the sends and recvs
1108  Vector<double*> classified_contributions_recv(nproc);
1109  Vector<unsigned*> classified_indices_recv(nproc);
1110  Vector<MPI_Request> classified_recv_requests;
1111  Vector<MPI_Request> classified_send_requests;
1112  Vector<unsigned> classified_recv_proc;
1113  for (unsigned p = 0; p < nproc; p++)
1114  {
1115  if (p != my_rank)
1116  {
1117  // recv
1118  if (n_classified_recv[p] > 0)
1119  {
1120  classified_contributions_recv[p] = new double[n_classified_recv[p]];
1121  classified_indices_recv[p] = new unsigned[n_classified_recv[p]];
1122 
1123  // data for the struct data type
1124  MPI_Datatype recv_types[2];
1125  MPI_Aint recv_displacements[2];
1126  int recv_sz[2];
1127 
1128  // contributions
1129  MPI_Type_contiguous(
1130  n_classified_recv[p], MPI_DOUBLE, &recv_types[0]);
1131  MPI_Type_commit(&recv_types[0]);
1132  MPI_Get_address(classified_contributions_recv[p],
1133  &recv_displacements[0]);
1134  recv_displacements[0] -= base_displacement;
1135  recv_sz[0] = 1;
1136 
1137  // indices
1138  MPI_Type_contiguous(
1139  n_classified_recv[p], MPI_UNSIGNED, &recv_types[1]);
1140  MPI_Type_commit(&recv_types[1]);
1141  MPI_Get_address(classified_indices_recv[p], &recv_displacements[1]);
1142  recv_displacements[1] -= base_displacement;
1143  recv_sz[1] = 1;
1144 
1145  // build the final recv type
1146  MPI_Datatype final_recv_type;
1147  MPI_Type_create_struct(
1148  2, recv_sz, recv_displacements, recv_types, &final_recv_type);
1149  MPI_Type_commit(&final_recv_type);
1150 
1151  // and recv
1152  MPI_Request req;
1153  MPI_Irecv(
1154  m_values, 1, final_recv_type, p, 0, comm_pt()->mpi_comm(), &req);
1155  classified_recv_requests.push_back(req);
1156  classified_recv_proc.push_back(p);
1157  MPI_Type_free(&recv_types[0]);
1158  MPI_Type_free(&recv_types[1]);
1159  MPI_Type_free(&final_recv_type);
1160  }
1161 
1162  // send
1163  if (n_classified_send[p] > 0)
1164  {
1165  // data for the struct data type
1166  MPI_Datatype send_types[2];
1167  MPI_Aint send_displacements[2];
1168  int send_sz[2];
1169 
1170  // contributions
1171  MPI_Type_contiguous(
1172  n_classified_send[p], MPI_DOUBLE, &send_types[0]);
1173  MPI_Type_commit(&send_types[0]);
1174  MPI_Get_address(&classified_contributions_send[p][0],
1175  &send_displacements[0]);
1176  send_displacements[0] -= base_displacement;
1177  send_sz[0] = 1;
1178 
1179  // indices
1180  MPI_Type_contiguous(
1181  n_classified_send[p], MPI_UNSIGNED, &send_types[1]);
1182  MPI_Type_commit(&send_types[1]);
1183  MPI_Get_address(&classified_indices_send[p][0],
1184  &send_displacements[1]);
1185  send_displacements[1] -= base_displacement;
1186  send_sz[1] = 1;
1187 
1188  // build the final send type
1189  MPI_Datatype final_send_type;
1190  MPI_Type_create_struct(
1191  2, send_sz, send_displacements, send_types, &final_send_type);
1192  MPI_Type_commit(&final_send_type);
1193 
1194  // and send
1195  MPI_Request req;
1196  MPI_Isend(
1197  m_values, 1, final_send_type, p, 0, comm_pt()->mpi_comm(), &req);
1198  classified_send_requests.push_back(req);
1199  MPI_Type_free(&send_types[0]);
1200  MPI_Type_free(&send_types[1]);
1201  MPI_Type_free(&final_send_type);
1202  }
1203  }
1204  }
1205 
1206  // next classify the data as it is received
1207  unsigned n_classified_recv_req = classified_recv_requests.size();
1208  while (n_classified_recv_req > 0)
1209  {
1210  // get the processor number and remove the completed request
1211  // for the vector of requests
1212  int req_num;
1213  MPI_Waitany(n_classified_recv_req,
1214  &classified_recv_requests[0],
1215  &req_num,
1216  MPI_STATUS_IGNORE);
1217  unsigned p = classified_recv_proc[req_num];
1218  classified_recv_requests.erase(classified_recv_requests.begin() +
1219  req_num);
1220  classified_recv_proc.erase(classified_recv_proc.begin() + req_num);
1221  n_classified_recv_req--;
1222 
1223  // next classify the dofs
1224  // and store them for sending to other processors if required
1225  unsigned n_recv = n_classified_recv[p];
1226  for (unsigned i = 0; i < n_recv; i++)
1227  {
1228  m_values[classified_indices_recv[p][i] - first_row] +=
1229  classified_contributions_recv[p][i];
1230  }
1231 
1232  // clean up
1233  delete[] classified_contributions_recv[p];
1234  delete[] classified_indices_recv[p];
1235  }
1236 
1237  // wait for the unclassified sends to complete
1238  unsigned n_unclassified_send_req = unclassified_send_requests.size();
1239  if (n_unclassified_send_req > 0)
1240  {
1241  MPI_Waitall(n_unclassified_send_req,
1242  &unclassified_send_requests[0],
1243  MPI_STATUS_IGNORE);
1244  }
1245  delete[] unclassified_contributions_send;
1246  delete[] unclassified_indices_send;
1247  delete[] n_unclassified_send;
1248 
1249  // wait for the classified sends to complete
1250  unsigned n_classified_send_req = classified_send_requests.size();
1251  if (n_classified_send_req > 0)
1252  {
1253  MPI_Waitall(n_classified_send_req,
1254  &classified_send_requests[0],
1255  MPI_STATUS_IGNORE);
1256  }
1257  delete[] classified_indices_send;
1258  delete[] classified_contributions_send;
1259  delete[] n_classified_recv;
1260  delete[] n_classified_send;
1261 #endif
1262  }
1263 
1264  // or if the problem is not distributed
1265  else
1266  {
1267  // find number of elements
1268  unsigned n_el = Solid_mesh_pt->nelement();
1269 
1270  // get the contribution for each element
1271  for (unsigned e = 0; e < n_el; e++)
1272  {
1273  // find number of degrees of freedom in the element
1274  // (this is slightly too big because it includes the
1275  // pressure dofs but this doesn't matter)
1276  unsigned el_dof = Solid_mesh_pt->element_pt(e)->ndof();
1277 
1278  // allocate local storage for the element's contribution to the
1279  // mass matrix diagonal
1280  Vector<double> el_vmm_diagonal(el_dof);
1281 
1282  SolidElementWithDiagonalMassMatrix* cast_el_pt =
1283  dynamic_cast<SolidElementWithDiagonalMassMatrix*>(
1284  Solid_mesh_pt->element_pt(e));
1285 
1286  if (cast_el_pt == 0)
1287  {
1288 #ifdef PARANOID
1289  // #pragma clang diagnostic push
1290  // #pragma clang diagnostic ignored
1291  // "-Wpotentially-evaluated-expression"
1292  std::ostringstream error_message;
1293  error_message << "Failed cast to "
1294  << "SolidElementWithDiagonalMassMatrix*\n"
1295  << "Element is of type: "
1296  << typeid(*(Solid_mesh_pt->element_pt(e))).name()
1297  << "\n"
1298  << typeid(Solid_mesh_pt->element_pt(e)).name()
1299  << std::endl;
1300  OomphLibWarning(error_message.str(),
1301  "PressureBasedSolidLSCPreconditioner::assemble_mass_"
1302  "matrix_diagonal()",
1303  OOMPH_EXCEPTION_LOCATION);
1304 //#pragma clang diagnostic pop
1305 #endif
1306  }
1307  else
1308  {
1309  cast_el_pt->get_mass_matrix_diagonal(el_vmm_diagonal);
1310  }
1311 
1312  // get the contribution for each dof
1313  for (unsigned i = 0; i < el_dof; i++)
1314  {
1315  // Get the equation number
1316  unsigned eqn_number = Solid_mesh_pt->element_pt(e)->eqn_number(i);
1317 
1318  // bypass non displacement/position DOFs
1319  if (this->block_number(eqn_number) == 0)
1320  {
1321  // get the index in the block
1322  unsigned index = this->index_in_block(eqn_number);
1323 
1324  // if it is required on this processor
1325  if (index >= first_row && index < first_row + nrow_local)
1326  {
1327  m_values[index - first_row] += el_vmm_diagonal[i];
1328  }
1329  }
1330  }
1331  }
1332  }
1333 
1334  // create column index and row start
1335  int* m_column_index = new int[nrow_local];
1336  int* m_row_start = new int[nrow_local + 1];
1337  for (unsigned i = 0; i < nrow_local; i++)
1338  {
1339  m_values[i] = 1 / m_values[i];
1340  m_column_index[i] = first_row + i;
1341  m_row_start[i] = i;
1342  }
1343  m_row_start[nrow_local] = nrow_local;
1344 
1345  // build the matrix
1346  CRDoubleMatrix* m_pt = new CRDoubleMatrix(this->block_distribution_pt(0));
1347  m_pt->build_without_copy(
1348  nrow, nrow_local, m_values, m_column_index, m_row_start);
1349 
1350  // return the matrix;
1351  return m_pt;
1352  }
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 >::any_mesh_distributed
bool any_mesh_distributed() const
Definition: block_preconditioner.h:3025
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::DistributableLinearAlgebraObject::nrow
unsigned nrow() const
access function to the number of global rows.
Definition: linear_algebra_distribution.h:463
oomph::DistributableLinearAlgebraObject::nrow_local
unsigned nrow_local() const
access function for the num of local rows on this processor.
Definition: linear_algebra_distribution.h:469
oomph::DistributableLinearAlgebraObject::first_row
unsigned first_row() const
access function for the first row on this processor
Definition: linear_algebra_distribution.h:481
oomph::GeneralisedElement::ndof
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
oomph::GeneralisedElement::eqn_number
unsigned long eqn_number(const unsigned &ieqn_local) const
Definition: elements.h:704
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
plotDoE.name
string name
Definition: plotDoE.py:33
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61

References oomph::BlockPreconditioner< CRDoubleMatrix >::any_mesh_distributed(), 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::DistributableLinearAlgebraObject::first_row(), oomph::SolidElementWithDiagonalMassMatrix::get_mass_matrix_diagonal(), i, oomph::BlockPreconditioner< CRDoubleMatrix >::index_in_block(), oomph::BlockPreconditioner< CRDoubleMatrix >::master_distribution_pt(), oomph::OomphCommunicator::my_rank(), plotDoE::name, oomph::GeneralisedElement::ndof(), oomph::Mesh::nelement(), oomph::OomphCommunicator::nproc(), oomph::LinearAlgebraDistribution::nrow(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::LinearAlgebraDistribution::nrow_local(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_EXCEPTION_LOCATION, p, Solid_mesh_pt, and v.

Referenced by setup().

◆ clean_up_memory()

void oomph::PressureBasedSolidLSCPreconditioner::clean_up_memory ( )
virtual

Helper function to delete preconditioner data.

Reimplemented from oomph::Preconditioner.

1358  {
1359  if (Preconditioner_has_been_setup)
1360  {
1361  // delete matvecs
1362  delete Bt_mat_vec_pt;
1363  Bt_mat_vec_pt = 0;
1364 
1365  if (!Form_BFBt_product)
1366  {
1367  delete F_mat_vec_pt;
1368  F_mat_vec_pt = 0;
1369  delete QBt_mat_vec_pt;
1370  QBt_mat_vec_pt = 0;
1371  }
1372  else
1373  {
1374  delete E_mat_vec_pt;
1375  E_mat_vec_pt = 0;
1376  }
1377 
1378  // delete stuff from displacement/position solve
1379  if (Using_default_f_preconditioner)
1380  {
1381  delete F_preconditioner_pt;
1382  F_preconditioner_pt = 0;
1383  }
1384 
1385  // delete stuff from Schur complement approx
1386  if (Using_default_p_preconditioner)
1387  {
1388  delete P_preconditioner_pt;
1389  P_preconditioner_pt = 0;
1390  }
1391  }
1392  }

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

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

◆ disable_doc_time()

void oomph::PressureBasedSolidLSCPreconditioner::disable_doc_time ( )
inline

Disable documentation of time.

265  {
266  Doc_time = false;
267  }

References Doc_time.

◆ disable_form_BFBt_product()

void oomph::PressureBasedSolidLSCPreconditioner::disable_form_BFBt_product ( )
inline

if this function is called then: in setup(...) : the matrices B, F are assembled and stored (the default behaviour) . in preconditioner_solve(...) : a sequence of matrix vector products with B, F, and Bt is performed. (Note: in this discussion no scaling was considered but B and Bt are replaced with BQ and QBt with scaling)

286  {
287  Form_BFBt_product = false;
288  }

References Form_BFBt_product.

◆ disable_p_matrix_scaling()

void oomph::PressureBasedSolidLSCPreconditioner::disable_p_matrix_scaling ( )
inline

Enable mass matrix diagonal scaling in the Schur complement approximation

197  {
198  P_matrix_using_scaling = false;
199  }

References P_matrix_using_scaling.

◆ enable_doc_time()

void oomph::PressureBasedSolidLSCPreconditioner::enable_doc_time ( )
inline

Enable documentation of time.

259  {
260  Doc_time = true;
261  }

References Doc_time.

◆ enable_form_BFBt_product()

void oomph::PressureBasedSolidLSCPreconditioner::enable_form_BFBt_product ( )
inline

If this function is called then: in setup(...) : BFBt is computed. in preconditioner_solve(...) : a single matrix vector product with BFBt is performed.

274  {
275  Form_BFBt_product = true;
276  }

References Form_BFBt_product.

◆ enable_p_matrix_scaling()

void oomph::PressureBasedSolidLSCPreconditioner::enable_p_matrix_scaling ( )
inline

Enable mass matrix diagonal scaling in the Schur complement approximation

190  {
191  P_matrix_using_scaling = true;
192  }

References P_matrix_using_scaling.

◆ is_p_matrix_using_scaling()

bool oomph::PressureBasedSolidLSCPreconditioner::is_p_matrix_using_scaling ( ) const
inline

Return whether the mass matrix is using diagonal scaling or not

204  {
205  return P_matrix_using_scaling;
206  }

References P_matrix_using_scaling.

◆ preconditioner_solve()

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

Apply preconditioner to Vector r.

Apply preconditioner to r.

Implements oomph::Preconditioner.

495  {
496 #ifdef PARANOID
497  if (Preconditioner_has_been_setup == false)
498  {
499  std::ostringstream error_message;
500  error_message << "setup must be called before using preconditioner_solve";
501  throw OomphLibError(
502  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
503  }
504  if (z.built())
505  {
506  if (z.nrow() != r.nrow())
507  {
508  std::ostringstream error_message;
509  error_message << "The vectors z and r must have the same number of "
510  << "of global rows";
511  throw OomphLibError(error_message.str(),
512  OOMPH_CURRENT_FUNCTION,
513  OOMPH_EXCEPTION_LOCATION);
514  }
515  }
516 #endif
517 
518  double t_start_overall = TimingHelpers::timer();
519  double t_start = TimingHelpers::timer();
520  double t_end = 0;
521 
522  // if z is not setup then give it the same distribution
523  if (!z.built())
524  {
525  z.build(r.distribution_pt(), 0.0);
526  }
527 
528  // Step 1 - apply approximate Schur inverse to pressure unknowns (block 1)
529  // -----------------------------------------------------------------------
530 
531  // Working vectors
532  DoubleVector temp_vec;
533  DoubleVector another_temp_vec;
534 
535  // Copy pressure values from residual vector to temp_vec:
536  // Loop over all entries in the global vector (this one
537  // includes displacement/position and pressure dofs in some random fashion)
538  this->get_block_vector(1, r, temp_vec);
539 
540 
541  if (Doc_time)
542  {
543  t_end = TimingHelpers::timer();
544  oomph_info << "LSC prec solve: Time for get block vector: "
545  << t_end - t_start << std::endl;
546  t_start = TimingHelpers::timer();
547  }
548 
549  // NOTE: The vector temp_vec now contains the vector r_p.
550 
551  // Solve first pressure Poisson system
552 #ifdef PARANOID
553  // check a solver has been set
554  if (P_preconditioner_pt == 0)
555  {
556  std::ostringstream error_message;
557  error_message << "P_preconditioner_pt has not been set.";
558  throw OomphLibError(
559  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
560  }
561 #endif
562 
563  // use some Preconditioner's preconditioner_solve function
564  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
565 
566 
567  if (Doc_time)
568  {
569  t_end = TimingHelpers::timer();
570  oomph_info << "LSC prec solve: First P solve [nrow="
571  << P_preconditioner_pt->nrow() << "]: " << t_end - t_start
572  << std::endl;
573  t_start = TimingHelpers::timer();
574  }
575 
576 
577  // NOTE: The vector another_temp_vec now contains the vector P^{-1} r_p
578 
579  // Multiply another_temp_vec by matrix E and stick the result into temp_vec
580  temp_vec.clear();
581  if (Form_BFBt_product)
582  {
583  E_mat_vec_pt->multiply(another_temp_vec, temp_vec);
584  }
585  else
586  {
587  QBt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
588  another_temp_vec.clear();
589  F_mat_vec_pt->multiply(temp_vec, another_temp_vec);
590  temp_vec.clear();
591  QBt_mat_vec_pt->multiply_transpose(another_temp_vec, temp_vec);
592  }
593 
594 
595  if (Doc_time)
596  {
597  t_end = TimingHelpers::timer();
598  oomph_info << "LSC prec solve: E matrix vector product: "
599  << t_end - t_start << std::endl;
600  t_start = TimingHelpers::timer();
601  }
602 
603  // NOTE: The vector temp_vec now contains E P^{-1} r_p
604 
605  // Solve second pressure Poisson system using preconditioner_solve
606  another_temp_vec.clear();
607  P_preconditioner_pt->preconditioner_solve(temp_vec, another_temp_vec);
608 
609 
610  if (Doc_time)
611  {
612  t_end = TimingHelpers::timer();
613  oomph_info << "LSC prec solve: Second P solve [nrow="
614  << P_preconditioner_pt->nrow() << "]: " << t_end - t_start
615  << std::endl;
616  t_start = TimingHelpers::timer();
617  }
618 
619 
620  // NOTE: The vector another_temp_vec now contains z_p = P^{-1} E P^{-1} r_p
621  // as required (apart from the sign which we'll fix in the
622  // next step.
623 
624  // Now copy another_temp_vec (i.e. z_p) back into the global vector z.
625  // Loop over all entries in the global results vector z:
626  temp_vec.build(another_temp_vec.distribution_pt(), 0.0);
627  temp_vec -= another_temp_vec;
628  return_block_vector(1, temp_vec, z);
629 
630  // Step 2 - apply preconditioner to displacement/positon unknowns (block 0)
631  // ------------------------------------------------------------------------
632 
633  // Recall that another_temp_vec (computed above) contains the
634  // negative of the solution of the Schur complement systen, -z_p.
635  // Multiply by G (stored in Block_matrix_pt(0,1) and store
636  // result in temp_vec (vector resizes itself).
637  temp_vec.clear();
638  Bt_mat_vec_pt->multiply(another_temp_vec, temp_vec);
639 
640 
641  if (Doc_time)
642  {
643  t_end = TimingHelpers::timer();
644  oomph_info << "LSC prec solve: G matrix vector product: "
645  << t_end - t_start << std::endl;
646  t_start = TimingHelpers::timer();
647  }
648 
649  // NOTE: temp_vec now contains -G z_p
650 
651  // The vector another_temp_vec is no longer needed -- re-use it to store
652  // displacement/position quantities:
653  another_temp_vec.clear();
654 
655  // Loop over all enries in the global vector and find the
656  // entries associated with the displacement/position:
657  get_block_vector(0, r, another_temp_vec);
658  another_temp_vec += temp_vec;
659 
660  // NOTE: The vector another_temp_vec now contains r_u - G z_p
661 
662  // Solve momentum system
663 #ifdef PARANOID
664  // check a solver has been set
665  if (F_preconditioner_pt == 0)
666  {
667  std::ostringstream error_message;
668  error_message << "F_preconditioner_pt has not been set.";
669  throw OomphLibError(
670  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
671  }
672 #endif
673 
674  // use some Preconditioner's preconditioner solve
675  // and return
676  if (F_preconditioner_is_block_preconditioner)
677  {
678  return_block_vector(0, another_temp_vec, z);
679  F_preconditioner_pt->preconditioner_solve(z, z);
680  }
681  else
682  {
683  F_preconditioner_pt->preconditioner_solve(another_temp_vec, temp_vec);
684  return_block_vector(0, temp_vec, z);
685  }
686 
687  if (Doc_time)
688  {
689  t_end = TimingHelpers::timer();
690  oomph_info << "LSC prec solve: F solve [nrow="
691  << P_preconditioner_pt->nrow() << "]: " << t_end - t_start
692  << std::endl;
693  oomph_info << "LSC prec solve: Overall " << t_end - t_start_overall
694  << std::endl;
695  }
696  }
DoubleVector
std::vector< double > DoubleVector
loads clump configuration
Definition: ClumpInput.h:26
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::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::Preconditioner::preconditioner_solve
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
oomph::PressureBasedSolidLSCPreconditioner::F_preconditioner_is_block_preconditioner
bool F_preconditioner_is_block_preconditioner
Definition: solid_preconditioners.h:327
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
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
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References Bt_mat_vec_pt, oomph::DoubleVector::build(), oomph::DoubleVector::built(), oomph::DoubleVector::clear(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_time, E_mat_vec_pt, F_mat_vec_pt, F_preconditioner_is_block_preconditioner, F_preconditioner_pt, Form_BFBt_product, oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vector(), oomph::MatrixVectorProduct::multiply(), oomph::MatrixVectorProduct::multiply_transpose(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, P_preconditioner_pt, Preconditioner_has_been_setup, oomph::Preconditioner::preconditioner_solve(), QBt_mat_vec_pt, UniformPSDSelfTest::r, oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vector(), and oomph::TimingHelpers::timer().

◆ set_f_preconditioner()

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

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

234  {
235  // If the default preconditioner has been used
236  // clean it up now...
237  if (Using_default_f_preconditioner)
238  {
239  delete F_preconditioner_pt;
240  }
241  F_preconditioner_pt = new_f_preconditioner_pt;
242  Using_default_f_preconditioner = false;
243  }

References F_preconditioner_pt, and Using_default_f_preconditioner.

◆ set_f_superlu_preconditioner()

void oomph::PressureBasedSolidLSCPreconditioner::set_f_superlu_preconditioner ( )
inline

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

248  {
249  if (!Using_default_f_preconditioner)
250  {
251  F_preconditioner_pt = new SuperLUPreconditioner;
252  Using_default_f_preconditioner = true;
253  }
254  }

References F_preconditioner_pt, and Using_default_f_preconditioner.

◆ set_p_preconditioner()

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

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

210  {
211  // If the default preconditioner has been used
212  // clean it up now...
213  if (Using_default_p_preconditioner)
214  {
215  delete P_preconditioner_pt;
216  }
217  P_preconditioner_pt = new_p_preconditioner_pt;
218  Using_default_p_preconditioner = false;
219  }

References P_preconditioner_pt, and Using_default_p_preconditioner.

◆ set_p_superlu_preconditioner()

void oomph::PressureBasedSolidLSCPreconditioner::set_p_superlu_preconditioner ( )
inline

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

224  {
225  if (!Using_default_p_preconditioner)
226  {
227  P_preconditioner_pt = new SuperLUPreconditioner;
228  Using_default_p_preconditioner = true;
229  }
230  }

References P_preconditioner_pt, and Using_default_p_preconditioner.

◆ set_solid_mesh()

void oomph::PressureBasedSolidLSCPreconditioner::set_solid_mesh ( Mesh *  mesh_pt)
inline

specify the mesh containing the mesh containing the block-preconditionable solid elements. The dimension of the problem must also be specified.

183  {
184  Solid_mesh_pt = mesh_pt;
185  }
oomph::BlockPreconditioner< CRDoubleMatrix >::mesh_pt
const Mesh * mesh_pt(const unsigned &i) const
Definition: block_preconditioner.h:1782

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

◆ setup()

void oomph::PressureBasedSolidLSCPreconditioner::setup ( )
virtual

Broken assignment operator.

Setup the preconditioner

Setup the least-squares commutator solid preconditioner. This extracts blocks corresponding to the position/displacement 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.

39  {
40  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41  // NOTE: In the interest of minimising memory usage, several containers
42  // are recycled, therefore their content/meaning changes
43  // throughout this function. The code is carefully annotated
44  // but you'll have to read it line by line!
45  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 
47  // make sure any old data is deleted
48  clean_up_memory();
49 
50 #ifdef PARANOID
51  // paranoid check that the solid mesh pt has been set
52  if (Solid_mesh_pt == 0)
53  {
54  std::ostringstream error_message;
55  error_message << "The solid mesh pointer must be set.\n"
56  << "Use method set_solid_mesh(...)";
57  throw OomphLibError(
58  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
59  }
60 #endif
61 
62  // set the mesh
63  this->set_mesh(0, Solid_mesh_pt);
64 
65  // Get blocks
66  // ----------
67 
68  // Set up block look up schemes (done automatically in the
69  // BlockPreconditioner base class, based on the information
70  // provided in the block-preconditionable elements in the problem)
71 
72  // this preconditioner has two types of block:
73  // type 0: displacement/positions - corresponding to DOFs 0 to n-2
74  // type 1: pressure - corresponding to DOF n-1
75  double t_block_start = TimingHelpers::timer();
76  unsigned ndof_types = 0;
77  if (this->is_subsidiary_block_preconditioner())
78  {
79  ndof_types = this->ndof_types();
80  }
81  else
82  {
83  ndof_types = this->ndof_types_in_mesh(0);
84  }
85  Vector<unsigned> dof_to_block_map(ndof_types);
86  dof_to_block_map[ndof_types - 1] = 1;
87  this->block_setup(dof_to_block_map);
88  double t_block_finish = TimingHelpers::timer();
89  double block_setup_time = t_block_finish - t_block_start;
90  if (Doc_time)
91  {
92  oomph_info << "Time for block_setup(...) [sec]: " << block_setup_time
93  << "\n";
94  }
95 
96  // determine whether the F preconditioner is a block preconditioner (and
97  // therefore a subsidiary preconditioner)
98  BlockPreconditioner<CRDoubleMatrix>* F_block_preconditioner_pt =
99  dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(F_preconditioner_pt);
100  F_preconditioner_is_block_preconditioner = true;
101  if (F_block_preconditioner_pt == 0)
102  {
103  F_preconditioner_is_block_preconditioner = false;
104  }
105 
106  // Get B (the divergence block)
107  double t_get_B_start = TimingHelpers::timer();
108  CRDoubleMatrix* b_pt = new CRDoubleMatrix;
109  this->get_block(1, 0, *b_pt);
110  double t_get_B_finish = TimingHelpers::timer();
111  if (Doc_time)
112  {
113  double get_B_time = t_get_B_finish - t_get_B_start;
114  oomph_info << "Time to get B [sec]: " << get_B_time << "\n";
115  }
116 
117  // the pointer for f
118  CRDoubleMatrix* f_pt = new CRDoubleMatrix;
119 
120  // the pointer for the p matrix
121  CRDoubleMatrix* p_matrix_pt = new CRDoubleMatrix;
122 
123  // the pointer for bt
124  CRDoubleMatrix* bt_pt = new CRDoubleMatrix;
125 
126  // if BFBt is to be formed
127  if (Form_BFBt_product)
128  {
129  // If using scaling then replace B with Bt
130  CRDoubleMatrix* ivmm_pt = 0;
131  if (P_matrix_using_scaling)
132  {
133  // Assemble the ivmm diagonal
134  double ivmm_assembly_start_t = TimingHelpers::timer();
135  ivmm_pt = assemble_mass_matrix_diagonal();
136  double ivmm_assembly_finish_t = TimingHelpers::timer();
137  if (Doc_time)
138  {
139  double ivmm_assembly_time =
140  ivmm_assembly_finish_t - ivmm_assembly_start_t;
141  oomph_info << "Time to assemble inverse mass matrix [sec]: "
142  << ivmm_assembly_time << "\n";
143  }
144 
145  // assemble BQ (stored in B)
146  double t_BQ_start = TimingHelpers::timer();
147  CRDoubleMatrix* temp_matrix_pt = new CRDoubleMatrix;
148  b_pt->multiply(*ivmm_pt, *temp_matrix_pt);
149  delete b_pt;
150  b_pt = 0;
151  b_pt = temp_matrix_pt;
152  double t_BQ_finish = TimingHelpers::timer();
153  if (Doc_time)
154  {
155  double t_BQ_time = t_BQ_finish - t_BQ_start;
156  oomph_info << "Time to generate BQ [sec]: " << t_BQ_time << std::endl;
157  }
158  }
159 
160  // Get Bt
161  double t_get_Bt_start = TimingHelpers::timer();
162  this->get_block(0, 1, *bt_pt);
163  double t_get_Bt_finish = TimingHelpers::timer();
164  if (Doc_time)
165  {
166  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
167  oomph_info << "Time to get Bt [sec]: " << t_get_Bt_time << std::endl;
168  }
169 
170  // now form the P matrix by multiplying B (which if using scaling will be
171  // BQ) with Bt
172  double t_P_start = TimingHelpers::timer();
173  b_pt->multiply(*bt_pt, *p_matrix_pt);
174  double t_P_finish = TimingHelpers::timer();
175  if (Doc_time)
176  {
177  double t_P_time = t_P_finish - t_P_start;
178  oomph_info << "Time to generate P matrix [sec]: " << t_P_time
179  << std::endl;
180  }
181 
182  // Multiply auxiliary matrix by diagonal of mass matrix if
183  // required
184  if (P_matrix_using_scaling)
185  {
186  CRDoubleMatrix* temp_matrix_pt = new CRDoubleMatrix;
187  double t_QBt_start = TimingHelpers::timer();
188  ivmm_pt->multiply(*bt_pt, *temp_matrix_pt);
189  delete bt_pt;
190  bt_pt = 0;
191  bt_pt = temp_matrix_pt;
192  double t_QBt_finish = TimingHelpers::timer();
193 
194  // Output times
195  if (Doc_time)
196  {
197  double t_QBt_time = t_QBt_finish - t_QBt_start;
198  oomph_info << "Time to generate QBt [sec]: " << t_QBt_time
199  << std::endl;
200  }
201  }
202 
203  // Clean up memory
204  delete ivmm_pt;
205 
206  // get block 0 0
207  double t_get_F_start = TimingHelpers::timer();
208  this->get_block(0, 0, *f_pt);
209  double t_get_F_finish = TimingHelpers::timer();
210  if (Doc_time)
211  {
212  double t_get_F_time = t_get_F_finish - t_get_F_start;
213  oomph_info << "Time to get F [sec]: " << t_get_F_time << std::endl;
214  }
215 
216  // Auxiliary matrix for intermediate results
217  double t_aux_matrix_start = TimingHelpers::timer();
218  CRDoubleMatrix* aux_matrix_pt = new CRDoubleMatrix;
219  f_pt->multiply(*bt_pt, *aux_matrix_pt);
220  double t_aux_matrix_finish = TimingHelpers::timer();
221  if (Doc_time)
222  {
223  double t_aux_time = t_aux_matrix_finish - t_aux_matrix_start;
224  oomph_info << "Time to generate FQBt [sec]: " << t_aux_time
225  << std::endl;
226  }
227 
228  // can now delete Bt (or QBt if using scaling)
229  delete bt_pt;
230  bt_pt = 0;
231 
232  // and if F requires a block preconditioner then we can delete F
233  if (F_preconditioner_is_block_preconditioner)
234  {
235  delete f_pt;
236  }
237 
238  // now form BFBt
239  double t_E_matrix_start = TimingHelpers::timer();
240  CRDoubleMatrix* e_matrix_pt = new CRDoubleMatrix;
241  b_pt->multiply(*aux_matrix_pt, *e_matrix_pt);
242  delete aux_matrix_pt;
243  delete b_pt;
244  double t_E_matrix_finish = TimingHelpers::timer();
245  if (Doc_time)
246  {
247  double t_E_time = t_E_matrix_finish - t_E_matrix_start;
248  oomph_info << "Time to generate E (B*(F*Bt)) [sec]: " << t_E_time
249  << std::endl;
250  }
251  double t_E_matvec_start = TimingHelpers::timer();
252  E_mat_vec_pt = new MatrixVectorProduct;
253  // E_mat_vec_pt->setup(e_matrix_pt);
254  this->setup_matrix_vector_product(E_mat_vec_pt, e_matrix_pt, 1);
255  double t_E_matvec_finish = TimingHelpers::timer();
256  delete e_matrix_pt;
257  if (Doc_time)
258  {
259  double t_E_time = t_E_matvec_finish - t_E_matvec_start;
260  oomph_info << "Time to build E (BFBt) matrix vector operator E [sec]: "
261  << t_E_time << std::endl;
262  }
263 
264  // and rebuild Bt
265  t_get_Bt_start = TimingHelpers::timer();
266  bt_pt = new CRDoubleMatrix;
267  this->get_block(0, 1, *bt_pt);
268  t_get_Bt_finish = TimingHelpers::timer();
269  if (Doc_time)
270  {
271  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
272  oomph_info << "Time to get Bt [sec]: " << t_get_Bt_time << std::endl;
273  }
274  }
275 
276 
278 
279 
280  // otherwise we are not forming BFBt
281  else
282  {
283  // get the inverse mass matrix (Q)
284  CRDoubleMatrix* ivmm_pt = 0;
285  if (P_matrix_using_scaling)
286  {
287  double ivmm_assembly_start_t = TimingHelpers::timer();
288  ivmm_pt = assemble_mass_matrix_diagonal();
289  double ivmm_assembly_finish_t = TimingHelpers::timer();
290  if (Doc_time)
291  {
292  double ivmm_assembly_time =
293  ivmm_assembly_finish_t - ivmm_assembly_start_t;
294  oomph_info << "Time to assemble Q (inverse diagonal "
295  << "mass matrix) [sec]: " << ivmm_assembly_time << "\n";
296  }
297  }
298 
299  // Get Bt
300  double t_get_Bt_start = TimingHelpers::timer();
301  this->get_block(0, 1, *bt_pt);
302  double t_get_Bt_finish = TimingHelpers::timer();
303  if (Doc_time)
304  {
305  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
306  oomph_info << "Time to get Bt [sec]: " << t_get_Bt_time << std::endl;
307  }
308 
309  // next QBt
310  if (P_matrix_using_scaling)
311  {
312  double t_QBt_matrix_start = TimingHelpers::timer();
313  CRDoubleMatrix* qbt_pt = new CRDoubleMatrix;
314  ivmm_pt->multiply(*bt_pt, *qbt_pt);
315  delete bt_pt;
316  bt_pt = 0;
317  bt_pt = qbt_pt;
318  double t_QBt_matrix_finish = TimingHelpers::timer();
319  if (Doc_time)
320  {
321  double t_QBt_time = t_QBt_matrix_finish - t_QBt_matrix_start;
322  oomph_info << "Time to generate QBt [sec]: " << t_QBt_time
323  << std::endl;
324  }
325  delete ivmm_pt;
326  }
327 
328  // form P
329  double t_p_matrix_start = TimingHelpers::timer();
330  b_pt->multiply(*bt_pt, *p_matrix_pt);
331  double t_p_matrix_finish = TimingHelpers::timer();
332  if (Doc_time)
333  {
334  double t_p_time = t_p_matrix_finish - t_p_matrix_start;
335  oomph_info << "Time to generate P [sec]: " << t_p_time << std::endl;
336  }
337  delete b_pt;
338  b_pt = 0;
339 
340  // build the matvec operator for QBt
341  double t_QBt_MV_start = TimingHelpers::timer();
342  QBt_mat_vec_pt = new MatrixVectorProduct;
343  // QBt_mat_vec_pt->setup(bt_pt);
344  this->setup_matrix_vector_product(QBt_mat_vec_pt, bt_pt, 1);
345  double t_QBt_MV_finish = TimingHelpers::timer();
346  if (Doc_time)
347  {
348  double t_p_time = t_QBt_MV_finish - t_QBt_MV_start;
349  oomph_info << "Time to build QBt matrix vector operator [sec]: "
350  << t_p_time << std::endl;
351  }
352  delete bt_pt;
353  bt_pt = 0;
354 
355  // get F
356  double t_get_F_start = TimingHelpers::timer();
357  this->get_block(0, 0, *f_pt);
358  double t_get_F_finish = TimingHelpers::timer();
359  if (Doc_time)
360  {
361  double t_get_F_time = t_get_F_finish - t_get_F_start;
362  oomph_info << "Time to get F [sec]: " << t_get_F_time << std::endl;
363  }
364 
365  // form the matrix vector product helper
366  double t_F_MV_start = TimingHelpers::timer();
367  F_mat_vec_pt = new MatrixVectorProduct;
368  // F_mat_vec_pt->setup(f_pt);
369  this->setup_matrix_vector_product(F_mat_vec_pt, f_pt, 0);
370  double t_F_MV_finish = TimingHelpers::timer();
371  if (Doc_time)
372  {
373  double t_F_MV_time = t_F_MV_finish - t_F_MV_start;
374  oomph_info << "Time to build F Matrix Vector Operator [sec]: "
375  << t_F_MV_time << std::endl;
376  }
377 
378  // if F is a block preconditioner then we can delete the F matrix
379  if (F_preconditioner_is_block_preconditioner)
380  {
381  delete f_pt;
382  f_pt = 0;
383  }
384 
385  // and rebuild Bt
386  t_get_Bt_start = TimingHelpers::timer();
387  bt_pt = new CRDoubleMatrix;
388  this->get_block(0, 1, *bt_pt);
389  t_get_Bt_finish = TimingHelpers::timer();
390  if (Doc_time)
391  {
392  double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
393  oomph_info << "Time to get Bt [sec]: " << t_get_Bt_time << std::endl;
394  }
395  }
396 
397 
399 
400 
401  // form the matrix vector operator for Bt
402  double t_Bt_MV_start = TimingHelpers::timer();
403  Bt_mat_vec_pt = new MatrixVectorProduct;
404  // Bt_mat_vec_pt->setup(bt_pt);
405  this->setup_matrix_vector_product(Bt_mat_vec_pt, bt_pt, 1);
406  double t_Bt_MV_finish = TimingHelpers::timer();
407  if (Doc_time)
408  {
409  double t_Bt_MV_time = t_Bt_MV_finish - t_Bt_MV_start;
410  oomph_info << "Time to build Bt Matrix Vector Operator [sec]: "
411  << t_Bt_MV_time << std::endl;
412  }
413  delete bt_pt;
414  bt_pt = 0;
415 
416  // if the P preconditioner has not been setup
417  if (P_preconditioner_pt == 0)
418  {
419  P_preconditioner_pt = new SuperLUPreconditioner;
420  Using_default_p_preconditioner = true;
421  }
422 
423  // std::stringstream junk;
424  // junk << "p_matrix" << comm_pt()->my_rank()
425  // << ".dat";
426  // p_matrix_pt->sparse_indexed_output_with_offset(junk.str());
427  // oomph_info << "Done output of " << junk.str() << std::endl;
428 
429  // Setup the preconditioner for the Pressure matrix
430  double t_p_prec_start = TimingHelpers::timer();
431  P_preconditioner_pt->setup(p_matrix_pt);
432  delete p_matrix_pt;
433  p_matrix_pt = 0;
434  p_matrix_pt = 0;
435  double t_p_prec_finish = TimingHelpers::timer();
436  if (Doc_time)
437  {
438  double t_p_prec_time = t_p_prec_finish - t_p_prec_start;
439  oomph_info << "P sub-preconditioner setup time [sec]: " << t_p_prec_time
440  << "\n";
441  }
442 
443  // Set up solver for solution of system with momentum matrix
444  // ----------------------------------------------------------
445 
446  // if the F preconditioner has not been setup
447  if (F_preconditioner_pt == 0)
448  {
449  F_preconditioner_pt = new SuperLUPreconditioner;
450  Using_default_f_preconditioner = true;
451  }
452 
453  // if F is a block preconditioner
454  double t_f_prec_start = TimingHelpers::timer();
455  if (F_preconditioner_is_block_preconditioner)
456  {
457  unsigned ndof_types = this->ndof_types();
458  ndof_types--;
459  Vector<unsigned> dof_map(ndof_types);
460  for (unsigned i = 0; i < ndof_types; i++)
461  {
462  dof_map[i] = i;
463  }
464  F_block_preconditioner_pt->turn_into_subsidiary_block_preconditioner(
465  this, dof_map);
466  F_block_preconditioner_pt->setup(matrix_pt());
467  }
468  // otherwise F is not a block preconditioner
469  else
470  {
471  F_preconditioner_pt->setup(f_pt);
472  delete f_pt;
473  f_pt = 0;
474  }
475  double t_f_prec_finish = TimingHelpers::timer();
476  if (Doc_time)
477  {
478  double t_f_prec_time = t_f_prec_finish - t_f_prec_start;
479  oomph_info << "F sub-preconditioner setup time [sec]: " << t_f_prec_time
480  << "\n";
481  }
482 
483  // Remember that the preconditioner has been setup so
484  // the stored information can be wiped when we
485  // come here next...
486  Preconditioner_has_been_setup = true;
487  }
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 >::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 >::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::Preconditioner::setup
void setup(DoubleMatrixBase *matrix_pt)
Definition: preconditioner.h:94
oomph::PressureBasedSolidLSCPreconditioner::assemble_mass_matrix_diagonal
CRDoubleMatrix * assemble_mass_matrix_diagonal()
Definition: solid_preconditioners.cc:705

References assemble_mass_matrix_diagonal(), oomph::BlockPreconditioner< CRDoubleMatrix >::block_setup(), Bt_mat_vec_pt, clean_up_memory(), Doc_time, E_mat_vec_pt, F_mat_vec_pt, F_preconditioner_is_block_preconditioner, F_preconditioner_pt, Form_BFBt_product, oomph::BlockPreconditioner< CRDoubleMatrix >::get_block(), i, oomph::BlockPreconditioner< CRDoubleMatrix >::is_subsidiary_block_preconditioner(), oomph::BlockPreconditioner< CRDoubleMatrix >::matrix_pt(), oomph::CRDoubleMatrix::multiply(), oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types(), oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types_in_mesh(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, P_matrix_using_scaling, P_preconditioner_pt, Preconditioner_has_been_setup, QBt_mat_vec_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::set_mesh(), oomph::Preconditioner::setup(), oomph::BlockPreconditioner< CRDoubleMatrix >::setup_matrix_vector_product(), Solid_mesh_pt, oomph::TimingHelpers::timer(), oomph::BlockPreconditioner< MATRIX >::turn_into_subsidiary_block_preconditioner(), Using_default_f_preconditioner, and Using_default_p_preconditioner.

Member Data Documentation

◆ Bt_mat_vec_pt

MatrixVectorProduct* oomph::PressureBasedSolidLSCPreconditioner::Bt_mat_vec_pt
private

MatrixVectorProduct operator for Bt;.

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

◆ Doc_time

bool oomph::PressureBasedSolidLSCPreconditioner::Doc_time
private

Set Doc_time to true for outputting results of timings.

Referenced by disable_doc_time(), enable_doc_time(), preconditioner_solve(), PressureBasedSolidLSCPreconditioner(), and setup().

◆ E_mat_vec_pt

MatrixVectorProduct* oomph::PressureBasedSolidLSCPreconditioner::E_mat_vec_pt
private

MatrixVectorProduct operator for E (BFBt) if BFBt is to be formed.

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

◆ F_mat_vec_pt

MatrixVectorProduct* oomph::PressureBasedSolidLSCPreconditioner::F_mat_vec_pt
private

MatrixVectorProduct operator for F if BFBt is not to be formed.

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

◆ F_preconditioner_is_block_preconditioner

bool oomph::PressureBasedSolidLSCPreconditioner::F_preconditioner_is_block_preconditioner
private

Boolean indicating whether the momentum system preconditioner is a block preconditioner

Referenced by preconditioner_solve(), and setup().

◆ F_preconditioner_pt

Preconditioner* oomph::PressureBasedSolidLSCPreconditioner::F_preconditioner_pt
private

Pointer to the 'preconditioner' for the F matrix.

Referenced by clean_up_memory(), preconditioner_solve(), PressureBasedSolidLSCPreconditioner(), set_f_preconditioner(), set_f_superlu_preconditioner(), and setup().

◆ Form_BFBt_product

bool oomph::PressureBasedSolidLSCPreconditioner::Form_BFBt_product
private

indicates whether BFBt should be formed or the component matrices should be retained. If true then: in setup(...) : BFBt is computed. in preconditioner_solve(...) : a single matrix vector product with BFBt is performed. if false then: in setup(...) : the matrices B, F are assembled and stored. in preconditioner_solve(...) : a sequence of matrix vector products with B, F, and Bt is performed. (Note: in this discussion no scaling was considered but B and Bt are replaced with BQ and QBt with scaling)

Referenced by clean_up_memory(), disable_form_BFBt_product(), enable_form_BFBt_product(), preconditioner_solve(), PressureBasedSolidLSCPreconditioner(), and setup().

◆ P_matrix_using_scaling

bool oomph::PressureBasedSolidLSCPreconditioner::P_matrix_using_scaling
private

Control flag is true if mass matrix diagonal scaling is used in the Schur complement approximation

Referenced by disable_p_matrix_scaling(), enable_p_matrix_scaling(), is_p_matrix_using_scaling(), PressureBasedSolidLSCPreconditioner(), and setup().

◆ P_preconditioner_pt

Preconditioner* oomph::PressureBasedSolidLSCPreconditioner::P_preconditioner_pt
private

Pointer to the 'preconditioner' for the pressure matrix.

Referenced by clean_up_memory(), preconditioner_solve(), PressureBasedSolidLSCPreconditioner(), set_p_preconditioner(), set_p_superlu_preconditioner(), and setup().

◆ Preconditioner_has_been_setup

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

◆ QBt_mat_vec_pt

MatrixVectorProduct* oomph::PressureBasedSolidLSCPreconditioner::QBt_mat_vec_pt
private

MatrixVectorProduct operator for QBt if BFBt is not to be formed.

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

◆ Solid_mesh_pt

Mesh* oomph::PressureBasedSolidLSCPreconditioner::Solid_mesh_pt
private

the pointer to the mesh of block preconditionable solid elements.

Referenced by assemble_mass_matrix_diagonal(), PressureBasedSolidLSCPreconditioner(), set_solid_mesh(), and setup().

◆ Using_default_f_preconditioner

bool oomph::PressureBasedSolidLSCPreconditioner::Using_default_f_preconditioner
private

flag indicating whether the default F preconditioner is used

Referenced by clean_up_memory(), PressureBasedSolidLSCPreconditioner(), set_f_preconditioner(), set_f_superlu_preconditioner(), and setup().

◆ Using_default_p_preconditioner

bool oomph::PressureBasedSolidLSCPreconditioner::Using_default_p_preconditioner
private

flag indicating whether the default P preconditioner is used

Referenced by clean_up_memory(), PressureBasedSolidLSCPreconditioner(), set_p_preconditioner(), set_p_superlu_preconditioner(), and setup().

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