Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
oomph::HelmholtzMGPreconditioner< DIM > Class Template Reference

#include <helmholtz_geometric_multigrid.h>

+ Inheritance diagram for oomph::HelmholtzMGPreconditioner< DIM >:

Public Types

typedef HelmholtzSmoother *(* PreSmootherFactoryFctPt) ()
 
typedef HelmholtzSmoother *(* PostSmootherFactoryFctPt) ()
 

Public Member Functions

void set_pre_smoother_factory_function (PreSmootherFactoryFctPt pre_smoother_fn)
 Access function to set the pre-smoother creation function. More...
 
void set_post_smoother_factory_function (PostSmootherFactoryFctPt post_smoother_fn)
 Access function to set the post-smoother creation function. More...
 
 HelmholtzMGPreconditioner (HelmholtzMGProblem *mg_problem_pt)
 
 ~HelmholtzMGPreconditioner ()
 Delete any dynamically allocated data. More...
 
void clean_up_memory ()
 Clean up anything that needs to be cleaned up. More...
 
double & tolerance ()
 Access function for the variable Tolerance (lvalue) More...
 
double & alpha_shift ()
 Function to change the value of the shift. More...
 
void disable_doc_time ()
 Disable time documentation. More...
 
void disable_v_cycle_output ()
 
void disable_output ()
 
void enable_doc_time ()
 Enable time documentation. More...
 
void enable_v_cycle_output ()
 Enable the output of the V-cycle timings and other output. More...
 
void enable_output ()
 Enable the output from anything that could have been suppressed. More...
 
void disable_smoother_and_superlu_doc_time ()
 Suppress the output of both smoothers and SuperLU. More...
 
unsigned & npost_smooth ()
 Return the number of post-smoothing iterations (lvalue) More...
 
unsigned & npre_smooth ()
 Return the number of pre-smoothing iterations (lvalue) More...
 
void pre_smooth (const unsigned &level)
 
void post_smooth (const unsigned &level)
 
double residual_norm (const unsigned &level)
 
double residual_norm (const unsigned &level, Vector< DoubleVector > &residual)
 Calculate the norm of the residual vector, r=b-Ax. More...
 
void setup_coarsest_level_structures ()
 
void direct_solve ()
 Call the direct solver (SuperLU) to solve the problem exactly. More...
 
void interpolation_matrix_set (const unsigned &level, double *value, int *col_index, int *row_st, unsigned &ncol, unsigned &nnz)
 
void interpolation_matrix_set (const unsigned &level, Vector< double > &value, Vector< int > &col_index, Vector< int > &row_st, unsigned &ncol, unsigned &nrow)
 
void set_restriction_matrices_as_interpolation_transposes ()
 
void restrict_residual (const unsigned &level)
 
void interpolate_and_correct (const unsigned &level)
 
void level_up_local_coord_of_node (const int &son_type, Vector< double > &s)
 
void setup_interpolation_matrices ()
 Setup the interpolation matrix on each level. More...
 
void setup_interpolation_matrices_unstructured ()
 Setup the interpolation matrices. More...
 
void setup_transfer_matrices ()
 Setup the transfer matrices on each level. More...
 
void full_setup ()
 Runs a full setup of the MG solver. More...
 
void preconditioner_solve (const DoubleVector &r, DoubleVector &z)
 Function applies MG to the vector r for a full solve. More...
 
unsigned iterations () const
 Number of iterations. More...
 
void level_up_local_coord_of_node (const int &son_type, Vector< double > &s)
 
void level_up_local_coord_of_node (const int &son_type, Vector< double > &s)
 
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 mg_solve (Vector< DoubleVector > &result)
 
void block_preconditioner_self_test ()
 
void setup ()
 
void setup_mg_hierarchy ()
 
void setup_mg_structures ()
 Set up the MG structures on each level. More...
 
void maximum_edge_width (const unsigned &level, double &h)
 
void setup_smoothers ()
 Set up the smoothers on all levels. More...
 
void maximum_edge_width (const unsigned &level, double &h)
 
void maximum_edge_width (const unsigned &level, double &h)
 

Private Attributes

PreSmootherFactoryFctPt Pre_smoother_factory_function_pt
 Function to create pre-smoothers. More...
 
PostSmootherFactoryFctPt Post_smoother_factory_function_pt
 Function to create post-smoothers. More...
 
HelmholtzMGProblem * Mg_problem_pt
 Pointer to the MG problem (deep copy) More...
 
Vector< HelmholtzMGProblem * > Mg_hierarchy_pt
 Vector containing pointers to problems in hierarchy. More...
 
Vector< Vector< CRDoubleMatrix * > > Mg_matrices_storage_pt
 
CRDoubleMatrix * Coarsest_matrix_mg_pt
 
DoubleVector Coarsest_x_mg
 
DoubleVector Coarsest_rhs_mg
 
Vector< CRDoubleMatrix * > Interpolation_matrices_storage_pt
 Vector to store the interpolation matrices. More...
 
Vector< CRDoubleMatrix * > Restriction_matrices_storage_pt
 Vector to store the restriction matrices. More...
 
Vector< Vector< DoubleVector > > X_mg_vectors_storage
 
Vector< Vector< DoubleVector > > Rhs_mg_vectors_storage
 
Vector< Vector< DoubleVector > > Residual_mg_vectors_storage
 
Vector< HelmholtzSmoother * > Pre_smoothers_storage_pt
 Vector to store the pre-smoothers. More...
 
Vector< HelmholtzSmoother * > Post_smoothers_storage_pt
 Vector to store the post-smoothers. More...
 
Vector< double > Max_edge_width
 
double Wavenumber
 The value of the wavenumber, k. More...
 
double Tolerance
 The tolerance of the multigrid preconditioner. More...
 
unsigned Nlevel
 The number of levels in the multigrid heirachy. More...
 
unsigned Npre_smooth
 Number of pre-smoothing steps. More...
 
unsigned Npost_smooth
 Number of post-smoothing steps. More...
 
unsigned Nvcycle
 Maximum number of V-cycles. More...
 
unsigned V_cycle_counter
 Pointer to counter for V-cycles. More...
 
bool Doc_time
 Indicates whether or not time documentation should be used. More...
 
bool Suppress_v_cycle_output
 Indicates whether or not the V-cycle output should be suppressed. More...
 
bool Suppress_all_output
 If this is set to true then all output from the solver is suppressed. More...
 
bool Has_been_setup
 Boolean variable to indicate whether or not the solver has been setup. More...
 
bool Has_been_solved
 
std::ostream * Stream_pt
 Pointer to the output stream – defaults to std::cout. More...
 
double Alpha_shift
 Temporary version of the shift – to run bash scripts. More...
 

Additional Inherited Members

- Protected Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
void set_nmesh (const unsigned &n)
 
void set_mesh (const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
 
void set_replacement_dof_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)
 
bool any_mesh_distributed () const
 
int internal_dof_number (const unsigned &i_dof) const
 
unsigned internal_index_in_dof (const unsigned &i_dof) const
 
unsigned internal_block_dimension (const unsigned &b) const
 
unsigned internal_dof_block_dimension (const unsigned &i) const
 
unsigned master_nrow () const
 
unsigned internal_master_dof_number (const unsigned &b) const
 
const LinearAlgebraDistribution * internal_preconditioner_matrix_distribution_pt () const
 
const LinearAlgebraDistribution * preconditioner_matrix_distribution_pt () const
 
- Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void clear_distribution ()
 
- Protected Attributes inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
MapMatrix< unsigned, CRDoubleMatrix * > Replacement_dof_block_pt
 The replacement dof-level blocks. More...
 
Vector< LinearAlgebraDistribution * > Block_distribution_pt
 The distribution for the blocks. More...
 
Vector< Vector< unsigned > > Block_to_dof_map_coarse
 
Vector< Vector< unsigned > > Block_to_dof_map_fine
 Mapping for the block types to the most fine grain dof types. More...
 
Vector< Vector< unsigned > > Doftype_coarsen_map_coarse
 
Vector< Vector< unsigned > > Doftype_coarsen_map_fine
 
Vector< LinearAlgebraDistribution * > Internal_block_distribution_pt
 Storage for the default distribution for each internal block. More...
 
Vector< LinearAlgebraDistribution * > Dof_block_distribution_pt
 
Vector< unsigned > Allow_multiple_element_type_in_mesh
 
Vector< const Mesh * > Mesh_pt
 
Vector< unsigned > Ndof_types_in_mesh
 
unsigned Internal_nblock_types
 
unsigned Internal_ndof_types
 
- Protected Attributes inherited from oomph::Preconditioner
bool Silent_preconditioner_setup
 Boolean to indicate whether or not the build should be done silently. More...
 
std::ostream * Stream_pt
 Pointer to the output stream – defaults to std::cout. More...
 

Detailed Description

template<unsigned DIM>
class oomph::HelmholtzMGPreconditioner< DIM >

/////////////////////////////////////////////////////// /////////////////////////////////////////////////////// ///////////////////////////////////////////////////////

Member Typedef Documentation

◆ PostSmootherFactoryFctPt

template<unsigned DIM>
typedef HelmholtzSmoother*(* oomph::HelmholtzMGPreconditioner< DIM >::PostSmootherFactoryFctPt) ()

typedef for a function that returns a pointer to an object of the class HelmholtzSmoother to be used as the post-smoother

◆ PreSmootherFactoryFctPt

template<unsigned DIM>
typedef HelmholtzSmoother*(* oomph::HelmholtzMGPreconditioner< DIM >::PreSmootherFactoryFctPt) ()

typedef for a function that returns a pointer to an object of the class HelmholtzSmoother to be used as the pre-smoother

Constructor & Destructor Documentation

◆ HelmholtzMGPreconditioner()

template<unsigned DIM>
oomph::HelmholtzMGPreconditioner< DIM >::HelmholtzMGPreconditioner ( HelmholtzMGProblem *  mg_problem_pt)
inline

Constructor: Set up default values for number of V-cycles and pre- and post-smoothing steps.

112  : BlockPreconditioner<CRDoubleMatrix>(),
113  Pre_smoother_factory_function_pt(0),
114  Post_smoother_factory_function_pt(0),
115  Mg_problem_pt(mg_problem_pt),
116  Tolerance(1.0e-09),
117  Npre_smooth(2),
118  Npost_smooth(2),
119  Nvcycle(1),
120  Doc_time(true),
121  Suppress_v_cycle_output(false),
122  Suppress_all_output(false),
123  Has_been_setup(false),
124  Has_been_solved(false),
125  Stream_pt(0),
126  Alpha_shift(0.0)
127  {
128  // Set the pointer to the finest level as the first
129  // entry in Mg_hierarchy_pt
130  Mg_hierarchy_pt.push_back(Mg_problem_pt);
131  } // End of HelmholtzMGPreconditioner
oomph::HelmholtzMGPreconditioner::Pre_smoother_factory_function_pt
PreSmootherFactoryFctPt Pre_smoother_factory_function_pt
Function to create pre-smoothers.
Definition: helmholtz_geometric_multigrid.h:551
oomph::HelmholtzMGPreconditioner::Tolerance
double Tolerance
The tolerance of the multigrid preconditioner.
Definition: helmholtz_geometric_multigrid.h:683
oomph::HelmholtzMGPreconditioner::Has_been_setup
bool Has_been_setup
Boolean variable to indicate whether or not the solver has been setup.
Definition: helmholtz_geometric_multigrid.h:710
oomph::HelmholtzMGPreconditioner::Npre_smooth
unsigned Npre_smooth
Number of pre-smoothing steps.
Definition: helmholtz_geometric_multigrid.h:689
oomph::HelmholtzMGPreconditioner::Post_smoother_factory_function_pt
PostSmootherFactoryFctPt Post_smoother_factory_function_pt
Function to create post-smoothers.
Definition: helmholtz_geometric_multigrid.h:554
oomph::HelmholtzMGPreconditioner::Doc_time
bool Doc_time
Indicates whether or not time documentation should be used.
Definition: helmholtz_geometric_multigrid.h:701
oomph::HelmholtzMGPreconditioner::Npost_smooth
unsigned Npost_smooth
Number of post-smoothing steps.
Definition: helmholtz_geometric_multigrid.h:692
oomph::HelmholtzMGPreconditioner::Mg_hierarchy_pt
Vector< HelmholtzMGProblem * > Mg_hierarchy_pt
Vector containing pointers to problems in hierarchy.
Definition: helmholtz_geometric_multigrid.h:607
oomph::HelmholtzMGPreconditioner::Suppress_v_cycle_output
bool Suppress_v_cycle_output
Indicates whether or not the V-cycle output should be suppressed.
Definition: helmholtz_geometric_multigrid.h:704
oomph::HelmholtzMGPreconditioner::Alpha_shift
double Alpha_shift
Temporary version of the shift – to run bash scripts.
Definition: helmholtz_geometric_multigrid.h:720
oomph::HelmholtzMGPreconditioner::Stream_pt
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout.
Definition: helmholtz_geometric_multigrid.h:717
oomph::HelmholtzMGPreconditioner::Mg_problem_pt
HelmholtzMGProblem * Mg_problem_pt
Pointer to the MG problem (deep copy)
Definition: helmholtz_geometric_multigrid.h:604
oomph::HelmholtzMGPreconditioner::Suppress_all_output
bool Suppress_all_output
If this is set to true then all output from the solver is suppressed.
Definition: helmholtz_geometric_multigrid.h:707
oomph::HelmholtzMGPreconditioner::Has_been_solved
bool Has_been_solved
Definition: helmholtz_geometric_multigrid.h:714
oomph::HelmholtzMGPreconditioner::Nvcycle
unsigned Nvcycle
Maximum number of V-cycles.
Definition: helmholtz_geometric_multigrid.h:695

References oomph::HelmholtzMGPreconditioner< DIM >::Mg_hierarchy_pt, and oomph::HelmholtzMGPreconditioner< DIM >::Mg_problem_pt.

◆ ~HelmholtzMGPreconditioner()

template<unsigned DIM>
oomph::HelmholtzMGPreconditioner< DIM >::~HelmholtzMGPreconditioner ( )
inline

Delete any dynamically allocated data.

135  {
136  // Run the function written to clean up the memory
137  clean_up_memory();
138  } // End of ~HelmholtzMGPreconditioner
oomph::HelmholtzMGPreconditioner::clean_up_memory
void clean_up_memory()
Clean up anything that needs to be cleaned up.
Definition: helmholtz_geometric_multigrid.h:141

References oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory().

Member Function Documentation

◆ alpha_shift()

template<unsigned DIM>
double& oomph::HelmholtzMGPreconditioner< DIM >::alpha_shift ( )
inline

Function to change the value of the shift.

212  {
213  // Return the alpha shift value
214  return Alpha_shift;
215  } // End of alpha_shift

References oomph::HelmholtzMGPreconditioner< DIM >::Alpha_shift.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ block_preconditioner_self_test()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::block_preconditioner_self_test
private

Function to ensure the block form of the Jacobian matches the form described, i.e. we should have: |--—|---—| | A_r | -A_c | A = |--—|---—|. | A_c | A_r | |--—|---—|

Check the block preconditioner framework returns the correct system matrix

874  {
875  // Start clock
876  clock_t t_bl_start = clock();
877 
878 #ifdef PARANOID
879  if (Mg_hierarchy_pt[0]->mesh_pt() == 0)
880  {
881  std::stringstream err;
882  err << "Please set pointer to mesh using set_bulk_helmholtz_mesh(...).\n";
883  throw OomphLibError(
884  err.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
885  }
886 #endif
887 
888  // The preconditioner works with one mesh; set it! Since we only use the
889  // block preconditioner on the finest level, we use the mesh from that level
890  this->set_nmesh(1);
891 
892  // Elements in actual pml layer are trivially wrapped versions of
893  // their bulk counterparts. Technically they are different
894  // elements so we have to allow different element types
895  bool allow_different_element_types_in_mesh = true;
896  this->set_mesh(
897  0, Mg_problem_pt->mesh_pt(), allow_different_element_types_in_mesh);
898 
899 #ifdef PARANOID
900  // This preconditioner only works for 2 dof types
901  unsigned n_dof_types = this->ndof_types();
902  if (n_dof_types != 2)
903  {
904  std::stringstream tmp;
905  tmp << "This preconditioner only works for problems with 2 dof types\n"
906  << "Yours has " << n_dof_types;
907  throw OomphLibError(
908  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
909  }
910 #endif
911 
912  // Set up the generic block look up scheme
913  this->block_setup();
914 
915  // Extract the number of blocks.
916  unsigned nblock_types = this->nblock_types();
917 #ifdef PARANOID
918  if (nblock_types != 2)
919  {
920  std::stringstream tmp;
921  tmp << "There are supposed to be two block types.\n"
922  << "Yours has " << nblock_types;
923  throw OomphLibError(
924  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
925  }
926 #endif
927 
928  // This is how the 2x2 block matrices are extracted. We retain the sanity
929  // check (i.e. the diagonals are the same and the off-diagonals are
930  // negatives of each other in PARANOID mode. Otherwise we only extract 2
931  // matrices
932  DenseMatrix<CRDoubleMatrix*> block_pt(nblock_types, nblock_types, 0);
933  for (unsigned i = 0; i < nblock_types; i++)
934  {
935  for (unsigned j = 0; j < nblock_types; j++)
936  {
937  // we want...
938  block_pt(i, j) = new CRDoubleMatrix;
939  this->get_block(i, j, *block_pt(i, j));
940  }
941  }
942 
943  // Check that diagonal matrices are the same
944  //------------------------------------------
945  {
946  unsigned nnz1 = block_pt(0, 0)->nnz();
947  unsigned nnz2 = block_pt(1, 1)->nnz();
948  if (nnz1 != nnz2)
949  {
950  std::stringstream tmp;
951  tmp << "nnz of diagonal blocks doesn't match: " << nnz1
952  << " != " << nnz2 << std::endl;
953  throw OomphLibError(
954  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
955  }
956  unsigned nrow1 = block_pt(0, 0)->nrow();
957  unsigned nrow2 = block_pt(1, 1)->nrow();
958  if (nrow1 != nrow2)
959  {
960  std::stringstream tmp;
961  tmp << "nrow of diagonal blocks doesn't match: " << nrow1
962  << " != " << nrow2 << std::endl;
963  throw OomphLibError(
964  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
965  }
966 
967  // Check entries
968  bool fail = false;
969  double tol = 1.0e-15;
970  std::stringstream tmp;
971 
972  // Check row starts
973  for (unsigned i = 0; i < nrow1 + 1; i++)
974  {
975  if (block_pt(0, 0)->row_start()[i] != block_pt(1, 1)->row_start()[i])
976  {
977  fail = true;
978  tmp << "Row starts of diag matrices don't match for row " << i
979  << " : " << block_pt(0, 0)->row_start()[i] << " "
980  << block_pt(1, 1)->row_start()[i] << " " << std::endl;
981  }
982  }
983 
984  // Check values and column indices
985  for (unsigned i = 0; i < nnz1; i++)
986  {
987  if (block_pt(0, 0)->column_index()[i] !=
988  block_pt(1, 1)->column_index()[i])
989  {
990  fail = true;
991  tmp << "Column of diag matrices indices don't match for entry " << i
992  << " : " << block_pt(0, 0)->column_index()[i] << " "
993  << block_pt(1, 1)->column_index()[i] << " " << std::endl;
994  }
995 
996  if (fabs(block_pt(0, 0)->value()[i] - block_pt(1, 1)->value()[i]) > tol)
997  {
998  fail = true;
999  tmp << "Values of diag matrices don't match for entry " << i << " : "
1000  << block_pt(0, 0)->value()[i] << " " << block_pt(1, 1)->value()[i]
1001  << " " << std::endl;
1002  }
1003  }
1004  if (fail)
1005  {
1006  throw OomphLibError(
1007  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1008  }
1009  }
1010 
1011  // Check that off-diagonal matrices are negatives
1012  //-----------------------------------------------
1013  {
1014  unsigned nnz1 = block_pt(0, 1)->nnz();
1015  unsigned nnz2 = block_pt(1, 0)->nnz();
1016  if (nnz1 != nnz2)
1017  {
1018  std::stringstream tmp;
1019  tmp << "nnz of diagonal blocks doesn't match: " << nnz1
1020  << " != " << nnz2 << std::endl;
1021  throw OomphLibError(
1022  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1023  }
1024  unsigned nrow1 = block_pt(0, 1)->nrow();
1025  unsigned nrow2 = block_pt(1, 0)->nrow();
1026  if (nrow1 != nrow2)
1027  {
1028  std::stringstream tmp;
1029  tmp << "nrow of off-diagonal blocks doesn't match: " << nrow1
1030  << " != " << nrow2 << std::endl;
1031  throw OomphLibError(
1032  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1033  }
1034 
1035 
1036  // Check entries
1037  bool fail = false;
1038  double tol = 1.0e-15;
1039  std::stringstream tmp;
1040 
1041  // Check row starts
1042  for (unsigned i = 0; i < nrow1 + 1; i++)
1043  {
1044  if (block_pt(0, 1)->row_start()[i] != block_pt(1, 0)->row_start()[i])
1045  {
1046  fail = true;
1047  tmp << "Row starts of off-diag matrices don't match for row " << i
1048  << " : " << block_pt(0, 1)->row_start()[i] << " "
1049  << block_pt(1, 0)->row_start()[i] << " " << std::endl;
1050  }
1051  }
1052 
1053  // Check values and column indices
1054  for (unsigned i = 0; i < nnz1; i++)
1055  {
1056  if (block_pt(0, 1)->column_index()[i] !=
1057  block_pt(1, 0)->column_index()[i])
1058  {
1059  fail = true;
1060  tmp << "Column indices of off-diag matrices don't match for entry "
1061  << i << " : " << block_pt(0, 1)->column_index()[i] << " "
1062  << block_pt(1, 0)->column_index()[i] << " " << std::endl;
1063  }
1064 
1065  if (fabs(block_pt(0, 1)->value()[i] + block_pt(1, 0)->value()[i]) > tol)
1066  {
1067  fail = true;
1068  tmp << "Values of off-diag matrices aren't negatives of "
1069  << "each other for entry " << i << " : "
1070  << block_pt(0, 1)->value()[i] << " " << block_pt(1, 0)->value()[i]
1071  << " " << std::endl;
1072  }
1073  }
1074  if (fail)
1075  {
1076  throw OomphLibError(
1077  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1078  }
1079  }
1080 
1081  // Clean up
1082  for (unsigned i = 0; i < nblock_types; i++)
1083  {
1084  for (unsigned j = 0; j < nblock_types; j++)
1085  {
1086  delete block_pt(i, j);
1087  block_pt(i, j) = 0;
1088  }
1089  }
1090 
1091  // Stop clock
1092  clock_t t_bl_end = clock();
1093  double total_setup_time = double(t_bl_end - t_bl_start) / CLOCKS_PER_SEC;
1094  oomph_info << "CPU time for block preconditioner self-test [sec]: "
1095  << total_setup_time << "\n"
1096  << std::endl;
1097  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
oomph::BlockPreconditioner< CRDoubleMatrix >::mesh_pt
const Mesh * mesh_pt(const unsigned &i) const
Definition: block_preconditioner.h:1782
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 >::nblock_types
unsigned nblock_types() const
Return the number of block types.
Definition: block_preconditioner.h:1670
oomph::BlockPreconditioner< CRDoubleMatrix >::ndof_types
unsigned ndof_types() const
Return the total number of DOF types.
Definition: block_preconditioner.h:1694
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::Problem::mesh_pt
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
tmp
Eigen::Matrix< Scalar, Dynamic, Dynamic, ColMajor > tmp
Definition: level3_impl.h:365
Eigen::value
squared absolute value
Definition: GlobalFunctions.h:87
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::oomph_info
OomphInfo oomph_info
Definition: oomph_definitions.cc:319
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 boost::multiprecision::fabs(), i, j, oomph::DenseMatrix< T >::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, tmp, and Eigen::value.

◆ clean_up_memory()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory ( )
inlinevirtual

Clean up anything that needs to be cleaned up.

Reimplemented from oomph::Preconditioner.

142  {
143  // We only need to destroy data if the solver has been set up and
144  // the data hasn't already been cleared
145  if (Has_been_setup)
146  {
147  // Loop over all of the levels in the hierarchy
148  for (unsigned i = 0; i < Nlevel - 1; i++)
149  {
150  // Delete the pre-smoother associated with this level
151  delete Pre_smoothers_storage_pt[i];
152 
153  // Make it a null pointer
154  Pre_smoothers_storage_pt[i] = 0;
155 
156  // Delete the post-smoother associated with this level
157  delete Post_smoothers_storage_pt[i];
158 
159  // Make it a null pointer
160  Post_smoothers_storage_pt[i] = 0;
161 
162  // Loop over the real and imaginary parts of the system matrix
163  // associated with the i-th level
164  for (unsigned j = 0; j < 2; j++)
165  {
166  // Delete the real/imaginary part of the system matrix
167  delete Mg_matrices_storage_pt[i][j];
168 
169  // Make it a null pointer
170  Mg_matrices_storage_pt[i][j] = 0;
171  }
172  }
173 
174  // Delete the expanded matrix associated with the problem on the
175  // coarsest level
176  delete Coarsest_matrix_mg_pt;
177 
178  // Make it a null pointer
179  Coarsest_matrix_mg_pt = 0;
180 
181  // Loop over all but the coarsest of the levels in the hierarchy
182  for (unsigned i = 0; i < Nlevel - 1; i++)
183  {
184  // Delete the interpolation matrix associated with the i-th level
185  delete Interpolation_matrices_storage_pt[i];
186 
187  // Make it a null pointer
188  Interpolation_matrices_storage_pt[i] = 0;
189 
190  // Delete the restriction matrix associated with the i-th level
191  delete Restriction_matrices_storage_pt[i];
192 
193  // Make it a null pointer
194  Restriction_matrices_storage_pt[i] = 0;
195  }
196 
197  // Everything has been deleted now so we need to indicate that the
198  // solver is not set up
199  Has_been_setup = false;
200  }
201  } // End of clean_up_memory
oomph::HelmholtzMGPreconditioner::Coarsest_matrix_mg_pt
CRDoubleMatrix * Coarsest_matrix_mg_pt
Definition: helmholtz_geometric_multigrid.h:625
oomph::HelmholtzMGPreconditioner::Restriction_matrices_storage_pt
Vector< CRDoubleMatrix * > Restriction_matrices_storage_pt
Vector to store the restriction matrices.
Definition: helmholtz_geometric_multigrid.h:651
oomph::HelmholtzMGPreconditioner::Interpolation_matrices_storage_pt
Vector< CRDoubleMatrix * > Interpolation_matrices_storage_pt
Vector to store the interpolation matrices.
Definition: helmholtz_geometric_multigrid.h:648
oomph::HelmholtzMGPreconditioner::Pre_smoothers_storage_pt
Vector< HelmholtzSmoother * > Pre_smoothers_storage_pt
Vector to store the pre-smoothers.
Definition: helmholtz_geometric_multigrid.h:669
oomph::HelmholtzMGPreconditioner::Post_smoothers_storage_pt
Vector< HelmholtzSmoother * > Post_smoothers_storage_pt
Vector to store the post-smoothers.
Definition: helmholtz_geometric_multigrid.h:672
oomph::HelmholtzMGPreconditioner::Mg_matrices_storage_pt
Vector< Vector< CRDoubleMatrix * > > Mg_matrices_storage_pt
Definition: helmholtz_geometric_multigrid.h:615
oomph::HelmholtzMGPreconditioner::Nlevel
unsigned Nlevel
The number of levels in the multigrid heirachy.
Definition: helmholtz_geometric_multigrid.h:686

References oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_matrix_mg_pt, oomph::HelmholtzMGPreconditioner< DIM >::Has_been_setup, i, oomph::HelmholtzMGPreconditioner< DIM >::Interpolation_matrices_storage_pt, j, oomph::HelmholtzMGPreconditioner< DIM >::Mg_matrices_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::Nlevel, oomph::HelmholtzMGPreconditioner< DIM >::Post_smoothers_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoothers_storage_pt, and oomph::HelmholtzMGPreconditioner< DIM >::Restriction_matrices_storage_pt.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::~HelmholtzMGPreconditioner().

◆ direct_solve()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::direct_solve ( )
inline

Call the direct solver (SuperLU) to solve the problem exactly.

376  {
377  // Concatenate the vectors in X_mg into one vector, coarsest_x_mg
378  DoubleVectorHelpers::concatenate(X_mg_vectors_storage[Nlevel - 1],
379  Coarsest_x_mg);
380 
381  // Concatenate the vectors in Rhs_mg into one vector, coarsest_rhs_mg
382  DoubleVectorHelpers::concatenate(Rhs_mg_vectors_storage[Nlevel - 1],
383  Coarsest_rhs_mg);
384 
385  // Get solution by direct solve:
386  Coarsest_matrix_mg_pt->solve(Coarsest_rhs_mg, Coarsest_x_mg);
387 
388  // Split the solution vector into a vector of DoubleVectors and store it
389  DoubleVectorHelpers::split(Coarsest_x_mg,
390  X_mg_vectors_storage[Nlevel - 1]);
391  } // End of direct_solve
oomph::DoubleMatrixBase::solve
void solve(DoubleVector &rhs)
Definition: matrices.cc:50
oomph::HelmholtzMGPreconditioner::Coarsest_rhs_mg
DoubleVector Coarsest_rhs_mg
Definition: helmholtz_geometric_multigrid.h:645
oomph::HelmholtzMGPreconditioner::X_mg_vectors_storage
Vector< Vector< DoubleVector > > X_mg_vectors_storage
Definition: helmholtz_geometric_multigrid.h:658
oomph::HelmholtzMGPreconditioner::Coarsest_x_mg
DoubleVector Coarsest_x_mg
Definition: helmholtz_geometric_multigrid.h:635
oomph::HelmholtzMGPreconditioner::Rhs_mg_vectors_storage
Vector< Vector< DoubleVector > > Rhs_mg_vectors_storage
Definition: helmholtz_geometric_multigrid.h:662
oomph::DoubleVectorHelpers::split
void split(const DoubleVector &in_vector, Vector< DoubleVector * > &out_vector_pt)
Definition: double_vector.cc:1413
oomph::DoubleVectorHelpers::concatenate
void concatenate(const Vector< DoubleVector * > &in_vector_pt, DoubleVector &out_vector)
Definition: double_vector.cc:993

References oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_matrix_mg_pt, oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_rhs_mg, oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_x_mg, oomph::DoubleVectorHelpers::concatenate(), oomph::HelmholtzMGPreconditioner< DIM >::Nlevel, oomph::HelmholtzMGPreconditioner< DIM >::Rhs_mg_vectors_storage, oomph::DoubleMatrixBase::solve(), oomph::DoubleVectorHelpers::split(), and oomph::HelmholtzMGPreconditioner< DIM >::X_mg_vectors_storage.

◆ disable_doc_time()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::disable_doc_time ( )
inline

Disable time documentation.

219  {
220  // Set the value of Doc_time to false
221  Doc_time = false;
222  } // End of disable_doc_time

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ disable_output()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::disable_output ( )
inline

Suppress anything that can be suppressed, i.e. any timings. Things like mesh adaptation can not however be silenced using this

238  {
239  // Set the value of Doc_time to false
240  Doc_time = false;
241 
242  // Enable the suppression of output from the V-cycle
243  Suppress_v_cycle_output = true;
244 
245  // Enable the suppression of everything
246  Suppress_all_output = true;
247 
248  // Store the output stream pointer
249  Stream_pt = oomph_info.stream_pt();
250 
251  // Now set the oomph_info stream pointer to the null stream to
252  // disable all possible output
253  oomph_info.stream_pt() = &oomph_nullstream;
254  } // End of disable_output
oomph::OomphInfo::stream_pt
std::ostream *& stream_pt()
Access function for the stream pointer.
Definition: oomph_definitions.h:464
oomph::oomph_nullstream
Nullstream oomph_nullstream
Single (global) instantiation of the Nullstream.
Definition: oomph_definitions.cc:313

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time, oomph::oomph_info, oomph::oomph_nullstream, oomph::OomphInfo::stream_pt(), oomph::HelmholtzMGPreconditioner< DIM >::Stream_pt, oomph::HelmholtzMGPreconditioner< DIM >::Suppress_all_output, and oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ disable_smoother_and_superlu_doc_time()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::disable_smoother_and_superlu_doc_time ( )
inline

Suppress the output of both smoothers and SuperLU.

288  {
289  // Loop over all levels of the hierarchy
290  for (unsigned i = 0; i < Nlevel - 1; i++)
291  {
292  // Disable time documentation on each level (for each pre-smoother)
293  Pre_smoothers_storage_pt[i]->disable_doc_time();
294 
295  // Disable time documentation on each level (for each post-smoother)
296  Post_smoothers_storage_pt[i]->disable_doc_time();
297  }
298 
299  // We only need a direct solve on the coarsest level so this is the
300  // only place we need to silence SuperLU
301  Coarsest_matrix_mg_pt->linear_solver_pt()->disable_doc_time();
302  } // End of disable_smoother_and_superlu_doc_time
oomph::DoubleMatrixBase::linear_solver_pt
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: matrices.h:296
oomph::LinearSolver::disable_doc_time
void disable_doc_time()
Disable documentation of solve times.
Definition: linear_solver.h:116

References oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_matrix_mg_pt, oomph::LinearSolver::disable_doc_time(), i, oomph::DoubleMatrixBase::linear_solver_pt(), oomph::HelmholtzMGPreconditioner< DIM >::Nlevel, oomph::HelmholtzMGPreconditioner< DIM >::Post_smoothers_storage_pt, and oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoothers_storage_pt.

◆ disable_v_cycle_output()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::disable_v_cycle_output ( )
inline

Disable all output from mg_solve apart from the number of V-cycles used to solve the problem

227  {
228  // Set the value of Doc_time to false
229  Doc_time = false;
230 
231  // Enable the suppression of output from the V-cycle
232  Suppress_v_cycle_output = true;
233  } // End of disable_v_cycle_output

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time, and oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ enable_doc_time()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::enable_doc_time ( )
inline

Enable time documentation.

258  {
259  // Set the value of Doc_time to true
260  Doc_time = true;
261  } // End of enable_doc_time

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time.

◆ enable_output()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::enable_output ( )
inline

Enable the output from anything that could have been suppressed.

275  {
276  // Enable time documentation
277  Doc_time = true;
278 
279  // Enable output from everything during the full setup of the solver
280  Suppress_all_output = false;
281 
282  // Enable output from the MG solver
283  Suppress_v_cycle_output = false;
284  } // End of enable_output

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time, oomph::HelmholtzMGPreconditioner< DIM >::Suppress_all_output, and oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output.

◆ enable_v_cycle_output()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::enable_v_cycle_output ( )
inline

Enable the output of the V-cycle timings and other output.

265  {
266  // Enable time documentation
267  Doc_time = true;
268 
269  // Enable output from the MG solver
270  Suppress_v_cycle_output = false;
271  } // End of enable_v_cycle_output

References oomph::HelmholtzMGPreconditioner< DIM >::Doc_time, and oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output.

◆ full_setup()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::full_setup

Runs a full setup of the MG solver.

Do a full setup (assumes everything will be setup around the HelmholtzMGProblem pointer given in the constructor)

1104  {
1105 #ifdef OOMPH_HAS_MPI
1106  // Make sure that this is running in serial. Can't guarantee it'll
1107  // work when the problem is distributed over several processors
1108  if (MPI_Helpers::communicator_pt()->nproc() > 1)
1109  {
1110  // Throw a warning
1111  OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
1112  OOMPH_CURRENT_FUNCTION,
1113  OOMPH_EXCEPTION_LOCATION);
1114  }
1115 #endif
1116 
1117  // Initialise the timer start variable
1118  double t_fs_start = 0.0;
1119 
1120  // If we're allowed to output
1121  if (!Suppress_all_output)
1122  {
1123  // Start the timer
1124  t_fs_start = TimingHelpers::timer();
1125 
1126  // Notify user that the hierarchy of levels is complete
1127  oomph_info
1128  << "\n========Starting Setup of Multigrid Preconditioner========"
1129  << std::endl;
1130 
1131  // Notify user that the hierarchy of levels is complete
1132  oomph_info
1133  << "\nStarting the full setup of the Helmholtz multigrid solver."
1134  << std::endl;
1135  }
1136 
1137 #ifdef PARANOID
1138  // PARANOID check - Make sure the dimension of the solver matches the
1139  // dimension of the elements used in the problems mesh
1140  if (dynamic_cast<FiniteElement*>(Mg_problem_pt->mesh_pt()->element_pt(0))
1141  ->dim() != DIM)
1142  {
1143  // Create an error message
1144  std::string err_strng = "The dimension of the elements used in the mesh ";
1145  err_strng += "does not match the dimension of the solver.";
1146 
1147  // Throw the error
1148  throw OomphLibError(
1149  err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1150  }
1151 
1152  // PARANOID check - The elements of the bulk mesh must all be refineable
1153  // elements otherwise we cannot deal with this
1154  if (Mg_problem_pt->mg_bulk_mesh_pt() != 0)
1155  {
1156  // Find the number of elements in the bulk mesh
1157  unsigned n_elements = Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
1158 
1159  // Loop over the elements in the mesh and ensure that they are
1160  // all refineable elements
1161  for (unsigned el_counter = 0; el_counter < n_elements; el_counter++)
1162  {
1163  // Upcast global mesh element to a refineable element
1164  RefineableElement* el_pt = dynamic_cast<RefineableElement*>(
1165  Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));
1166 
1167  // Check if the upcast worked or not; if el_pt is a null pointer the
1168  // element is not refineable
1169  if (el_pt == 0)
1170  {
1171  // Create an output steam
1172  std::ostringstream error_message_stream;
1173 
1174  // Store the error message
1175  error_message_stream << "Element in bulk mesh could not be upcast to "
1176  << "a refineable element." << std::endl;
1177 
1178  // Throw an error
1179  throw OomphLibError(error_message_stream.str(),
1180  OOMPH_CURRENT_FUNCTION,
1181  OOMPH_EXCEPTION_LOCATION);
1182  }
1183  } // for (unsigned el_counter=0;el_counter<n_elements;el_counter++)
1184  }
1185  // If the provided bulk mesh pointer is a null pointer
1186  else
1187  {
1188  // Create an output steam
1189  std::ostringstream error_message_stream;
1190 
1191  // Store the error message
1192  error_message_stream
1193  << "The provided bulk mesh pointer is a null pointer. " << std::endl;
1194 
1195  // Throw an error
1196  throw OomphLibError(error_message_stream.str(),
1197  OOMPH_CURRENT_FUNCTION,
1198  OOMPH_EXCEPTION_LOCATION);
1199  } // if (Mg_problem_pt->mg_bulk_mesh_pt()!=0)
1200 #endif
1201 
1202  // If this is not the first Newton step then we will already have things
1203  // in storage. If this is the case, delete them
1204  clean_up_memory();
1205 
1206  // Resize the Mg_hierarchy_pt vector
1207  Mg_hierarchy_pt.resize(1, 0);
1208 
1209  // Set the pointer to the finest level as the first entry in Mg_hierarchy_pt
1210  Mg_hierarchy_pt[0] = Mg_problem_pt;
1211 
1212  // Create the hierarchy of levels
1213  setup_mg_hierarchy();
1214 
1215  // Set up the interpolation and restriction matrices
1216  setup_transfer_matrices();
1217 
1218  // Set up the data structures on each level, i.e. the system matrix,
1219  // LHS and RHS vectors and smoothers
1220  setup_mg_structures();
1221 
1222  // Set up the smoothers on all of the levels
1223  setup_smoothers();
1224 
1225  // Loop over all of the coarser levels
1226  for (unsigned i = 1; i < Nlevel; i++)
1227  {
1228  // Delete the i-th coarse-grid HelmholtzMGProblem object
1229  delete Mg_hierarchy_pt[i];
1230 
1231  // Set it to be a null pointer
1232  Mg_hierarchy_pt[i] = 0;
1233  }
1234 
1235  // Indicate that the full setup has been completed
1236  Has_been_setup = true;
1237 
1238  // If we're allowed to output to the screen
1239  if (!Suppress_all_output)
1240  {
1241  // Output the time taken to complete the full setup
1242  double t_fs_end = TimingHelpers::timer();
1243  double full_setup_time = t_fs_end - t_fs_start;
1244 
1245  // Output the CPU time
1246  oomph_info << "\nCPU time for full setup [sec]: " << full_setup_time
1247  << std::endl;
1248 
1249  // Notify user that the hierarchy of levels is complete
1250  oomph_info << "\n==============="
1251  << "Multigrid Full Setup Complete"
1252  << "==============\n"
1253  << std::endl;
1254  } // if (!Suppress_all_output)
1255  } // End of full_setup
oomph::HelmholtzMGPreconditioner::setup_mg_structures
void setup_mg_structures()
Set up the MG structures on each level.
Definition: helmholtz_geometric_multigrid.h:1476
oomph::HelmholtzMGPreconditioner::setup_transfer_matrices
void setup_transfer_matrices()
Setup the transfer matrices on each level.
Definition: helmholtz_geometric_multigrid.h:1422
oomph::HelmholtzMGPreconditioner::setup_smoothers
void setup_smoothers()
Set up the smoothers on all levels.
Definition: helmholtz_geometric_multigrid.h:2541
oomph::HelmholtzMGPreconditioner::setup_mg_hierarchy
void setup_mg_hierarchy()
Definition: helmholtz_geometric_multigrid.h:1263
oomph::HelmholtzMGProblem::mg_bulk_mesh_pt
virtual TreeBasedRefineableMeshBase * mg_bulk_mesh_pt()=0
oomph::MPI_Helpers::communicator_pt
static OomphCommunicator * communicator_pt()
access to the global oomph-lib communicator
Definition: oomph_utilities.cc:1046
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
DIM
#define DIM
Definition: linearised_navier_stokes_elements.h:44
oomph::Global_string_for_annotation::string
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286
oomph::TimingHelpers::timer
double timer()
returns the time in seconds after some point in past
Definition: oomph_utilities.cc:1295

References oomph::TerminateHelper::clean_up_memory(), oomph::MPI_Helpers::communicator_pt(), oomph::FiniteElement::dim(), DIM, i, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Global_string_for_annotation::string(), and oomph::TimingHelpers::timer().

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::setup().

◆ interpolate_and_correct()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::interpolate_and_correct ( const unsigned &  level)

Interpolate solution at current level onto next finer mesh and correct the solution x at that level

3597  {
3598 #ifdef PARANOID
3599  // Check to make sure we can actually restrict the vector
3600  if (!(level > 0))
3601  {
3602  // Throw an error
3603  throw OomphLibError("Input level exceeds the possible parameter choice.",
3604  OOMPH_CURRENT_FUNCTION,
3605  OOMPH_EXCEPTION_LOCATION);
3606  }
3607 #endif
3608 
3609  // Build distribution of a temporary vector (real part)
3610  DoubleVector temp_soln_r(
3611  X_mg_vectors_storage[level - 1][0].distribution_pt());
3612 
3613  // Build distribution of a temporary vector (imaginary part)
3614  DoubleVector temp_soln_c(
3615  X_mg_vectors_storage[level - 1][1].distribution_pt());
3616 
3617  // Interpolate the solution vector
3618  Interpolation_matrices_storage_pt[level - 1]->multiply(
3619  X_mg_vectors_storage[level][0], temp_soln_r);
3620 
3621  // Interpolate the solution vector
3622  Interpolation_matrices_storage_pt[level - 1]->multiply(
3623  X_mg_vectors_storage[level][1], temp_soln_c);
3624 
3625  // Update the real part of the solution
3626  X_mg_vectors_storage[level - 1][0] += temp_soln_r;
3627 
3628  // Update the imaginary part of the solution
3629  X_mg_vectors_storage[level - 1][1] += temp_soln_c;
3630 
3631  } // End of interpolate_and_correct
DoubleVector
std::vector< double > DoubleVector
loads clump configuration
Definition: ClumpInput.h:26
oomph::DistributableLinearAlgebraObject::distribution_pt
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
Definition: linear_algebra_distribution.h:457

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ interpolation_matrix_set() [1/2]

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::interpolation_matrix_set ( const unsigned &  level,
double *  value,
int *  col_index,
int *  row_st,
unsigned &  ncol,
unsigned &  nnz 
)
inline

Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be used as the full weighting restriction.

402  {
403  // Dynamically allocate the interpolation matrix
404  Interpolation_matrices_storage_pt[level] = new CRDoubleMatrix;
405 
406  // Build the matrix
407  Interpolation_matrices_storage_pt[level]->build_without_copy(
408  ncol, nnz, value, col_index, row_st);
409 
410  } // End of interpolation_matrix_set

References oomph::HelmholtzMGPreconditioner< DIM >::Interpolation_matrices_storage_pt, and Eigen::value.

◆ interpolation_matrix_set() [2/2]

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::interpolation_matrix_set ( const unsigned &  level,
Vector< double > &  value,
Vector< int > &  col_index,
Vector< int > &  row_st,
unsigned &  ncol,
unsigned &  nrow 
)
inline

Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be used as the full weighting restriction.

421  {
422  // Dynamically allocate the interpolation matrix
423  Interpolation_matrices_storage_pt[level] = new CRDoubleMatrix;
424 
425  // Make the distribution pointer
426  LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
427  Mg_hierarchy_pt[level]->communicator_pt(), nrow, false);
428 
429 #ifdef PARANOID
430 #ifdef OOMPH_HAS_MPI
431  // Set up the warning messages
432  std::string warning_message =
433  "Setup of interpolation matrix distribution ";
434  warning_message += "has not been tested with MPI.";
435 
436  // If we're not running the code in serial
437  if (dist_pt->communicator_pt()->nproc() > 1)
438  {
439  // Throw a warning
440  OomphLibWarning(
441  warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
442  }
443 #endif
444 #endif
445 
446  // Build the matrix itself
447  Interpolation_matrices_storage_pt[level]->build(
448  dist_pt, ncol, value, col_index, row_st);
449 
450  // Delete the newly created distribution pointer
451  delete dist_pt;
452 
453  // Make it a null pointer
454  dist_pt = 0;
455 
456  } // End of interpolation_matrix_set
oomph::DistributableLinearAlgebraObject::nrow
unsigned nrow() const
access function to the number of global rows.
Definition: linear_algebra_distribution.h:463

References oomph::LinearAlgebraDistribution::communicator_pt(), oomph::HelmholtzMGPreconditioner< DIM >::Interpolation_matrices_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::Mg_hierarchy_pt, oomph::OomphCommunicator::nproc(), oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Global_string_for_annotation::string(), and Eigen::value.

◆ iterations()

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::iterations ( ) const
inline

Number of iterations.

540  {
541  // Return the number of V-cycles which have been done
542  return V_cycle_counter;
543  } // End of iterations
oomph::HelmholtzMGPreconditioner::V_cycle_counter
unsigned V_cycle_counter
Pointer to counter for V-cycles.
Definition: helmholtz_geometric_multigrid.h:698

References oomph::HelmholtzMGPreconditioner< DIM >::V_cycle_counter.

◆ level_up_local_coord_of_node() [1/3]

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::level_up_local_coord_of_node ( const int &  son_type,
Vector< double > &  s 
)

Given the son_type of an element and a local node number j in that element with nnode_1d nodes per coordinate direction, return the local coordinate s in its father element. Needed in the setup of the interpolation matrices

◆ level_up_local_coord_of_node() [2/3]

void oomph::HelmholtzMGPreconditioner< 3 >::level_up_local_coord_of_node ( const int &  son_type,
Vector< double > &  s 
)

Given the son type of the element and the local coordinate s of a given node in the son element, return the local coordinate s in its father element. 3D case.

3443  {
3444  // If the element is unrefined between the levels the local coordinate
3445  // of the node in one element is the same as that in the other element
3446  // therefore we only need to perform calculations if the levels are
3447  // different (i.e. son_type is not OMEGA)
3448  if (son_type != Tree::OMEGA)
3449  {
3450  // Scale the local coordinate from the range [-1,1]x[-1,1]x[-1,1]
3451  // to the range [0,1]x[0,1]x[0,1] to match the width of the local
3452  // coordinate range of the fine element from the perspective of
3453  // the father element. This then simply requires a shift of the
3454  // coordinates to match which type of son element we're dealing with
3455  s[0] = (s[0] + 1.0) / 2.0;
3456  s[1] = (s[1] + 1.0) / 2.0;
3457  s[2] = (s[2] + 1.0) / 2.0;
3458 
3459  // Cases: The son_type determines how the local coordinates should be
3460  // shifted to give the local coordinates in the coarse mesh element
3461  switch (son_type)
3462  {
3463  case OcTreeNames::LDF:
3464  s[0] -= 1;
3465  s[1] -= 1;
3466  break;
3467 
3468  case OcTreeNames::LDB:
3469  s[0] -= 1;
3470  s[1] -= 1;
3471  s[2] -= 1;
3472  break;
3473 
3474  case OcTreeNames::LUF:
3475  s[0] -= 1;
3476  break;
3477 
3478  case OcTreeNames::LUB:
3479  s[0] -= 1;
3480  s[2] -= 1;
3481  break;
3482 
3483  case OcTreeNames::RDF:
3484  s[1] -= 1;
3485  break;
3486 
3487  case OcTreeNames::RDB:
3488  s[1] -= 1;
3489  s[2] -= 1;
3490  break;
3491 
3492  case OcTreeNames::RUF:
3493  break;
3494 
3495  case OcTreeNames::RUB:
3496  s[2] -= 1;
3497  break;
3498  }
3499  } // if son_type!=Tree::OMEGA
3500  } // End of level_up_local_coord_of_node
oomph::Tree::OMEGA
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
s
RealScalar s
Definition: level1_cplx_impl.h:130
oomph::OcTreeNames::RDF
@ RDF
Definition: octree.h:54
oomph::OcTreeNames::RUB
@ RUB
Definition: octree.h:52
oomph::OcTreeNames::LUF
@ LUF
Definition: octree.h:55
oomph::OcTreeNames::LDF
@ LDF
Definition: octree.h:53
oomph::OcTreeNames::RDB
@ RDB
Definition: octree.h:50
oomph::OcTreeNames::LUB
@ LUB
Definition: octree.h:51
oomph::OcTreeNames::RUF
@ RUF
Definition: octree.h:56
oomph::OcTreeNames::LDB
@ LDB
Definition: octree.h:49

References oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, oomph::Tree::OMEGA, oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, and s.

◆ level_up_local_coord_of_node() [3/3]

void oomph::HelmholtzMGPreconditioner< 2 >::level_up_local_coord_of_node ( const int &  son_type,
Vector< double > &  s 
)

Given the son type of the element and the local coordinate s of a given node in the son element, return the local coordinate s in its father element. 2D case.

3510  {
3511  // If the element is unrefined between the levels the local coordinate
3512  // of the node in one element is the same as that in the other element
3513  // therefore we only need to perform calculations if the levels are
3514  // different (i.e. son_type is not OMEGA)
3515  if (son_type != Tree::OMEGA)
3516  {
3517  // Scale the local coordinate from the range [-1,1]x[-1,1] to the range
3518  // [0,1]x[0,1] to match the width of the local coordinate range of the
3519  // fine element from the perspective of the father element. This
3520  // then simply requires a shift of the coordinates to match which type
3521  // of son element we're dealing with
3522  s[0] = (s[0] + 1.0) / 2.0;
3523  s[1] = (s[1] + 1.0) / 2.0;
3524 
3525  // Cases: The son_type determines how the local coordinates should be
3526  // shifted to give the local coordinates in the coarse mesh element
3527  switch (son_type)
3528  {
3529  // If we're dealing with the bottom-left element we need to shift
3530  // the range [0,1]x[0,1] to [-1,0]x[-1,0]
3531  case QuadTreeNames::SW:
3532  s[0] -= 1;
3533  s[1] -= 1;
3534  break;
3535 
3536  // If we're dealing with the bottom-right element we need to shift
3537  // the range [0,1]x[0,1] to [0,1]x[-1,0]
3538  case QuadTreeNames::SE:
3539  s[1] -= 1;
3540  break;
3541 
3542  // If we're dealing with the top-right element we need to shift the
3543  // range [0,1]x[0,1] to [0,1]x[0,1], i.e. nothing needs to be done
3544  case QuadTreeNames::NE:
3545  break;
3546 
3547  // If we're dealing with the top-left element we need to shift
3548  // the range [0,1]x[0,1] to [-1,0]x[0,1]
3549  case QuadTreeNames::NW:
3550  s[0] -= 1;
3551  break;
3552  }
3553  } // if son_type!=Tree::OMEGA
3554  } // End of level_up_local_coord_of_node
oomph::QuadTreeNames::SE
@ SE
Definition: quadtree.h:57
oomph::QuadTreeNames::NW
@ NW
Definition: quadtree.h:58
oomph::QuadTreeNames::NE
@ NE
Definition: quadtree.h:59
oomph::QuadTreeNames::SW
@ SW
Definition: quadtree.h:56

References oomph::QuadTreeNames::NE, oomph::QuadTreeNames::NW, oomph::Tree::OMEGA, s, oomph::QuadTreeNames::SE, and oomph::QuadTreeNames::SW.

◆ maximum_edge_width() [1/3]

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::maximum_edge_width ( const unsigned &  level,
double &  h 
)
private

Estimate the value of the parameter h on the level-th problem in the hierarchy.

◆ maximum_edge_width() [2/3]

void oomph::HelmholtzMGPreconditioner< 2 >::maximum_edge_width ( const unsigned &  level,
double &  h 
)
private

Find the value of the parameters h on the level-th problem in the hierarchy. The value of h is determined by looping over each element in the mesh and calculating the length of each side and take the maximum value.Note, this is a heuristic calculation; if the mesh is nonuniform or adaptive refinement is used then the value of h, is not the same everywhere so we find the maximum edge width instead. If, however, uniform refinement is used on a uniform mesh (using quad elements) then this will return the correct value of h.

This is the explicit template specialisation of the case DIM=2.

2224  {
2225  // Create a pointer to the "bulk" mesh
2226  TreeBasedRefineableMeshBase* bulk_mesh_pt =
2227  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt();
2228 
2229  // Reset the value of h
2230  h = 0.0;
2231 
2232  // Find out how many nodes there are along one edge of the first element.
2233  // We assume here that all elements have the same number of nodes
2234  unsigned nnode_1d =
2235  dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(0))->nnode_1d();
2236 
2237  // Sort out corner node IDs:
2238  //--------------------------
2239  // Initialise a vector to store local node IDs of the corners
2240  Vector<unsigned> corner_node_id(4, 0);
2241 
2242  // Identify the local node ID of the first corner
2243  corner_node_id[0] = 0;
2244 
2245  // Identify the local node ID of the second corner
2246  corner_node_id[1] = nnode_1d - 1;
2247 
2248  // Identify the local node ID of the third corner
2249  corner_node_id[2] = nnode_1d * nnode_1d - 1;
2250 
2251  // Identify the local node ID of the fourth corner
2252  corner_node_id[3] = nnode_1d * (nnode_1d - 1);
2253 
2254  // Create storage for the nodal information:
2255  //------------------------------------------
2256  // Pointer to the first corner node on the j-th edge
2257  Node* first_node_pt = 0;
2258 
2259  // Pointer to the second corner node on the j-th edge
2260  Node* second_node_pt = 0;
2261 
2262  // Vector to store the (Eulerian) position of the first corner node
2263  // along a given edge of the element
2264  Vector<double> first_node_x(2, 0.0);
2265 
2266  // Vector to store the (Eulerian) position of the second corner node
2267  // along a given edge of the element
2268  Vector<double> second_node_x(2, 0.0);
2269 
2270  // Calculate h:
2271  //-------------
2272  // Find out how many elements there are in the bulk mesh
2273  unsigned n_element = bulk_mesh_pt->nelement();
2274 
2275  // Store a pointer which will point to a given element in the bulk mesh
2276  FiniteElement* el_pt = 0;
2277 
2278  // Initialise a dummy variable to compare with h
2279  double test_h = 0.0;
2280 
2281  // Store the number of edges in a 2D quad element
2282  unsigned n_edge = 4;
2283 
2284  // Loop over all of the elements in the bulk mesh
2285  for (unsigned i = 0; i < n_element; i++)
2286  {
2287  // Upcast the pointer to the i-th element to a FiniteElement pointer
2288  el_pt = dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(i));
2289 
2290  // Loop over the edges of the element
2291  for (unsigned j = 0; j < n_edge; j++)
2292  {
2293  // Get the local node ID of the first corner node on the j-th edge
2294  first_node_pt = el_pt->node_pt(corner_node_id[j]);
2295 
2296  // Get the local node ID of the second corner node on the j-th edge
2297  second_node_pt = el_pt->node_pt(corner_node_id[(j + 1) % 4]);
2298 
2299  // Obtain the (Eulerian) position of the first corner node
2300  for (unsigned n = 0; n < 2; n++)
2301  {
2302  // Grab the n-th coordinate of this node
2303  first_node_x[n] = first_node_pt->x(n);
2304  }
2305 
2306  // Obtain the (Eulerian) position of the second corner node
2307  for (unsigned n = 0; n < 2; n++)
2308  {
2309  // Grab the n-th coordinate of this node
2310  second_node_x[n] = second_node_pt->x(n);
2311  }
2312 
2313  // Reset the value of test_h
2314  test_h = 0.0;
2315 
2316  // Calculate the norm of (second_node_x-first_node_x)
2317  for (unsigned n = 0; n < 2; n++)
2318  {
2319  // Update the value of test_h
2320  test_h += pow(second_node_x[n] - first_node_x[n], 2.0);
2321  }
2322 
2323  // Square root the value of h
2324  test_h = sqrt(test_h);
2325 
2326  // Check if the value of test_h is greater than h
2327  if (test_h > h)
2328  {
2329  // If the value of test_h is greater than h then store it
2330  h = test_h;
2331  }
2332  } // for (unsigned j=0;j<n_edge;j++)
2333  } // for (unsigned i=0;i<n_element;i++)
2334  } // End of maximum_edge_width
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625

References oomph::Mesh::element_pt(), i, j, n, oomph::Mesh::nelement(), oomph::FiniteElement::node_pt(), Eigen::bfloat16_impl::pow(), sqrt(), and oomph::Node::x().

◆ maximum_edge_width() [3/3]

void oomph::HelmholtzMGPreconditioner< 3 >::maximum_edge_width ( const unsigned &  level,
double &  h 
)
private

Find the value of the parameters h on the level-th problem in the hierarchy. The value of h is determined by looping over each element in the mesh and calculating the length of each side and take the maximum value. Note, this is a heuristic calculation; if the mesh is non-uniform or adaptive refinement is used then the value of h, is not the same everywhere so we find the maximum edge width instead. If, however, uniform refinement is used on a uniform mesh (using quad elements) then this will return the correct value of h.

This is the explicit template specialisation of the case DIM=3. The calculation of h is different here. In 2D we were able to loop over each pair of nodes in an anti-clockwise manner since the only node pairs were {(C0,C1),(C1,C2),(C2,C3),(C3,C0)} where CN denotes the N-th corner in the element. In 3D this method cannot be used since we have 12 edges to consider.

2356  {
2357  // Create a pointer to the "bulk" mesh
2358  TreeBasedRefineableMeshBase* bulk_mesh_pt =
2359  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt();
2360 
2361  //--------------------------------
2362  // Find the maximum edge width, h:
2363  //--------------------------------
2364  // Reset the value of h
2365  h = 0.0;
2366 
2367  // Find out how many nodes there are along one edge of the first element.
2368  // We assume here that all elements have the same number of nodes
2369  unsigned nnode_1d =
2370  dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(0))->nnode_1d();
2371 
2372  // Sort out corner node IDs:
2373  //--------------------------
2374  // Grab the octree pointer from the first element in the bulk mesh
2375  OcTree* octree_pt =
2376  dynamic_cast<RefineableQElement<3>*>(bulk_mesh_pt->element_pt(0))
2377  ->octree_pt();
2378 
2379  // Initialise a vector to store local node IDs of the corners
2380  Vector<unsigned> corner_node_id(8, 0);
2381 
2382  // Identify the local node ID of the first corner
2383  corner_node_id[0] =
2384  octree_pt->vertex_to_node_number(OcTreeNames::LDB, nnode_1d);
2385 
2386  // Identify the local node ID of the second corner
2387  corner_node_id[1] =
2388  octree_pt->vertex_to_node_number(OcTreeNames::RDB, nnode_1d);
2389 
2390  // Identify the local node ID of the third corner
2391  corner_node_id[2] =
2392  octree_pt->vertex_to_node_number(OcTreeNames::LUB, nnode_1d);
2393 
2394  // Identify the local node ID of the fourth corner
2395  corner_node_id[3] =
2396  octree_pt->vertex_to_node_number(OcTreeNames::RUB, nnode_1d);
2397 
2398  // Identify the local node ID of the fifth corner
2399  corner_node_id[4] =
2400  octree_pt->vertex_to_node_number(OcTreeNames::LDF, nnode_1d);
2401 
2402  // Identify the local node ID of the sixth corner
2403  corner_node_id[5] =
2404  octree_pt->vertex_to_node_number(OcTreeNames::RDF, nnode_1d);
2405 
2406  // Identify the local node ID of the seventh corner
2407  corner_node_id[6] =
2408  octree_pt->vertex_to_node_number(OcTreeNames::LUF, nnode_1d);
2409 
2410  // Identify the local node ID of the eighth corner
2411  corner_node_id[7] =
2412  octree_pt->vertex_to_node_number(OcTreeNames::RUF, nnode_1d);
2413 
2414  // Sort out the local node ID pairs:
2415  //----------------------------------
2416  // Store the number of edges in a 3D cubic element
2417  unsigned n_edge = 12;
2418 
2419  // Create a vector to store the pairs of adjacent nodes
2420  Vector<std::pair<unsigned, unsigned>> edge_node_pair(n_edge);
2421 
2422  // First edge
2423  edge_node_pair[0] = std::make_pair(corner_node_id[0], corner_node_id[1]);
2424 
2425  // Second edge
2426  edge_node_pair[1] = std::make_pair(corner_node_id[0], corner_node_id[2]);
2427 
2428  // Third edge
2429  edge_node_pair[2] = std::make_pair(corner_node_id[0], corner_node_id[4]);
2430 
2431  // Fourth edge
2432  edge_node_pair[3] = std::make_pair(corner_node_id[1], corner_node_id[3]);
2433 
2434  // Fifth edge
2435  edge_node_pair[4] = std::make_pair(corner_node_id[1], corner_node_id[5]);
2436 
2437  // Sixth edge
2438  edge_node_pair[5] = std::make_pair(corner_node_id[2], corner_node_id[3]);
2439 
2440  // Seventh edge
2441  edge_node_pair[6] = std::make_pair(corner_node_id[2], corner_node_id[6]);
2442 
2443  // Eighth edge
2444  edge_node_pair[7] = std::make_pair(corner_node_id[3], corner_node_id[7]);
2445 
2446  // Ninth edge
2447  edge_node_pair[8] = std::make_pair(corner_node_id[4], corner_node_id[5]);
2448 
2449  // Tenth edge
2450  edge_node_pair[9] = std::make_pair(corner_node_id[4], corner_node_id[6]);
2451 
2452  // Eleventh edge
2453  edge_node_pair[10] = std::make_pair(corner_node_id[5], corner_node_id[7]);
2454 
2455  // Twelfth edge
2456  edge_node_pair[11] = std::make_pair(corner_node_id[6], corner_node_id[7]);
2457 
2458  // Create storage for the nodal information:
2459  //------------------------------------------
2460  // Pointer to the first corner node on the j-th edge
2461  Node* first_node_pt = 0;
2462 
2463  // Pointer to the second corner node on the j-th edge
2464  Node* second_node_pt = 0;
2465 
2466  // Vector to store the (Eulerian) position of the first corner node
2467  // along a given edge of the element
2468  Vector<double> first_node_x(3, 0.0);
2469 
2470  // Vector to store the (Eulerian) position of the second corner node
2471  // along a given edge of the element
2472  Vector<double> second_node_x(3, 0.0);
2473 
2474  // Calculate h:
2475  //-------------
2476  // Find out how many elements there are in the bulk mesh
2477  unsigned n_element = bulk_mesh_pt->nelement();
2478 
2479  // Store a pointer which will point to a given element in the bulk mesh
2480  FiniteElement* el_pt = 0;
2481 
2482  // Initialise a dummy variable to compare with h
2483  double test_h = 0.0;
2484 
2485  // Loop over all of the elements in the bulk mesh
2486  for (unsigned i = 0; i < n_element; i++)
2487  {
2488  // Upcast the pointer to the i-th element to a FiniteElement pointer
2489  el_pt = dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(i));
2490 
2491  // Loop over the edges of the element
2492  for (unsigned j = 0; j < n_edge; j++)
2493  {
2494  // Get the local node ID of the first corner node on the j-th edge
2495  first_node_pt = el_pt->node_pt(edge_node_pair[j].first);
2496 
2497  // Get the local node ID of the second corner node on the j-th edge
2498  second_node_pt = el_pt->node_pt(edge_node_pair[j].second);
2499 
2500  // Obtain the (Eulerian) position of the first corner node
2501  for (unsigned n = 0; n < 3; n++)
2502  {
2503  // Grab the n-th coordinate of this node
2504  first_node_x[n] = first_node_pt->x(n);
2505  }
2506 
2507  // Obtain the (Eulerian) position of the second corner node
2508  for (unsigned n = 0; n < 3; n++)
2509  {
2510  // Grab the n-th coordinate of this node
2511  second_node_x[n] = second_node_pt->x(n);
2512  }
2513 
2514  // Reset the value of test_h
2515  test_h = 0.0;
2516 
2517  // Calculate the norm of (second_node_x-first_node_x)
2518  for (unsigned n = 0; n < 3; n++)
2519  {
2520  // Update the value of test_h
2521  test_h += pow(second_node_x[n] - first_node_x[n], 2.0);
2522  }
2523 
2524  // Square root the value of h
2525  test_h = sqrt(test_h);
2526 
2527  // Check if the value of test_h is greater than h
2528  if (test_h > h)
2529  {
2530  // If the value of test_h is greater than h then store it
2531  h = test_h;
2532  }
2533  } // for (unsigned j=0;j<n_edge;j++)
2534  } // for (unsigned i=0;i<n_element;i++)
2535  } // End of maximum_edge_width

References oomph::Mesh::element_pt(), i, j, oomph::OcTreeNames::LDB, oomph::OcTreeNames::LDF, oomph::OcTreeNames::LUB, oomph::OcTreeNames::LUF, n, oomph::Mesh::nelement(), oomph::FiniteElement::node_pt(), Eigen::bfloat16_impl::pow(), oomph::OcTreeNames::RDB, oomph::OcTreeNames::RDF, oomph::OcTreeNames::RUB, oomph::OcTreeNames::RUF, sqrt(), oomph::OcTree::vertex_to_node_number(), and oomph::Node::x().

◆ mg_solve()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::mg_solve ( Vector< DoubleVector > &  result)
private

Do the actual solve – this is called through the pure virtual solve function in the LinearSolver base class. The function is stored as protected to allow the HelmholtzMGPreconditioner derived class to use the solver

Linear solver. This is where the general V-cycle algorithm is implemented

3639  {
3640  // If we're allowed to time things
3641  double t_mg_start = 0.0;
3642  if (!Suppress_v_cycle_output)
3643  {
3644  // Start the clock!
3645  t_mg_start = TimingHelpers::timer();
3646  }
3647 
3648  // Current level
3649  unsigned finest_level = 0;
3650 
3651  // Initialise the V-cycle counter
3652  V_cycle_counter = 0;
3653 
3654  // Calculate the norm of the residual then output
3655  double normalised_residual_norm = residual_norm(finest_level);
3656  if (!Suppress_v_cycle_output)
3657  {
3658  oomph_info << "\nResidual on finest level for V-cycle: "
3659  << normalised_residual_norm << std::endl;
3660  }
3661 
3662  // Outer loop over V-cycles
3663  //-------------------------
3664  // While the tolerance is not met and the maximum number of
3665  // V-cycles has not been completed
3666  while ((normalised_residual_norm > Tolerance) &&
3667  (V_cycle_counter != Nvcycle))
3668  {
3669  // If the user did not wish to suppress the V-cycle output
3670  if (!Suppress_v_cycle_output)
3671  {
3672  // Output the V-cycle counter
3673  oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
3674  }
3675 
3676  //---------------------------------------------------------------------
3677  // Loop downwards over all levels that have coarser levels beneath them
3678  //---------------------------------------------------------------------
3679  for (unsigned i = 0; i < Nlevel - 1; i++)
3680  {
3681  // Initialise X_mg and Residual_mg to 0.0 except for the finest level
3682  // since X_mg contains the current approximation to the solution and
3683  // Residual_mg contains the RHS vector on the finest level
3684  if (i != 0)
3685  {
3686  // Initialise the real part of the solution vector
3687  X_mg_vectors_storage[i][0].initialise(0.0);
3688 
3689  // Initialise the imaginary part of the solution vector
3690  X_mg_vectors_storage[i][1].initialise(0.0);
3691 
3692  // Initialise the real part of the residual vector
3693  Residual_mg_vectors_storage[i][0].initialise(0.0);
3694 
3695  // Initialise the imaginary part of the residual vector
3696  Residual_mg_vectors_storage[i][1].initialise(0.0);
3697  }
3698 
3699  // Perform a few pre-smoothing steps and return vector that contains
3700  // the residuals of the linear system at this level.
3701  pre_smooth(i);
3702 
3703  // Restrict the residual to the next coarser mesh and
3704  // assign it to the RHS vector at that level
3705  restrict_residual(i);
3706 
3707  } // Moving down the V-cycle
3708 
3709  //-----------------------------------------------------------
3710  // Reached the lowest level: Do a direct solve, using the RHS
3711  // vector obtained by restriction from above.
3712  //-----------------------------------------------------------
3713  direct_solve();
3714 
3715  //---------------------------------------------------------------
3716  // Loop upwards over all levels that have finer levels above them
3717  //---------------------------------------------------------------
3718  for (unsigned i = Nlevel - 1; i > 0; i--)
3719  {
3720  // Interpolate solution at current level onto
3721  // next finer mesh and correct the solution x at that level
3722  interpolate_and_correct(i);
3723 
3724  // Perform a few post-smoothing steps (ignore
3725  // vector that contains the residuals of the linear system
3726  // at this level)
3727  post_smooth(i - 1);
3728  }
3729 
3730  // Set counter for number of cycles (for doc)
3731  V_cycle_counter++;
3732 
3733  // Calculate the new residual norm then output (if allowed)
3734  normalised_residual_norm = residual_norm(finest_level);
3735 
3736  // Print the residual on the finest level
3737  if (!Suppress_v_cycle_output)
3738  {
3739  oomph_info << "Residual on finest level of V-cycle: "
3740  << normalised_residual_norm << std::endl;
3741  }
3742  } // End of the V-cycles
3743 
3744  // Copy the solution into the result vector
3745  result = X_mg_vectors_storage[finest_level];
3746 
3747  // Need an extra line space if V-cycle output is suppressed
3748  if (!Suppress_v_cycle_output)
3749  {
3750  oomph_info << std::endl;
3751  }
3752 
3753  // If all output is to be suppressed
3754  if (!Suppress_all_output)
3755  {
3756  // Output number of V-cycles taken to solve
3757  if (normalised_residual_norm < Tolerance)
3758  {
3759  Has_been_solved = true;
3760  }
3761  } // if (!Suppress_all_output)
3762 
3763  // If the V-cycle output isn't suppressed
3764  if (!Suppress_v_cycle_output)
3765  {
3766  // Stop the clock
3767  double t_mg_end = TimingHelpers::timer();
3768  double total_G_setup_time = double(t_mg_end - t_mg_start);
3769  oomph_info << "CPU time for MG solver [sec]: " << total_G_setup_time
3770  << std::endl;
3771  }
3772  } // end of mg_solve
oomph::HelmholtzMGPreconditioner::restrict_residual
void restrict_residual(const unsigned &level)
Definition: helmholtz_geometric_multigrid.h:3563
oomph::HelmholtzMGPreconditioner::pre_smooth
void pre_smooth(const unsigned &level)
Definition: helmholtz_geometric_multigrid.h:325
oomph::HelmholtzMGPreconditioner::direct_solve
void direct_solve()
Call the direct solver (SuperLU) to solve the problem exactly.
Definition: helmholtz_geometric_multigrid.h:375
oomph::HelmholtzMGPreconditioner::post_smooth
void post_smooth(const unsigned &level)
Definition: helmholtz_geometric_multigrid.h:340
oomph::HelmholtzMGPreconditioner::interpolate_and_correct
void interpolate_and_correct(const unsigned &level)
Definition: helmholtz_geometric_multigrid.h:3595
oomph::HelmholtzMGPreconditioner::Residual_mg_vectors_storage
Vector< Vector< DoubleVector > > Residual_mg_vectors_storage
Definition: helmholtz_geometric_multigrid.h:666
oomph::HelmholtzMGPreconditioner::residual_norm
double residual_norm(const unsigned &level)
Definition: helmholtz_geometric_multigrid.h:352

References i, oomph::oomph_info, and oomph::TimingHelpers::timer().

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve().

◆ npost_smooth()

template<unsigned DIM>
unsigned& oomph::HelmholtzMGPreconditioner< DIM >::npost_smooth ( )
inline

Return the number of post-smoothing iterations (lvalue)

306  {
307  // Return the number of post-smoothing iterations to be done on each
308  // level of the hierarchy
309  return Npost_smooth;
310  } // End of npost_smooth

References oomph::HelmholtzMGPreconditioner< DIM >::Npost_smooth.

◆ npre_smooth()

template<unsigned DIM>
unsigned& oomph::HelmholtzMGPreconditioner< DIM >::npre_smooth ( )
inline

Return the number of pre-smoothing iterations (lvalue)

314  {
315  // Return the number of pre-smoothing iterations to be done on each
316  // level of the hierarchy
317  return Npre_smooth;
318  } // End of npre_smooth

References oomph::HelmholtzMGPreconditioner< DIM >::Npre_smooth.

◆ post_smooth()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::post_smooth ( const unsigned &  level)
inline

Post-smoother: Perform max_iter smoothing steps on the linear system Ax=b with current RHS vector, b, starting with current solution vector, x. Uses the default smoother (set in the HelmholtzMGProblem constructor) which can be overloaded for specific problem.

341  {
342  // Run post-smoother 'max_iter' times
343  Post_smoothers_storage_pt[level]->complex_smoother_solve(
344  Rhs_mg_vectors_storage[level], X_mg_vectors_storage[level]);
345 
346  // Calculate the residual vector on this level
347  residual_norm(level);
348  } // End of post_smooth

References oomph::HelmholtzMGPreconditioner< DIM >::Post_smoothers_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::residual_norm(), oomph::HelmholtzMGPreconditioner< DIM >::Rhs_mg_vectors_storage, and oomph::HelmholtzMGPreconditioner< DIM >::X_mg_vectors_storage.

◆ pre_smooth()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::pre_smooth ( const unsigned &  level)
inline

Pre-smoother: Perform 'max_iter' smoothing steps on the linear system Ax=b with current RHS vector, b, starting with current solution vector, x. Return the residual vector r=b-Ax. Uses the default smoother (set in the HelmholtzMGProblem constructor) which can be overloaded for a specific problem.

326  {
327  // Run pre-smoother 'max_iter' times
328  Pre_smoothers_storage_pt[level]->complex_smoother_solve(
329  Rhs_mg_vectors_storage[level], X_mg_vectors_storage[level]);
330 
331  // Calculate the residual vector on this level
332  residual_norm(level);
333  } // End of pre_smooth

References oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoothers_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::residual_norm(), oomph::HelmholtzMGPreconditioner< DIM >::Rhs_mg_vectors_storage, and oomph::HelmholtzMGPreconditioner< DIM >::X_mg_vectors_storage.

◆ preconditioner_solve()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve ( const DoubleVector &  r,
DoubleVector &  z 
)
inlinevirtual

Function applies MG to the vector r for a full solve.

Implements oomph::Preconditioner.

505  {
506  // Split up the RHS vector into DoubleVectors, whose entries are
507  // arranged to match the matrix blocks and assign it
508  this->get_block_vectors(r, Rhs_mg_vectors_storage[0]);
509 
510  // Split up the solution vector into DoubleVectors, whose entries are
511  // arranged to match the matrix blocks and assign it
512  this->get_block_vectors(z, X_mg_vectors_storage[0]);
513 
514  // Run the MG method and assign the solution to z
515  this->mg_solve(X_mg_vectors_storage[0]);
516 
517  // Copy solution in block vectors block_z back to z
518  this->return_block_vectors(X_mg_vectors_storage[0], z);
519 
520  // Only output if the V-cycle output isn't suppressed
521  if (!(this->Suppress_v_cycle_output))
522  {
523  // Notify user that the hierarchy of levels is complete
524  oomph_info << "\n=========="
525  << "Multigrid Preconditioner Solve Complete"
526  << "=========" << std::endl;
527  }
528 
529  // Only enable and assign the stream pointer again if we originally
530  // suppressed everything otherwise it won't be set yet
531  if (this->Suppress_all_output)
532  {
533  // Now enable the stream pointer again
534  oomph_info.stream_pt() = this->Stream_pt;
535  }
536  } // End of preconditioner_solve
oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vectors
void return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
Definition: block_preconditioner.cc:3443
oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vectors
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Definition: block_preconditioner.cc:2939
oomph::HelmholtzMGPreconditioner::mg_solve
void mg_solve(Vector< DoubleVector > &result)
Definition: helmholtz_geometric_multigrid.h:3638
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20

References oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vectors(), oomph::HelmholtzMGPreconditioner< DIM >::mg_solve(), oomph::oomph_info, oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vectors(), oomph::HelmholtzMGPreconditioner< DIM >::Rhs_mg_vectors_storage, oomph::OomphInfo::stream_pt(), oomph::HelmholtzMGPreconditioner< DIM >::Stream_pt, oomph::HelmholtzMGPreconditioner< DIM >::Suppress_all_output, oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output, and oomph::HelmholtzMGPreconditioner< DIM >::X_mg_vectors_storage.

◆ residual_norm() [1/2]

template<unsigned DIM>
double oomph::HelmholtzMGPreconditioner< DIM >::residual_norm ( const unsigned &  level)
inline

Return norm of residual r=b-Ax and the residual vector itself on the level-th level

353  {
354  // Loop over the real and imaginary part of the residual vector on
355  // the given level
356  for (unsigned j = 0; j < 2; j++)
357  {
358  // And zero the entries of residual
359  Residual_mg_vectors_storage[level][j].initialise(0.0);
360  }
361 
362  // Return the norm of the residual
363  return residual_norm(level, Residual_mg_vectors_storage[level]);
364  } // End of residual_norm

References j, and oomph::HelmholtzMGPreconditioner< DIM >::Residual_mg_vectors_storage.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::post_smooth(), and oomph::HelmholtzMGPreconditioner< DIM >::pre_smooth().

◆ residual_norm() [2/2]

template<unsigned DIM>
double oomph::HelmholtzMGPreconditioner< DIM >::residual_norm ( const unsigned &  level,
Vector< DoubleVector > &  residual 
)

Calculate the norm of the residual vector, r=b-Ax.

Calculating the residual r=b-Ax in the complex case requires more care than the real case. To calculate the residual vector we split A, x and b into their complex components: r = b - A*x, = (b_r + i*b_c) - (A_r + i*A_c)*(x_r + i*x_c), = [b_r - A_r*x_r + A_c*x_c] + i*[b_c - A_r*x_c - A_c*x_r], ==> real(r) = b_r - A_r*x_r + A_c*x_c, & imag(r) = b_c - A_r*x_c - A_c*x_r.

736  {
737  // Number of rows in each block vector
738  unsigned n_rows = X_mg_vectors_storage[level][0].nrow();
739 
740 #ifdef PARANOID
741  // PARANOID check - if the residual vector doesn't have length 2 it cannot
742  // be used here since we need two vectors corresponding to the real and
743  // imaginary part of the residual
744  if (residual.size() != 2)
745  {
746  // Throw an error if the residual vector doesn't have the correct length
747  throw OomphLibError("This residual vector must have length 2. ",
748  OOMPH_CURRENT_FUNCTION,
749  OOMPH_EXCEPTION_LOCATION);
750  }
751  if (residual[0].nrow() != residual[1].nrow())
752  {
753  // Create an output stream
754  std::ostringstream error_message_stream;
755 
756  // Store the error message
757  error_message_stream << "Residual[0] has length: " << residual[0].nrow()
758  << "\nResidual[1] has length: " << residual[1].nrow()
759  << "\nThis method requires that the constituent "
760  << "DoubleVectors in residual have the same length. "
761  << std::endl;
762 
763  // Throw an error
764  throw OomphLibError(error_message_stream.str(),
765  OOMPH_CURRENT_FUNCTION,
766  OOMPH_EXCEPTION_LOCATION);
767  }
768 #endif
769 
770  // Loop over the block vectors
771  for (unsigned i = 0; i < 2; i++)
772  {
773  // Start by setting the distribution of the residuals vector if
774  // it is not set up
775  if (!residual[i].distribution_built())
776  {
777  // Set up distribution
778  LinearAlgebraDistribution dist(
779  Mg_hierarchy_pt[level]->communicator_pt(), n_rows, false);
780 
781  // Build the distribution
782  residual[i].build(&dist, 0.0);
783  }
784  // Otherwise just zero the entries of residual
785  else
786  {
787 #ifdef PARANOID
788  // PARANOID check - if the residuals are distributed then this method
789  // cannot be used, a distributed residuals can only be assembled by
790  // get_jacobian(...) for CRDoubleMatrices
791  if (residual[i].distributed())
792  {
793  throw OomphLibError(
794  "This method can only assemble a non-distributed residuals vector ",
795  OOMPH_CURRENT_FUNCTION,
796  OOMPH_EXCEPTION_LOCATION);
797  }
798 #endif
799 
800  // And zero the entries of residual
801  residual[i].initialise(0.0);
802  }
803  }
804 
805  // Store the pointer to the distribution of Matrix_real_pt (the same as
806  // Matrix_imag_pt presumably so we only need to use one)
807  LinearAlgebraDistribution* dist_pt =
808  Mg_matrices_storage_pt[level][0]->distribution_pt();
809 
810  // Create 4 temporary vectors to store the various matrix-vector products
811  // required. The appropriate combination has been appended to the name of
812  // the vector:
813  // Vector to store A_r*x_r:
814  DoubleVector temp_vec_rr(dist_pt, 0.0);
815 
816  // Vector to store A_r*x_c:
817  DoubleVector temp_vec_rc(dist_pt, 0.0);
818 
819  // Vector to store A_c*x_r:
820  DoubleVector temp_vec_cr(dist_pt, 0.0);
821 
822  // Vector to store A_c*x_c:
823  DoubleVector temp_vec_cc(dist_pt, 0.0);
824 
825  // We can't delete the distribution pointer because the Jacobian on the
826  // finest level is using it but we can make it a null pointer
827  dist_pt = 0;
828 
829  // Calculate the different combinations of A*x (or A*e depending on the
830  // level in the heirarchy) in the complex case:
831  // A_r*x_r:
832  Mg_matrices_storage_pt[level][0]->multiply(X_mg_vectors_storage[level][0],
833  temp_vec_rr);
834 
835  // A_r*x_c:
836  Mg_matrices_storage_pt[level][0]->multiply(X_mg_vectors_storage[level][1],
837  temp_vec_rc);
838 
839  // A_c*x_r:
840  Mg_matrices_storage_pt[level][1]->multiply(X_mg_vectors_storage[level][0],
841  temp_vec_cr);
842 
843  // A_c*x_c:
844  Mg_matrices_storage_pt[level][1]->multiply(X_mg_vectors_storage[level][1],
845  temp_vec_cc);
846 
847  // Real part of the residual:
848  residual[0] = Rhs_mg_vectors_storage[level][0];
849  residual[0] -= temp_vec_rr;
850  residual[0] += temp_vec_cc;
851 
852  // Imaginary part of the residual:
853  residual[1] = Rhs_mg_vectors_storage[level][1];
854  residual[1] -= temp_vec_rc;
855  residual[1] -= temp_vec_cr;
856 
857  // Get the residual norm of the real part of the residual vector
858  double norm_r = residual[0].norm();
859 
860  // Get the residual norm of the imaginary part of the residual vector
861  double norm_c = residual[1].norm();
862 
863  // Return the true norm of the residual vector which depends on the
864  // norm of the real part and the norm of the imaginary part
865  return sqrt(norm_r * norm_r + norm_c * norm_c);
866  }
oomph::DistributableLinearAlgebraObject::distributed
bool distributed() const
distribution is serial or distributed
Definition: linear_algebra_distribution.h:493
oomph::DistributableLinearAlgebraObject::distribution_built
bool distribution_built() const
Definition: linear_algebra_distribution.h:500

References i, oomph::Vector< _Tp >::initialise(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and sqrt().

◆ restrict_residual()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::restrict_residual ( const unsigned &  level)

Restrict residual (computed on level-th MG level) to the next coarser mesh and stick it into the coarse mesh RHS vector.

Restrict residual (computed on current MG level) to next coarser mesh and stick it into the coarse mesh RHS vector using the restriction matrix (if restrict_flag=1) or the transpose of the interpolation matrix (if restrict_flag=2)

3564  {
3565 #ifdef PARANOID
3566  // Check to make sure we can actually restrict the vector
3567  if (!(level < Nlevel - 1))
3568  {
3569  // Throw an error
3570  throw OomphLibError("Input level exceeds the possible parameter choice.",
3571  OOMPH_CURRENT_FUNCTION,
3572  OOMPH_EXCEPTION_LOCATION);
3573  }
3574 #endif
3575 
3576  // Multiply the real part of the residual vector by the restriction
3577  // matrix on the level-th level
3578  Restriction_matrices_storage_pt[level]->multiply(
3579  Residual_mg_vectors_storage[level][0],
3580  Rhs_mg_vectors_storage[level + 1][0]);
3581 
3582  // Multiply the imaginary part of the residual vector by the restriction
3583  // matrix on the level-th level
3584  Restriction_matrices_storage_pt[level]->multiply(
3585  Residual_mg_vectors_storage[level][1],
3586  Rhs_mg_vectors_storage[level + 1][1]);
3587 
3588  } // End of restrict_residual

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ set_post_smoother_factory_function()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::set_post_smoother_factory_function ( PostSmootherFactoryFctPt  post_smoother_fn)
inline

Access function to set the post-smoother creation function.

104  {
105  // Assign the function pointer
106  Post_smoother_factory_function_pt = post_smoother_fn;
107  }

References oomph::HelmholtzMGPreconditioner< DIM >::Post_smoother_factory_function_pt.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ set_pre_smoother_factory_function()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::set_pre_smoother_factory_function ( PreSmootherFactoryFctPt  pre_smoother_fn)
inline

Access function to set the pre-smoother creation function.

96  {
97  // Assign the function pointer
98  Pre_smoother_factory_function_pt = pre_smoother_fn;
99  }

References oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoother_factory_function_pt.

Referenced by PMLStructuredCubicHelmholtz< ELEMENT >::set_gmres_multigrid_solver(), and PMLHelmholtzMGProblem< ELEMENT >::set_gmres_multigrid_solver().

◆ set_restriction_matrices_as_interpolation_transposes()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::set_restriction_matrices_as_interpolation_transposes ( )
inline

Builds a CRDoubleMatrix on each level that is used to restrict the residual between levels. The transpose can be used as the interpolation matrix

462  {
463  for (unsigned i = 0; i < Nlevel - 1; i++)
464  {
465  // Dynamically allocate the restriction matrix
466  Restriction_matrices_storage_pt[i] = new CRDoubleMatrix;
467 
468  // Set the restriction matrix to be the transpose of the
469  // interpolation matrix
470  Interpolation_matrices_storage_pt[i]->get_matrix_transpose(
471  Restriction_matrices_storage_pt[i]);
472  }
473  } // End of set_restriction_matrices_as_interpolation_transposes

References i, oomph::HelmholtzMGPreconditioner< DIM >::Interpolation_matrices_storage_pt, oomph::HelmholtzMGPreconditioner< DIM >::Nlevel, and oomph::HelmholtzMGPreconditioner< DIM >::Restriction_matrices_storage_pt.

◆ setup() [1/4]

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup ( )
inlineprivatevirtual

Function to set up the hierachy of levels. Creates a vector of pointers to each MG level

Implements oomph::Preconditioner.

574  {
575  // Run the full setup
576  this->full_setup();
577 
578  // Only enable and assign the stream pointer again if we originally
579  // suppressed everything otherwise it won't be set yet
580  if (this->Suppress_all_output)
581  {
582  // Now enable the stream pointer again
583  oomph_info.stream_pt() = this->Stream_pt;
584  }
585  } // End of setup
oomph::HelmholtzMGPreconditioner::full_setup
void full_setup()
Runs a full setup of the MG solver.
Definition: helmholtz_geometric_multigrid.h:1103

References oomph::HelmholtzMGPreconditioner< DIM >::full_setup(), oomph::oomph_info, oomph::OomphInfo::stream_pt(), oomph::HelmholtzMGPreconditioner< DIM >::Stream_pt, and oomph::HelmholtzMGPreconditioner< DIM >::Suppress_all_output.

◆ setup() [2/4]

template<unsigned DIM>
virtual void oomph::Preconditioner::setup
virtual

Use the version in the Preconditioner base class for the alternative setup function that takes a matrix pointer as an argument.

Implements oomph::Preconditioner.

◆ setup() [3/4]

template<unsigned DIM>
void oomph::Preconditioner::setup
inlinevirtual

Use the version in the Preconditioner base class for the alternative setup function that takes a matrix pointer as an argument.

Implements oomph::Preconditioner.

121  {
122  ObsoleteCode::obsolete();
123  setup(matrix_pt);
124  }
oomph::BlockPreconditioner< CRDoubleMatrix >::matrix_pt
CRDoubleMatrix * matrix_pt() const
Definition: block_preconditioner.h:520
oomph::HelmholtzMGPreconditioner::setup
void setup()
Definition: helmholtz_geometric_multigrid.h:573
oomph::ObsoleteCode::obsolete
void obsolete()
Output warning message.
Definition: oomph_utilities.cc:1102

◆ setup() [4/4]

template<unsigned DIM>
void oomph::Preconditioner::setup
inlinevirtual

Use the version in the Preconditioner base class for the alternative setup function that takes a matrix pointer as an argument.

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

◆ setup_coarsest_level_structures()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_coarsest_level_structures

Function to create the fully expanded system matrix on the coarsest level

Function to set up structures on the coarsest level of the MG hierarchy. This includes setting up the CRDoubleMatrix version of the coarsest level system matrix. This would otherwise be stored as a vector of pointers to the constituent CRDoubleMatrix objects which has the form: |--—| | A_r | Matrix_mg_pt = |--—| | A_i | |--—| and we want to construct: |--—|---—| | A_r | -A_c | Coarse_matrix_mg_pt = |--—|---—| | A_c | A_r | |--—|---—| Once this is done we have to set up the distributions of the vectors associated with Coarse_matrix_mg_pt

1981  {
1982  // Start timer
1983  double t_cm_start = TimingHelpers::timer();
1984 
1985  //---------------------------------------------------
1986  // Extract information from the constituent matrices:
1987  //---------------------------------------------------
1988 
1989  // Grab the real and imaginary matrix parts from storage
1990  CRDoubleMatrix* real_matrix_pt = Mg_matrices_storage_pt[Nlevel - 1][0];
1991  CRDoubleMatrix* imag_matrix_pt = Mg_matrices_storage_pt[Nlevel - 1][1];
1992 
1993  // Number of nonzero entries in A_r
1994  unsigned nnz_r = real_matrix_pt->nnz();
1995  unsigned nnz_c = imag_matrix_pt->nnz();
1996 
1997  // Calculate the total number of rows (and thus columns) in the real matrix
1998  unsigned n_rows_r = real_matrix_pt->nrow();
1999 
2000  // Acquire access to the value, row_start and column_index arrays from
2001  // the real and imaginary portions of the full matrix.
2002  // Real part:
2003  const double* value_r_pt = real_matrix_pt->value();
2004  const int* row_start_r_pt = real_matrix_pt->row_start();
2005  const int* column_index_r_pt = real_matrix_pt->column_index();
2006 
2007  // Imaginary part:
2008  const double* value_c_pt = imag_matrix_pt->value();
2009  const int* row_start_c_pt = imag_matrix_pt->row_start();
2010  const int* column_index_c_pt = imag_matrix_pt->column_index();
2011 
2012 #ifdef PARANOID
2013  // PARANOID check - make sure the matrices have the same number of rows
2014  // to ensure they are compatible
2015  if (real_matrix_pt->nrow() != imag_matrix_pt->nrow())
2016  {
2017  std::ostringstream error_message_stream;
2018  error_message_stream << "The real and imaginary matrices do not have "
2019  << "compatible sizes";
2020  throw OomphLibError(error_message_stream.str(),
2021  OOMPH_CURRENT_FUNCTION,
2022  OOMPH_EXCEPTION_LOCATION);
2023  }
2024  // PARANOID check - make sure the matrices have the same number of columns
2025  // to ensure they are compatible
2026  if (real_matrix_pt->ncol() != imag_matrix_pt->ncol())
2027  {
2028  std::ostringstream error_message_stream;
2029  error_message_stream << "The real and imaginary matrices do not have "
2030  << "compatible sizes";
2031  throw OomphLibError(error_message_stream.str(),
2032  OOMPH_CURRENT_FUNCTION,
2033  OOMPH_EXCEPTION_LOCATION);
2034  }
2035 #endif
2036 
2037  // Calculate the total number of nonzeros in the full matrix
2038  unsigned nnz = 2 * (nnz_r + nnz_c);
2039 
2040  // Allocate space for the row_start and column_index vectors associated with
2041  // the complete matrix
2042  Vector<double> value(nnz, 0.0);
2043  Vector<int> column_index(nnz, 0);
2044  Vector<int> row_start(2 * n_rows_r + 1, 0);
2045 
2046  //----------------------------
2047  // Build the row start vector:
2048  //----------------------------
2049 
2050  // Loop over the rows of the row_start vector. This is decomposed into
2051  // two parts. The first (n_rows_r+1) entries are found by computing
2052  // the entry-wise addition of row_start_r+row_start_c. The remaining
2053  // entries are almost the same as the first (n_rows_r+1). The only
2054  // distinction is that we need to shift the values of the entries by
2055  // the number of nonzeros in the top half of A. This is obvious from
2056  // observing the structure of the complete matrix.
2057 
2058  // Loop over the rows in the top half:
2059  for (unsigned i = 0; i < n_rows_r; i++)
2060  {
2061  // Set the i-th entry in the row start vector
2062  row_start[i] = *(row_start_r_pt + i) + *(row_start_c_pt + i);
2063  }
2064 
2065  // Loop over the rows in the bottom half:
2066  for (unsigned i = n_rows_r; i < 2 * n_rows_r; i++)
2067  {
2068  // Set the i-th entry in the row start vector (bottom half)
2069  row_start[i] = row_start[i - n_rows_r] + (nnz_r + nnz_c);
2070  }
2071 
2072  // Set the last entry in the row start vector
2073  row_start[2 * n_rows_r] = nnz;
2074 
2075  //-----------------------------------------
2076  // Build the column index and value vector:
2077  //-----------------------------------------
2078 
2079  // Variable to store the index of the nonzero
2080  unsigned i_nnz = 0;
2081 
2082  // Loop over the top half of the complete matrix
2083  for (unsigned i = 0; i < n_rows_r; i++)
2084  {
2085  // Calculate the number of nonzeros in the i-th row of A_r
2086  unsigned i_row_r_nnz = *(row_start_r_pt + i + 1) - *(row_start_r_pt + i);
2087 
2088  // Calculate the number of nonzeros in the i-th row of A_c
2089  unsigned i_row_c_nnz = *(row_start_c_pt + i + 1) - *(row_start_c_pt + i);
2090 
2091  // The index of the first entry in the i-th row of A_r
2092  unsigned i_first_entry_r = *(row_start_r_pt + i);
2093 
2094  // The index of the first entry in the i-th row of A_c
2095  unsigned i_first_entry_c = *(row_start_c_pt + i);
2096 
2097  // Loop over the number of nonzeros in the row associated with A_r
2098  for (unsigned j = 0; j < i_row_r_nnz; j++)
2099  {
2100  // Assign the column index of the j-th entry in the i-th row of A_r
2101  // to the column_index vector
2102  column_index[i_nnz] = *(column_index_r_pt + i_first_entry_r + j);
2103 
2104  // Assign the corresponding entry in the value vector
2105  value[i_nnz] = *(value_r_pt + i_first_entry_r + j);
2106 
2107  // Increment the value of i_nnz
2108  i_nnz++;
2109  }
2110 
2111  // Loop over the number of nonzeros in the row associated with -A_c
2112  for (unsigned j = 0; j < i_row_c_nnz; j++)
2113  {
2114  // Assign the column index of the j-th entry in the i-th row of -A_c
2115  // to the column_index vector
2116  column_index[i_nnz] =
2117  *(column_index_c_pt + i_first_entry_c + j) + n_rows_r;
2118 
2119  // Assign the corresponding entry in the value vector
2120  value[i_nnz] = -*(value_c_pt + i_first_entry_c + j);
2121 
2122  // Increment the value of i_nnz
2123  i_nnz++;
2124  }
2125  } // for (unsigned i=0;i<n_rows_r;i++)
2126 
2127  // Loop over the bottom half of the complete matrix
2128  for (unsigned i = n_rows_r; i < 2 * n_rows_r; i++)
2129  {
2130  // Calculate the number of nonzeros in row i of A_r
2131  unsigned i_row_r_nnz =
2132  *(row_start_r_pt + i - n_rows_r + 1) - *(row_start_r_pt + i - n_rows_r);
2133 
2134  // Calculate the number of nonzeros in row i of A_c
2135  unsigned i_row_c_nnz =
2136  *(row_start_c_pt + i - n_rows_r + 1) - *(row_start_c_pt + i - n_rows_r);
2137 
2138  // Get the index of the first entry in the i-th row of A_r
2139  unsigned i_first_entry_r = *(row_start_r_pt + i - n_rows_r);
2140 
2141  // Get the index of the first entry in the i-th row of A_c
2142  unsigned i_first_entry_c = *(row_start_c_pt + i - n_rows_r);
2143 
2144  // Loop over the number of nonzeros in the row associated with A_c
2145  for (unsigned j = 0; j < i_row_c_nnz; j++)
2146  {
2147  // Assign the column index of the j-th entry in the i-th row of A_c
2148  // to the column_index vector
2149  column_index[i_nnz] = *(column_index_c_pt + i_first_entry_c + j);
2150 
2151  // Assign the corresponding entry in the value vector
2152  value[i_nnz] = *(value_c_pt + i_first_entry_c + j);
2153 
2154  // Increment the value of i_nnz
2155  i_nnz++;
2156  }
2157 
2158  // Loop over the number of nonzeros in the row associated with A_r
2159  for (unsigned j = 0; j < i_row_r_nnz; j++)
2160  {
2161  // Assign the column index of the j-th entry in the i-th row of A_r
2162  // to the column_index vector
2163  column_index[i_nnz] =
2164  *(column_index_r_pt + i_first_entry_r + j) + n_rows_r;
2165 
2166  // Assign the corresponding entry in the value vector
2167  value[i_nnz] = *(value_r_pt + i_first_entry_r + j);
2168 
2169  // Increment the value of i_nnz
2170  i_nnz++;
2171  }
2172  } // for (unsigned i=n_rows_r;i<2*n_rows_r;i++)
2173 
2174  //-----------------------
2175  // Build the full matrix:
2176  //-----------------------
2177 
2178  // Allocate space for a CRDoubleMatrix
2179  Coarsest_matrix_mg_pt = new CRDoubleMatrix;
2180 
2181  // Make the distribution pointer
2182  LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
2183  Mg_hierarchy_pt[Nlevel - 1]->communicator_pt(), 2 * n_rows_r, false);
2184 
2185  // First, we need to build the matrix. Making use of its particular
2186  // structure we know that there are 2*n_rows_r columns in this matrix.
2187  // The remaining information has just been sorted out
2188  Coarsest_matrix_mg_pt->build(
2189  dist_pt, 2 * n_rows_r, value, column_index, row_start);
2190 
2191  // Build the distribution of associated solution vector
2192  Coarsest_x_mg.build(dist_pt);
2193 
2194  // Build the distribution of associated RHS vector
2195  Coarsest_rhs_mg.build(dist_pt);
2196 
2197  // Delete the associated distribution pointer
2198  delete dist_pt;
2199 
2200  // Summarise setup
2201  double t_cm_end = TimingHelpers::timer();
2202  double total_setup_time = double(t_cm_end - t_cm_start);
2203  oomph_info << " - Time for formation of the full matrix "
2204  << "on the coarsest level [sec]: " << total_setup_time << "\n"
2205  << std::endl;
2206  } // End of setup_coarsest_matrix_mg
oomph::CRDoubleMatrix::build
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
Definition: matrices.cc:1672
oomph::DoubleVector::build
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
Definition: double_vector.cc:35

References oomph::CRDoubleMatrix::column_index(), i, j, oomph::CRDoubleMatrix::ncol(), oomph::CRDoubleMatrix::nnz(), oomph::CRDoubleMatrix::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::CRDoubleMatrix::row_start(), oomph::TimingHelpers::timer(), Eigen::value, and oomph::CRDoubleMatrix::value().

◆ setup_interpolation_matrices()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_interpolation_matrices

Setup the interpolation matrix on each level.

Set up the interpolation matrices.

2745  {
2746  // Variable to hold the number of sons in an element
2747  unsigned n_sons;
2748 
2749  // Number of son elements
2750  if (DIM == 2)
2751  {
2752  n_sons = 4;
2753  }
2754  else if (DIM == 3)
2755  {
2756  n_sons = 8;
2757  }
2758  else
2759  {
2760  std::ostringstream error_message_stream;
2761  error_message_stream << "DIM should be 2 or 3 not " << DIM << std::endl;
2762  throw OomphLibError(error_message_stream.str(),
2763  OOMPH_CURRENT_FUNCTION,
2764  OOMPH_EXCEPTION_LOCATION);
2765  }
2766 
2767 #ifdef PARANOID
2768  // PARANOID check - for our implementation we demand that the first
2769  // submesh is the refineable bulk mesh that is provided by the function
2770  // mg_bulk_mesh_pt. Since each mesh in the hierarchy will be set up
2771  // in the same manner through the problem constructor we only needed
2772  // to check the finest level mesh
2773  if (Mg_hierarchy_pt[0]->mesh_pt(0) != Mg_hierarchy_pt[0]->mg_bulk_mesh_pt())
2774  {
2775  // Create an output stream
2776  std::ostringstream error_message_stream;
2777 
2778  // Create the error message
2779  error_message_stream << "MG solver requires the first submesh be the "
2780  << "refineable bulk mesh provided by the user."
2781  << std::endl;
2782 
2783  // Throw the error message
2784  throw OomphLibError(error_message_stream.str(),
2785  OOMPH_CURRENT_FUNCTION,
2786  OOMPH_EXCEPTION_LOCATION);
2787  }
2788 #endif
2789 
2790  // Vector of local coordinates in the element
2791  Vector<double> s(DIM, 0.0);
2792 
2793  // Vector to contain the (Eulerian) spatial location of the fine node.
2794  // This will only be used if we have a PML layer attached to the mesh
2795  // which will require the use of locate_zeta functionality
2796  Vector<double> fine_node_position(DIM, 0.0);
2797 
2798  // Find the number of nodes in an element in the mesh. Since this
2799  // quantity will be the same in all meshes we can precompute it here
2800  // using the original problem pointer
2801  unsigned nnod_element =
2802  dynamic_cast<FiniteElement*>(Mg_problem_pt->mesh_pt()->element_pt(0))
2803  ->nnode();
2804 
2805  // Initialise a null pointer which will be used to point to an object
2806  // of the class MeshAsGeomObject
2807  MeshAsGeomObject* coarse_mesh_from_obj_pt = 0;
2808 
2809  // Loop over each level (apart from the coarsest level; an interpolation
2810  // matrix and thus a restriction matrix is not needed here), with 0 being
2811  // the finest level and Nlevel-1 being the coarsest level
2812  for (unsigned level = 0; level < Nlevel - 1; level++)
2813  {
2814  // Store the ID of the fine level
2815  unsigned fine_level = level;
2816 
2817  // Store the ID of the coarse level
2818  unsigned coarse_level = level + 1;
2819 
2820  // Make a pointer to the mesh on the finer level and dynamic_cast
2821  // it as an object of the refineable mesh class.
2822  Mesh* ref_fine_mesh_pt = Mg_hierarchy_pt[fine_level]->mesh_pt();
2823 
2824  // Make a pointer to the mesh on the coarse level and dynamic_cast
2825  // it as an object of the refineable mesh class
2826  Mesh* ref_coarse_mesh_pt = Mg_hierarchy_pt[coarse_level]->mesh_pt();
2827 
2828  // Access information about the number of elements in the fine mesh
2829  // from the pointer to the fine mesh (to loop over the rows of the
2830  // interpolation matrix)
2831  unsigned fine_n_element = ref_fine_mesh_pt->nelement();
2832 
2833  // Find the number of elements in the bulk mesh
2834  unsigned n_bulk_mesh_element =
2835  Mg_hierarchy_pt[fine_level]->mg_bulk_mesh_pt()->nelement();
2836 
2837  // If there are elements apart from those in the bulk mesh
2838  // then we require locate_zeta functionality. For this reason
2839  // we have to assign a value to the MeshAsGeomObject pointer
2840  // which we do by creating a MeshAsGeomObject from the coarse
2841  // mesh pointer
2842  if (fine_n_element > n_bulk_mesh_element)
2843  {
2844  // Assign the pointer to coarse_mesh_from_obj_pt
2845  coarse_mesh_from_obj_pt =
2846  new MeshAsGeomObject(Mg_hierarchy_pt[coarse_level]->mesh_pt());
2847  }
2848 
2849  // Find the number of dofs in the fine mesh
2850  unsigned n_fine_dof = Mg_hierarchy_pt[fine_level]->ndof();
2851 
2852  // Find the number of dofs in the coarse mesh
2853  unsigned n_coarse_dof = Mg_hierarchy_pt[coarse_level]->ndof();
2854 
2855  // To get the number of rows in the interpolation matrix we divide
2856  // the number of dofs on each level by 2 since there are 2 dofs at
2857  // each node in the mesh corresponding to the real and imaginary
2858  // part of the solution:
2859 
2860  // Get the numbers of rows in the interpolation matrix.
2861  unsigned n_row = n_fine_dof / 2.0;
2862 
2863  // Get the numbers of columns in the interpolation matrix
2864  unsigned n_col = n_coarse_dof / 2.0;
2865 
2866  // Mapping relating the pointers to related elements in the coarse
2867  // and fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]. If
2868  // the element in the fine mesh has been unrefined between these two
2869  // levels, this map returns the father element in the coarsened mesh.
2870  // If this element in the fine mesh has not been unrefined, the map
2871  // returns the pointer to the same-sized equivalent element in the
2872  // coarsened mesh.
2873  std::map<RefineableQElement<DIM>*, RefineableQElement<DIM>*>
2874  coarse_mesh_reference_element_pt;
2875 
2876  // Variable to count the number of elements in the fine mesh
2877  unsigned e_fine = 0;
2878 
2879  // Variable to count the number of elements in the coarse mesh
2880  unsigned e_coarse = 0;
2881 
2882  // While loop over fine elements (while because we're incrementing the
2883  // counter internally if the element was unrefined...). We only bother
2884  // going until we've covered all of the elements in the bulk mesh
2885  // because once we reach the PML layer (if there is one) there will be
2886  // no tree structure to match in between levels as PML elements are
2887  // simply regenerated once the bulk mesh has been refined.
2888  while (e_fine < n_bulk_mesh_element)
2889  {
2890  // Pointer to element in fine mesh
2891  RefineableQElement<DIM>* el_fine_pt =
2892  dynamic_cast<RefineableQElement<DIM>*>(
2893  ref_fine_mesh_pt->finite_element_pt(e_fine));
2894 
2895 #ifdef PARANOID
2896  // Make sure the coarse level element pointer is not a null pointer. If
2897  // it is something has gone wrong
2898  if (e_coarse > ref_coarse_mesh_pt->nelement() - 1)
2899  {
2900  // Create an output stream
2901  std::ostringstream error_message_stream;
2902 
2903  // Create an error message
2904  error_message_stream
2905  << "The coarse level mesh has " << ref_coarse_mesh_pt->nelement()
2906  << " elements but the coarse\nelement loop "
2907  << "is looking at the " << e_coarse << "-th element!" << std::endl;
2908 
2909  // Throw the error message
2910  throw OomphLibError(error_message_stream.str(),
2911  OOMPH_CURRENT_FUNCTION,
2912  OOMPH_EXCEPTION_LOCATION);
2913  }
2914 #endif
2915 
2916  // Pointer to element in coarse mesh
2917  RefineableQElement<DIM>* el_coarse_pt =
2918  dynamic_cast<RefineableQElement<DIM>*>(
2919  ref_coarse_mesh_pt->finite_element_pt(e_coarse));
2920 
2921  // If the levels are different then the element in the fine
2922  // mesh has been unrefined between these two levels
2923  if (el_fine_pt->tree_pt()->level() > el_coarse_pt->tree_pt()->level())
2924  {
2925  // The element in the fine mesh has been unrefined between these two
2926  // levels. Hence it and its three brothers (ASSUMED to be stored
2927  // consecutively in the Mesh's vector of pointers to its constituent
2928  // elements -- we'll check this!) share the same father element in
2929  // the coarse mesh, currently pointed to by el_coarse_pt.
2930  for (unsigned i = 0; i < n_sons; i++)
2931  {
2932  // Set mapping to father element in coarse mesh
2933  coarse_mesh_reference_element_pt
2934  [dynamic_cast<RefineableQElement<DIM>*>(
2935  ref_fine_mesh_pt->finite_element_pt(e_fine))] = el_coarse_pt;
2936 
2937  // Increment counter for elements in fine mesh
2938  e_fine++;
2939  }
2940  }
2941  // The element in the fine mesh has not been unrefined between
2942  // these two levels, so the reference element is the same-sized
2943  // equivalent element in the coarse mesh
2944  else if (el_fine_pt->tree_pt()->level() ==
2945  el_coarse_pt->tree_pt()->level())
2946  {
2947  // The element in the fine mesh has not been unrefined between these
2948  // two levels
2949  coarse_mesh_reference_element_pt
2950  [dynamic_cast<RefineableQElement<DIM>*>(
2951  ref_fine_mesh_pt->finite_element_pt(e_fine))] = el_coarse_pt;
2952 
2953  // Increment counter for elements in fine mesh
2954  e_fine++;
2955  }
2956  // If the element has been unrefined between levels. Not something
2957  // we can deal with at the moment (although it would be relatively
2958  // simply to implement...).
2959  // Option 1: Find the son elements (from the coarse mesh) associated
2960  // with the father element (from the fine mesh) and extract the
2961  // appropriate nodal values to find the nodal values in the father
2962  // element.
2963  // Option 2: Use the function
2964  // unrefine_uniformly();
2965  // to ensure that the coarser meshes really only have father elements
2966  // or the element itself.
2967  else
2968  {
2969  // Create an output stream
2970  std::ostringstream error_message_stream;
2971 
2972  // Create an error message
2973  error_message_stream << "Element unrefined between levels! Can't "
2974  << "handle this case yet..." << std::endl;
2975 
2976  // Throw the error message
2977  throw OomphLibError(error_message_stream.str(),
2978  OOMPH_CURRENT_FUNCTION,
2979  OOMPH_EXCEPTION_LOCATION);
2980  } // if (el_fine_pt->tree_pt()->level()!=...)
2981 
2982  // Increment counter for elements in coarse mesh
2983  e_coarse++;
2984  } // while(e_fine<n_bulk_mesh_element)
2985 
2986  // To allow an update of a row only once we use STL vectors for bools
2987  std::vector<bool> contribution_made(n_row, false);
2988 
2989  // Create the row start vector whose i-th entry will contain the index
2990  // in column_start where the entries in the i-th row of the matrix start
2991  Vector<int> row_start(n_row + 1);
2992 
2993  // Create the column index vector whose entries will store the column
2994  // index of each contribution, i.e. the global equation of the coarse
2995  // mesh node which supplies a contribution. We don't know how many
2996  // entries this will have so we dynamically allocate entries at run time
2997  Vector<int> column_index;
2998 
2999  // Create the value vector which will be used to store the nonzero
3000  // entries in the CRDoubleMatrix. We don't know how many entries this
3001  // will have either so we dynamically allocate its entries at run time
3002  Vector<double> value;
3003 
3004  // The value of index will tell us which row of the interpolation matrix
3005  // we're working on in the following for loop
3006  // DRAIG: Should this worry us? We assume that every node we cover is
3007  // the next node in the mesh (we loop over the elements and the nodes
3008  // inside that). This does work but it may not work for some meshes
3009  unsigned index = 0;
3010 
3011  // New loop to go over each element in the fine mesh
3012  for (unsigned k = 0; k < fine_n_element; k++)
3013  {
3014  // Set a pointer to the element in the fine mesh
3015  RefineableQElement<DIM>* el_fine_pt =
3016  dynamic_cast<RefineableQElement<DIM>*>(
3017  ref_fine_mesh_pt->finite_element_pt(k));
3018 
3019  // Initialise a null pointer which will point to the corresponding
3020  // coarse mesh element. If we're in the bulk mesh, this will point
3021  // to the parent element of the fine mesh element or itself (if
3022  // there was no refinement between the levels). If, however, we're
3023  // in the PML layer then this will point to element that encloses
3024  // the fine mesh PML element
3025  RefineableQElement<DIM>* el_coarse_pt = 0;
3026 
3027  // Variable to store the son type of the fine mesh element with
3028  // respect to the corresponding coarse mesh element (set to OMEGA
3029  // if no unrefinement took place)
3030  int son_type = 0;
3031 
3032  // If we're in the bulk mesh we can assign the coarse mesh element
3033  // pointer right now using the map coarse_mesh_reference_element_pt
3034  if (k < n_bulk_mesh_element)
3035  {
3036  // Get the reference element (either the father element or the
3037  // same-sized element) in the coarse mesh
3038  el_coarse_pt = coarse_mesh_reference_element_pt[el_fine_pt];
3039 
3040  // Check if the elements are different on both levels (i.e. to check
3041  // if any unrefinement took place)
3042  if (el_fine_pt->tree_pt()->level() !=
3043  el_coarse_pt->tree_pt()->level())
3044  {
3045  // If the element was refined then we need the son type
3046  son_type = el_fine_pt->tree_pt()->son_type();
3047  }
3048  else
3049  {
3050  // If the element is the same in both meshes
3051  son_type = Tree::OMEGA;
3052  }
3053  } // if (k<n_bulk_mesh_element)
3054 
3055  // Loop through all the nodes in an element in the fine mesh
3056  for (unsigned i = 0; i < nnod_element; i++)
3057  {
3058  // Row number in interpolation matrix: Global equation number
3059  // of the d.o.f. stored at this node in the fine element
3060  int ii = el_fine_pt->node_pt(i)->eqn_number(0);
3061 
3062  // Check whether or not the node is a proper d.o.f.
3063  if (ii >= 0)
3064  {
3065  // Only assign values to the given row of the interpolation matrix
3066  // if they haven't already been assigned. Notice, we check if the
3067  // (ii/2)-th entry has been set. This is because there are two dofs
3068  // at each node so dividing by two will give us the row number
3069  // associated with this node
3070  if (contribution_made[ii / 2] == false)
3071  {
3072  // The value of index was initialised when we allocated space
3073  // for the three vectors to store the CRDoubleMatrix. We use
3074  // index to go through the entries of the row_start vector.
3075  // The next row starts at the value.size() (draw out the entries
3076  // in value if this doesn't make sense noting that the storage
3077  // for the vector 'value' is dynamically allocated)
3078  row_start[index] = value.size();
3079 
3080  // If we're in the PML layer then we need to find the parent
3081  // element associated with a given node using locate_zeta.
3082  // Unfortunately, if this is the case and a given local node in
3083  // the fine mesh element lies on the interface between two
3084  // elements (E1) and (E2) in the coarse mesh then either element
3085  // will be returned. Suppose, for instance, a pointer to (E1) is
3086  // returned. If the next local node lies in the (E2) then the
3087  // contributions which we pick up from the local nodes in (E1)
3088  // will most likely be incorrect while the column index (the
3089  // global equation number of the coarse mesh node) corresponding
3090  // to the given contribution will definitely be incorrect. If
3091  // we're not using linear quad elements we can get around this by
3092  // using locate_zeta carefully by identifying a node which is
3093  // internal to the fine mesh element. This will allow us to
3094  // correctly obtain the corresponding element in the coarse mesh
3095  if (k >= n_bulk_mesh_element)
3096  {
3097  // Get the (Eulerian) spatial location of the i-th local node
3098  // in the fine mesh element
3099  el_fine_pt->node_pt(i)->position(fine_node_position);
3100 
3101  // Initalise a null pointer to point to an object of the class
3102  // GeomObject
3103  GeomObject* el_pt = 0;
3104 
3105  // Get the reference element (either the father element or the
3106  // same-sized element) in the coarse-mesh PML layer using
3107  // the locate_zeta functionality
3108  coarse_mesh_from_obj_pt->locate_zeta(
3109  fine_node_position, el_pt, s);
3110 
3111  // Upcast the GeomElement object to a RefineableQElement object
3112  el_coarse_pt = dynamic_cast<RefineableQElement<DIM>*>(el_pt);
3113  }
3114  else
3115  {
3116  // Calculate the local coordinates of the given node
3117  el_fine_pt->local_coordinate_of_node(i, s);
3118 
3119  // Find the local coordinates s, of the present node, in the
3120  // reference element in the coarse mesh, given the element's
3121  // son_type (e.g. SW,SE... )
3122  level_up_local_coord_of_node(son_type, s);
3123  }
3124 
3125  // Allocate space for shape functions in the coarse mesh
3126  Shape psi(el_coarse_pt->nnode());
3127 
3128  // Set the shape function in the reference element
3129  el_coarse_pt->shape(s, psi);
3130 
3131  // Auxiliary storage
3132  std::map<unsigned, double> contribution;
3133  Vector<unsigned> keys;
3134 
3135  // Find the number of nodes in an element in the coarse mesh
3136  unsigned nnod_coarse = el_coarse_pt->nnode();
3137 
3138  // Loop through all the nodes in the reference element
3139  for (unsigned j = 0; j < nnod_coarse; j++)
3140  {
3141  // Column number in interpolation matrix: Global equation
3142  // number of the d.o.f. stored at this node in the coarse
3143  // element
3144  int jj = el_coarse_pt->node_pt(j)->eqn_number(0);
3145 
3146  // If the value stored at this node is pinned or hanging
3147  if (jj < 0)
3148  {
3149  // Hanging node: In this case we need to accumulate the
3150  // contributions from the master nodes
3151  if (el_coarse_pt->node_pt(j)->is_hanging())
3152  {
3153  // Find the number of master nodes of the hanging
3154  // the node in the reference element
3155  HangInfo* hang_info_pt =
3156  el_coarse_pt->node_pt(j)->hanging_pt();
3157  unsigned nmaster = hang_info_pt->nmaster();
3158 
3159  // Loop over the master nodes
3160  for (unsigned i_master = 0; i_master < nmaster; i_master++)
3161  {
3162  // Set up a pointer to the master node
3163  Node* master_node_pt =
3164  hang_info_pt->master_node_pt(i_master);
3165 
3166  // The column number in the interpolation matrix: the
3167  // global equation number of the d.o.f. stored at this
3168  // master node for the coarse element
3169  int master_jj = master_node_pt->eqn_number(0);
3170 
3171  // Is the master node a proper d.o.f.?
3172  if (master_jj >= 0)
3173  {
3174  // If the weight of the master node is non-zero
3175  if (psi(j) * hang_info_pt->master_weight(i_master) !=
3176  0.0)
3177  {
3178  contribution[master_jj / 2] +=
3179  psi(j) * hang_info_pt->master_weight(i_master);
3180  }
3181  } // if (master_jj>=0)
3182  } // for (unsigned i_master=0;i_master<nmaster;i_master++)
3183  } // if (el_coarse_pt->node_pt(j)->is_hanging())
3184  }
3185  // In the case that the node is not pinned or hanging
3186  else
3187  {
3188  // If we can get a nonzero contribution from the shape
3189  // function
3190  if (psi(j) != 0.0)
3191  {
3192  contribution[jj / 2] += psi(j);
3193  }
3194  } // if (jj<0) else
3195  } // for (unsigned j=0;j<nnod_coarse;j++)
3196 
3197  // Put the contributions into the value vector
3198  for (std::map<unsigned, double>::iterator it =
3199  contribution.begin();
3200  it != contribution.end();
3201  ++it)
3202  {
3203  if (it->second != 0)
3204  {
3205  // The first entry in contribution is the column index
3206  column_index.push_back(it->first);
3207 
3208  // The second entry in contribution is the value of the
3209  // contribution
3210  value.push_back(it->second);
3211  }
3212  } // for (std::map<unsigned,double>::iterator it=...)
3213 
3214  // Increment the value of index by 1 to indicate that we have
3215  // finished the row, or equivalently, covered another global
3216  // node in the fine mesh
3217  index++;
3218 
3219  // Change the entry in contribution_made to true now to indicate
3220  // that the row has been filled
3221  contribution_made[ii / 2] = true;
3222  } // if(contribution_made[ii/2]==false)
3223  } // if (ii>=0)
3224  } // for(unsigned i=0;i<nnod_element;i++)
3225  } // for (unsigned k=0;k<fine_n_element;k++)
3226 
3227  // Set the last entry in the row_start vector
3228  row_start[n_row] = value.size();
3229 
3230  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
3231  // using the vectors value, row_start, column_index and the value
3232  // of fine_n_unknowns and coarse_n_unknowns
3233  interpolation_matrix_set(
3234  level, value, column_index, row_start, n_col, n_row);
3235  } // for (unsigned level=0;level<Nlevel-1;level++)
3236  } // End of setup_interpolation_matrices
oomph::HelmholtzMGPreconditioner::interpolation_matrix_set
void interpolation_matrix_set(const unsigned &level, double *value, int *col_index, int *row_st, unsigned &ncol, unsigned &nnz)
Definition: helmholtz_geometric_multigrid.h:396
oomph::HelmholtzMGPreconditioner::level_up_local_coord_of_node
void level_up_local_coord_of_node(const int &son_type, Vector< double > &s)
k
char char char int int * k
Definition: level2_impl.h:374

References DIM, oomph::Data::eqn_number(), oomph::Mesh::finite_element_pt(), i, j, k, oomph::MeshAsGeomObject::locate_zeta(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::Mesh::nelement(), oomph::HangInfo::nmaster(), oomph::Tree::OMEGA, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, s, and Eigen::value.

◆ setup_interpolation_matrices_unstructured()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_interpolation_matrices_unstructured

Setup the interpolation matrices.

Setup the interpolation matrix on each level (used for unstructured meshes)

3244  {
3245 #ifdef PARANOID
3246  if ((DIM != 2) && (DIM != 3))
3247  {
3248  std::ostringstream error_message_stream;
3249  error_message_stream << "DIM should be 2 or 3 not " << DIM << std::endl;
3250  throw OomphLibError(error_message_stream.str(),
3251  OOMPH_CURRENT_FUNCTION,
3252  OOMPH_EXCEPTION_LOCATION);
3253  }
3254 #endif
3255 
3256  // Vector of local coordinates in the element
3257  Vector<double> s(DIM, 0.0);
3258 
3259  // Loop over each level (apart from the coarsest level; an interpolation
3260  // matrix and thus a restriction matrix is not needed here), with 0 being
3261  // the finest level and Nlevel-1 being the coarsest level
3262  for (unsigned level = 0; level < Nlevel - 1; level++)
3263  {
3264  // Assign values to a couple of variables to help out
3265  unsigned coarse_level = level + 1;
3266  unsigned fine_level = level;
3267 
3268  // Make a pointer to the mesh on the finer level and dynamic_cast
3269  // it as an object of the refineable mesh class
3270  Mesh* ref_fine_mesh_pt = Mg_hierarchy_pt[fine_level]->mg_bulk_mesh_pt();
3271 
3272  // To use the locate zeta functionality the coarse mesh must be of the
3273  // type MeshAsGeomObject
3274  MeshAsGeomObject* coarse_mesh_from_obj_pt =
3275  new MeshAsGeomObject(Mg_hierarchy_pt[coarse_level]->mg_bulk_mesh_pt());
3276 
3277  // Access information about the number of degrees of freedom
3278  // from the pointers to the problem on each level
3279  unsigned coarse_n_unknowns = Mg_hierarchy_pt[coarse_level]->ndof();
3280  unsigned fine_n_unknowns = Mg_hierarchy_pt[fine_level]->ndof();
3281 
3282  // Make storage vectors to form the interpolation matrix using a
3283  // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
3284  // row_start contains entries which tells us at which entry of the
3285  // vector column_start we start on the i-th row of the actual matrix.
3286  // The entries of column_start indicate the column position of each
3287  // non-zero entry in every row. This runs through the first row
3288  // from left to right then the second row (again, left to right)
3289  // and so on until the end. The entries in value are the entries in
3290  // the interpolation matrix. The vector column_start has the same length
3291  // as the vector value.
3292  Vector<double> value;
3293  Vector<int> column_index;
3294  Vector<int> row_start(fine_n_unknowns + 1);
3295 
3296  // Vector to contain the (Eulerian) spatial location of the fine node
3297  Vector<double> fine_node_position(DIM);
3298 
3299  // Find the number of nodes in the mesh
3300  unsigned n_node_fine_mesh = ref_fine_mesh_pt->nnode();
3301 
3302  // Loop over the unknowns in the mesh
3303  for (unsigned i_fine_node = 0; i_fine_node < n_node_fine_mesh;
3304  i_fine_node++)
3305  {
3306  // Set a pointer to the i_fine_unknown-th node in the fine mesh
3307  Node* fine_node_pt = ref_fine_mesh_pt->node_pt(i_fine_node);
3308 
3309  // Get the global equation number
3310  int i_fine = fine_node_pt->eqn_number(0);
3311 
3312  // If the node is a proper d.o.f.
3313  if (i_fine >= 0)
3314  {
3315  // Row number in interpolation matrix: Global equation number
3316  // of the d.o.f. stored at this node in the fine element
3317  row_start[i_fine] = value.size();
3318 
3319  // Get the (Eulerian) spatial location of the fine node
3320  fine_node_pt->position(fine_node_position);
3321 
3322  // Create a null pointer to the GeomObject class
3323  GeomObject* el_pt = 0;
3324 
3325  // Get the reference element (either the father element or the
3326  // same-sized element) in the coarse mesh using locate_zeta
3327  coarse_mesh_from_obj_pt->locate_zeta(fine_node_position, el_pt, s);
3328 
3329  // Upcast GeomElement as a FiniteElement
3330  FiniteElement* el_coarse_pt = dynamic_cast<FiniteElement*>(el_pt);
3331 
3332  // Find the number of nodes in the element
3333  unsigned n_node = el_coarse_pt->nnode();
3334 
3335  // Allocate space for shape functions in the coarse mesh
3336  Shape psi(n_node);
3337 
3338  // Calculate the geometric shape functions at local coordinate s
3339  el_coarse_pt->shape(s, psi);
3340 
3341  // Auxiliary storage
3342  std::map<unsigned, double> contribution;
3343  Vector<unsigned> keys;
3344 
3345  // Loop through all the nodes in the (coarse mesh) element containing
3346  // the node pointed to by fine_node_pt (fine mesh)
3347  for (unsigned j_node = 0; j_node < n_node; j_node++)
3348  {
3349  // Get the j_coarse_unknown-th node in the coarse element
3350  Node* coarse_node_pt = el_coarse_pt->node_pt(j_node);
3351 
3352  // Column number in interpolation matrix: Global equation number of
3353  // the d.o.f. stored at this node in the coarse element
3354  int j_coarse = coarse_node_pt->eqn_number(0);
3355 
3356  // If the value stored at this node is pinned or hanging
3357  if (j_coarse < 0)
3358  {
3359  // Hanging node: In this case we need to accumulate the
3360  // contributions from the master nodes
3361  if (el_coarse_pt->node_pt(j_node)->is_hanging())
3362  {
3363  // Find the number of master nodes of the hanging
3364  // the node in the reference element
3365  HangInfo* hang_info_pt = coarse_node_pt->hanging_pt();
3366  unsigned nmaster = hang_info_pt->nmaster();
3367 
3368  // Loop over the master nodes
3369  for (unsigned i_master = 0; i_master < nmaster; i_master++)
3370  {
3371  // Set up a pointer to the master node
3372  Node* master_node_pt = hang_info_pt->master_node_pt(i_master);
3373 
3374  // The column number in the interpolation matrix: the
3375  // global equation number of the d.o.f. stored at this master
3376  // node for the coarse element
3377  int master_jj = master_node_pt->eqn_number(0);
3378 
3379  // Is the master node a proper d.o.f.?
3380  if (master_jj >= 0)
3381  {
3382  // If the weight of the master node is non-zero
3383  if (psi(j_node) * hang_info_pt->master_weight(i_master) !=
3384  0.0)
3385  {
3386  contribution[master_jj] +=
3387  psi(j_node) * hang_info_pt->master_weight(i_master);
3388  }
3389  } // End of if statement (check it's not a boundary node)
3390  } // End of the loop over the master nodes
3391  } // End of the if statement for only hanging nodes
3392  } // End of the if statement for pinned or hanging nodes
3393  // In the case that the node is not pinned or hanging
3394  else
3395  {
3396  // If we can get a nonzero contribution from the shape function
3397  // at the j_node-th node in the element
3398  if (psi(j_node) != 0.0)
3399  {
3400  contribution[j_coarse] += psi(j_node);
3401  }
3402  } // End of the if-else statement (check if the node was
3403  // pinned/hanging)
3404  } // Finished loop over the nodes j in the reference element (coarse)
3405 
3406  // Put the contributions into the value vector
3407  for (std::map<unsigned, double>::iterator it = contribution.begin();
3408  it != contribution.end();
3409  ++it)
3410  {
3411  if (it->second != 0)
3412  {
3413  value.push_back(it->second);
3414  column_index.push_back(it->first);
3415  }
3416  } // End of putting contributions into the value vector
3417  } // End check (whether or not the fine node was a d.o.f.)
3418  } // End of the for-loop over nodes in the fine mesh
3419 
3420  // Set the last entry of row_start
3421  row_start[fine_n_unknowns] = value.size();
3422 
3423  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
3424  // using the vectors value, row_start, column_index and the value
3425  // of fine_n_unknowns and coarse_n_unknowns
3426  interpolation_matrix_set(level,
3427  value,
3428  column_index,
3429  row_start,
3430  coarse_n_unknowns,
3431  fine_n_unknowns);
3432  } // End of loop over each level
3433  } // End of setup_interpolation_matrices_unstructured

References DIM, oomph::Data::eqn_number(), oomph::Node::hanging_pt(), oomph::Node::is_hanging(), oomph::MeshAsGeomObject::locate_zeta(), oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::Mesh::nnode(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::Node::position(), s, oomph::FiniteElement::shape(), and Eigen::value.

◆ setup_mg_hierarchy()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_mg_hierarchy
private

Function to set up the hierachy of levels. Creates a vector of pointers to each MG level

Set up the MG hierarchy. Creates a vector of pointers to each MG level and resizes internal storage for multigrid data

1264  {
1265  // Initialise the timer start variable
1266  double t_m_start = 0.0;
1267 
1268  // Notify the user if it is allowed
1269  if (!Suppress_all_output)
1270  {
1271  // Notify user of progress
1272  oomph_info << "\n==============="
1273  << "Creating Multigrid Hierarchy"
1274  << "===============\n"
1275  << std::endl;
1276 
1277  // Start clock
1278  t_m_start = TimingHelpers::timer();
1279  }
1280 
1281  // Create a bool to indicate whether or not we could create an unrefined
1282  // copy. This bool will be assigned the value FALSE when the current copy
1283  // is the last level of the multigrid hierarchy
1284  bool managed_to_create_unrefined_copy = true;
1285 
1286  // Now keep making copies and try to make an unrefined copy of
1287  // the mesh
1288  unsigned level = 0;
1289 
1290  // Set up all of the levels by making a completely unrefined copy
1291  // of the problem using the function make_new_problem
1292  while (managed_to_create_unrefined_copy)
1293  {
1294  // Make a new object of the same type as the derived problem
1295  HelmholtzMGProblem* new_problem_pt = Mg_problem_pt->make_new_problem();
1296 
1297  // Do anything that needs to be done before we can refine the mesh
1298  new_problem_pt->actions_before_adapt();
1299 
1300  // To create the next level in the hierarchy we need to create a mesh
1301  // which matches the refinement of the current problem and then unrefine
1302  // the mesh. This can alternatively be done using the function
1303  // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
1304  // reference mesh to do precisely the above with
1305  managed_to_create_unrefined_copy =
1306  new_problem_pt->mg_bulk_mesh_pt()
1307  ->refine_base_mesh_as_in_reference_mesh_minus_one(
1308  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt());
1309 
1310  // If we were able to unrefine the problem on the current level
1311  // then add the unrefined problem to a vector of the levels
1312  if (managed_to_create_unrefined_copy)
1313  {
1314  // Another level has been created so increment the level counter
1315  level++;
1316 
1317  // If the documentation of everything has not been suppressed
1318  // then tell the user we managed to create another level
1319  if (!Suppress_all_output)
1320  {
1321  // Notify user that unrefinement was successful
1322  oomph_info << "Success! Level " << level
1323  << " has been created:" << std::endl;
1324  }
1325 
1326  // Do anything that needs to be done after refinement
1327  new_problem_pt->actions_after_adapt();
1328 
1329  // Do the equation numbering for the new problem
1330  oomph_info << "\n - Number of equations: "
1331  << new_problem_pt->assign_eqn_numbers() << "\n"
1332  << std::endl;
1333 
1334  // Add the new problem pointer onto the vector of MG levels
1335  // and increment the value of level by 1
1336  Mg_hierarchy_pt.push_back(new_problem_pt);
1337  }
1338  // If we weren't able to create an unrefined copy
1339  else
1340  {
1341  // Delete the new problem
1342  delete new_problem_pt;
1343 
1344  // Make it a null pointer
1345  new_problem_pt = 0;
1346 
1347  // Assign the number of levels to Nlevel
1348  Nlevel = Mg_hierarchy_pt.size();
1349 
1350  // If we're allowed to document then tell the user we've reached
1351  // the coarsest level of the hierarchy
1352  if (!Suppress_all_output)
1353  {
1354  // Notify the user
1355  oomph_info << "Reached the coarsest level! "
1356  << "Number of levels: " << Nlevel << std::endl;
1357  }
1358  } // if (managed_to_create_unrefined_copy)
1359  } // while (managed_to_create_unrefined_copy)
1360 
1361  //------------------------------------------------------------------
1362  // Given that we know the number of levels in the hierarchy we can
1363  // resize the vectors which will store all the information required
1364  // for our solver:
1365  //------------------------------------------------------------------
1366  // Resize the vector storing all of the system matrices
1367  Mg_matrices_storage_pt.resize(Nlevel);
1368 
1369  // Resize the vector storing all of the solution vectors (X_mg)
1370  X_mg_vectors_storage.resize(Nlevel);
1371 
1372  // Resize the vector storing all of the RHS vectors (Rhs_mg)
1373  Rhs_mg_vectors_storage.resize(Nlevel);
1374 
1375  // Resize the vector storing all of the residual vectors
1376  Residual_mg_vectors_storage.resize(Nlevel);
1377 
1378  // Allocate space for the Max_edge_width vector, we only need the value
1379  // of h on the levels where we use a smoother
1380  Max_edge_width.resize(Nlevel - 1, 0.0);
1381 
1382  // Allocate space for the pre-smoother storage vector (remember, we do
1383  // not need a smoother on the coarsest level; we use a direct solve there)
1384  Pre_smoothers_storage_pt.resize(Nlevel - 1, 0);
1385 
1386  // Allocate space for the post-smoother storage vector (remember, we do
1387  // not need a smoother on the coarsest level; we use a direct solve there)
1388  Post_smoothers_storage_pt.resize(Nlevel - 1, 0);
1389 
1390  // Resize the vector storing all of the interpolation matrices
1391  Interpolation_matrices_storage_pt.resize(Nlevel - 1, 0);
1392 
1393  // Resize the vector storing all of the restriction matrices
1394  Restriction_matrices_storage_pt.resize(Nlevel - 1, 0);
1395 
1396  // If we're allowed to output information to the screen
1397  if (!Suppress_all_output)
1398  {
1399  // Stop clock
1400  double t_m_end = TimingHelpers::timer();
1401  double total_setup_time = double(t_m_end - t_m_start);
1402  oomph_info
1403  << "\nCPU time for creation of hierarchy of MG problems [sec]: "
1404  << total_setup_time << std::endl;
1405 
1406  // Notify user that the hierarchy of levels is complete
1407  oomph_info << "\n==============="
1408  << "Hierarchy Creation Complete"
1409  << "================\n"
1410  << std::endl;
1411  }
1412  } // End of setup_mg_hierarchy
oomph::HelmholtzMGPreconditioner::Max_edge_width
Vector< double > Max_edge_width
Definition: helmholtz_geometric_multigrid.h:677
oomph::HelmholtzMGProblem::make_new_problem
virtual HelmholtzMGProblem * make_new_problem()=0
oomph::Problem::actions_before_adapt
virtual void actions_before_adapt()
Definition: problem.h:1022

References oomph::Problem::actions_after_adapt(), oomph::Problem::actions_before_adapt(), oomph::Problem::assign_eqn_numbers(), oomph::HelmholtzMGProblem::make_new_problem(), oomph::HelmholtzMGProblem::mg_bulk_mesh_pt(), oomph::oomph_info, oomph::TreeBasedRefineableMeshBase::refine_base_mesh_as_in_reference_mesh_minus_one(), and oomph::TimingHelpers::timer().

◆ setup_mg_structures()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_mg_structures
private

Set up the MG structures on each level.

Function to set up the hierachy of levels. Creates a vector of pointers to each MG level

1477  {
1478  // Initialise the timer start variable
1479  double t_m_start = 0.0;
1480 
1481  // Start the clock (if we're allowed to time things)
1482  if (!Suppress_all_output)
1483  {
1484  // Start the clock
1485  t_m_start = TimingHelpers::timer();
1486  }
1487 
1488  // Calculate the wavenumber value:
1489  //--------------------------------
1490  // Reset the value of Wavenumber
1491  Wavenumber = 0.0;
1492 
1493  // Upcast the first element in the bulk mesh
1494  PMLHelmholtzEquations<DIM>* pml_helmholtz_el_pt =
1495  dynamic_cast<PMLHelmholtzEquations<DIM>*>(
1496  Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(0));
1497 
1498  // Grab the k_squared value from the element pointer and square root.
1499  // Note, we assume the wavenumber is the same everywhere in the mesh
1500  // and it is also the same on every level.
1501  Wavenumber = sqrt(pml_helmholtz_el_pt->k_squared());
1502 
1503  // We don't need the pointer anymore so make it a null pointer but don't
1504  // delete the underlying element data
1505  pml_helmholtz_el_pt = 0;
1506 
1507  // Set up the system matrix and accompanying vectors on each level:
1508  //-----------------------------------------------------------------
1509  // Loop over each level and extract the system matrix, solution vector
1510  // right-hand side vector and residual vector (to store the value of r=b-Ax)
1511  for (unsigned i = 0; i < Nlevel; i++)
1512  {
1513  // If we're allowed to output
1514  if (!Suppress_all_output)
1515  {
1516  // Output the level we're working on
1517  oomph_info << "Setting up MG structures on level: " << i << "\n"
1518  << std::endl;
1519  }
1520 
1521  // Resize the solution and RHS vector
1522  unsigned n_dof_per_block = Mg_hierarchy_pt[i]->ndof() / 2;
1523 
1524  // Create the linear algebra distribution to build the vectors
1525  LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
1526  Mg_hierarchy_pt[i]->communicator_pt(), n_dof_per_block, false);
1527 
1528 #ifdef PARANOID
1529 #ifdef OOMPH_HAS_MPI
1530  // Set up the warning messages
1531  std::string warning_message = "Setup of distribution has not been ";
1532  warning_message += "tested with MPI.";
1533 
1534  // If we're not running the code in serial
1535  if (dist_pt->communicator_pt()->nproc() > 1)
1536  {
1537  // Throw a warning
1538  OomphLibWarning(
1539  warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1540  }
1541 #endif
1542 #endif
1543 
1544  // Create storage for the i-th level system matrix:
1545  //-------------------------------------------------
1546  // Resize the i-th entry of the matrix storage vector to contain two
1547  // CRDoubleMatrix pointers
1548  Mg_matrices_storage_pt[i].resize(2, 0);
1549 
1550  // Loop over the real and imaginary part
1551  for (unsigned j = 0; j < 2; j++)
1552  {
1553  // Dynamically allocate a new CRDoubleMatrix
1554  Mg_matrices_storage_pt[i][j] = new CRDoubleMatrix;
1555  }
1556 
1557  // Build the matrix distribution (real part)
1558  Mg_matrices_storage_pt[i][0]->build(dist_pt);
1559 
1560  // Build the matrix distribution (imaginary part)
1561  Mg_matrices_storage_pt[i][1]->build(dist_pt);
1562 
1563  // Create storage for the i-th level solution vector:
1564  //---------------------------------------------------
1565  // Resize the i-th level solution vector to contain two DoubleVectors
1566  X_mg_vectors_storage[i].resize(2);
1567 
1568  // Build the approximate solution (real part)
1569  X_mg_vectors_storage[i][0].build(dist_pt);
1570 
1571  // Build the approximate solution (imaginary part)
1572  X_mg_vectors_storage[i][1].build(dist_pt);
1573 
1574  // Create storage for the i-th level RHS vector:
1575  //----------------------------------------------
1576  // Resize the i-th level RHS vector to contain two DoubleVectors
1577  Rhs_mg_vectors_storage[i].resize(2);
1578 
1579  // Build the point source function (real part)
1580  Rhs_mg_vectors_storage[i][0].build(dist_pt);
1581 
1582  // Build the point source function (imaginary part)
1583  Rhs_mg_vectors_storage[i][1].build(dist_pt);
1584 
1585  // Create storage for the i-th level residual vector:
1586  //---------------------------------------------------
1587  // Resize the i-th level residual vector to contain two DoubleVectors
1588  Residual_mg_vectors_storage[i].resize(2);
1589 
1590  // Build the residual vector, r=b-Ax (real part)
1591  Residual_mg_vectors_storage[i][0].build(dist_pt);
1592 
1593  // Build the residual vector, r=b-Ax (imaginary part)
1594  Residual_mg_vectors_storage[i][1].build(dist_pt);
1595 
1596  // Delete the distribution pointer
1597  delete dist_pt;
1598 
1599  // Make it a null pointer
1600  dist_pt = 0;
1601 
1602  // Compute system matrix on the current level. On the finest level of the
1603  // hierarchy the system matrix is given by the complex-shifted Laplacian
1604  // preconditioner. On the subsequent levels the Galerkin approximation
1605  // is used to give us a coarse-grid representation of the problem
1606  if (i == 0)
1607  {
1608  // Initialise the timer start variable
1609  double t_jac_start = 0.0;
1610 
1611  // If we're allowed to output things
1612  if (!Suppress_all_output)
1613  {
1614  // Start timer for Jacobian setup
1615  t_jac_start = TimingHelpers::timer();
1616  }
1617 
1618 #ifdef PARANOID
1619  // Make sure the block preconditioner returns the correct sort of matrix
1620  block_preconditioner_self_test();
1621 #endif
1622 
1623  // The preconditioner works with one mesh; set it!
1624  this->set_nmesh(1);
1625 
1626  // Elements in actual pml layer are trivially wrapped versions of their
1627  // bulk counterparts. Technically they are different elements so we have
1628  // to allow different element types
1629  bool allow_different_element_types_in_mesh = true;
1630  this->set_mesh(0,
1631  Mg_hierarchy_pt[0]->mesh_pt(),
1632  allow_different_element_types_in_mesh);
1633 
1634 #ifdef PARANOID
1635  // Find the number of dof types we're dealing with
1636  unsigned n_dof_types = this->ndof_types();
1637 
1638  // This preconditioner only works for 2 dof types
1639  if (n_dof_types != 2)
1640  {
1641  // Create an error message
1642  std::stringstream tmp;
1643  tmp
1644  << "This preconditioner only works for problems with 2 dof types\n"
1645  << "Yours has " << n_dof_types;
1646 
1647  // Throw an error
1648  throw OomphLibError(
1649  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1650  }
1651 #endif
1652 
1653  // If we're not using a value of zero for the alpha shift
1654  if (Alpha_shift != 0.0)
1655  {
1656  //------------------------------------------------------------------
1657  // Set the damping in all of the PML elements to create the complex-
1658  // shifted Laplacian preconditioner:
1659  //------------------------------------------------------------------
1660  // Create a pointer which will contain the value of the shift given
1661  // by Alpha_shift
1662  double* alpha_shift_pt = new double(Alpha_shift);
1663 
1664  // Calculate the number of elements in the mesh
1665  unsigned n_element = Mg_hierarchy_pt[0]->mesh_pt()->nelement();
1666 
1667  // To grab the static variable used to set the value of alpha we first
1668  // need to get an element of type PMLHelmholtzEquations (we
1669  // arbitrarily chose the first element in the mesh)
1670  PMLHelmholtzEquations<DIM>* el_pt =
1671  dynamic_cast<PMLHelmholtzEquations<DIM>*>(
1672  Mg_hierarchy_pt[0]->mesh_pt()->element_pt(0));
1673 
1674  // Now grab the pointer from the element
1675  static double* default_physical_constant_value_pt = el_pt->alpha_pt();
1676 
1677  // Loop over all of the elements
1678  for (unsigned i_el = 0; i_el < n_element; i_el++)
1679  {
1680  // Upcast from a GeneralisedElement to a PmlHelmholtzElement
1681  PMLHelmholtzEquations<DIM>* el_pt =
1682  dynamic_cast<PMLHelmholtzEquations<DIM>*>(
1683  Mg_hierarchy_pt[0]->mesh_pt()->element_pt(i_el));
1684 
1685  // Make the internal element alpha pointer point to Alpha_shift (the
1686  // chosen value of the shift to be applied to the problem)
1687  el_pt->alpha_pt() = alpha_shift_pt;
1688  }
1689 
1690  //---------------------------------------------------------------------
1691  // We want to solve the problem with the modified Jacobian but the
1692  // Preconditioner will have a handle to pointer which points to the
1693  // system matrix. If we wish to use the block preconditioner to
1694  // extract the appropriate blocks of the matrix we need to assign it.
1695  // To avoid losing the original system matrix we will create a
1696  // temporary pointer to it which will be reassigned after we have use
1697  // the block preconditioner to extract the blocks of the shifted
1698  // matrix:
1699  //---------------------------------------------------------------------
1700  // Create a temporary pointer to the Jacobian
1701  CRDoubleMatrix* jacobian_pt = this->matrix_pt();
1702 
1703  // Create a new CRDoubleMatrix to hold the shifted Jacobian matrix
1704  CRDoubleMatrix* shifted_jacobian_pt = new CRDoubleMatrix;
1705 
1706  // Allocate space for the residuals
1707  DoubleVector residuals;
1708 
1709  // Get the residuals vector and the shifted Jacobian. Note, we do
1710  // not need to assign the residuals vector; we're simply using
1711  // MG as a preconditioner
1712  Mg_hierarchy_pt[0]->get_jacobian(residuals, *shifted_jacobian_pt);
1713 
1714  // Replace the current matrix used in Preconditioner by the new matrix
1715  this->set_matrix_pt(shifted_jacobian_pt);
1716 
1717  // Set up the generic block look up scheme
1718  this->block_setup();
1719 
1720  // Extract the number of blocks.
1721  unsigned nblock_types = this->nblock_types();
1722 
1723 #ifdef PARANOID
1724  // PARANOID check - there must only be two block types
1725  if (nblock_types != 2)
1726  {
1727  // Create the error message
1728  std::stringstream tmp;
1729  tmp << "There are supposed to be two block types.\nYours has "
1730  << nblock_types << std::endl;
1731 
1732  // Throw an error
1733  throw OomphLibError(
1734  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1735  }
1736 #endif
1737 
1738  // Store the level
1739  unsigned level = 0;
1740 
1741  // Loop over the rows of the block matrix
1742  for (unsigned i_row = 0; i_row < nblock_types; i_row++)
1743  {
1744  // Fix the column index
1745  unsigned j_col = 0;
1746 
1747  // Extract the required blocks, i.e. the first column
1748  this->get_block(
1749  i_row, j_col, *Mg_matrices_storage_pt[level][i_row]);
1750  }
1751 
1752  // The blocks have been extracted from the shifted Jacobian therefore
1753  // we no longer have any use for it
1754  delete shifted_jacobian_pt;
1755 
1756  // Now the appropriate blocks have been extracted from the shifted
1757  // Jacobian we reset the matrix pointer in Preconditioner to the
1758  // original Jacobian so the linear solver isn't affected
1759  this->set_matrix_pt(jacobian_pt);
1760 
1761  //--------------------------------------------------------
1762  // Reassign the damping factor in all of the PML elements:
1763  //--------------------------------------------------------
1764  // Loop over all of the elements
1765  for (unsigned i_el = 0; i_el < n_element; i_el++)
1766  {
1767  // Upcast from a GeneralisedElement to a PmlHelmholtzElement
1768  PMLHelmholtzEquations<DIM>* el_pt =
1769  dynamic_cast<PMLHelmholtzEquations<DIM>*>(
1770  Mg_hierarchy_pt[0]->mesh_pt()->element_pt(i_el));
1771 
1772  // Set the value of alpha
1773  el_pt->alpha_pt() = default_physical_constant_value_pt;
1774  }
1775 
1776  // We've finished using the alpha_shift_pt pointer so delete it
1777  // as it was dynamically allocated
1778  delete alpha_shift_pt;
1779 
1780  // Make it a null pointer
1781  alpha_shift_pt = 0;
1782  }
1783  // If the value of the shift is zero then we use the original
1784  // Jacobian matrix
1785  else
1786  {
1787  // The Jacobian has already been provided so now we just need to set
1788  // up the generic block look up scheme
1789  this->block_setup();
1790 
1791  // Extract the number of blocks.
1792  unsigned nblock_types = this->nblock_types();
1793 
1794 #ifdef PARANOID
1795  // If there are not only two block types we have a problem
1796  if (nblock_types != 2)
1797  {
1798  std::stringstream tmp;
1799  tmp << "There are supposed to be two block types.\n"
1800  << "Yours has " << nblock_types;
1801  throw OomphLibError(
1802  tmp.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1803  }
1804 #endif
1805 
1806  // Store the level
1807  unsigned level = 0;
1808 
1809  // Loop over the rows of the block matrix
1810  for (unsigned i_row = 0; i_row < nblock_types; i_row++)
1811  {
1812  // Fix the column index (since we only want the first column)
1813  unsigned j_col = 0;
1814 
1815  // Extract the required blocks
1816  this->get_block(
1817  i_row, j_col, *Mg_matrices_storage_pt[level][i_row]);
1818  }
1819  } // if (Alpha_shift!=0.0) else
1820 
1821  if (!Suppress_all_output)
1822  {
1823  // Document the time taken
1824  double t_jac_end = TimingHelpers::timer();
1825  double jacobian_setup_time = t_jac_end - t_jac_start;
1826  oomph_info << " - Time for setup of Jacobian block matrix [sec]: "
1827  << jacobian_setup_time << "\n"
1828  << std::endl;
1829  }
1830  }
1831  // If we're not dealing with the finest level we're dealing with a
1832  // coarse-grid problem
1833  else
1834  {
1835  // Initialise the timer start variable
1836  double t_gal_start = 0.0;
1837 
1838  // If we're allowed
1839  if (!Suppress_all_output)
1840  {
1841  // Start timer for Galerkin matrix calculation
1842  t_gal_start = TimingHelpers::timer();
1843  }
1844 
1845  //---------------------------------------------------------------------
1846  // The system matrix on the coarser levels must be formed using the
1847  // Galerkin approximation which we do by calculating the product
1848  // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
1849  // finer grid system matrix is formed by multiplying by the (fine grid)
1850  // restriction matrix from the left and the (fine grid) interpolation
1851  // matrix from the left. Fortunately, since the restriction and
1852  // interpolation acts on the real and imaginary parts separately we
1853  // can calculate the real and imaginary parts of the Galerkin
1854  // approximation separately.
1855  //---------------------------------------------------------------------
1856 
1857  //----------------------------------------------
1858  // Real component of the Galerkin approximation:
1859  //----------------------------------------------
1860  // First we need to calculate A^h * I^h_2h which we store as A^2h
1861  Mg_matrices_storage_pt[i - 1][0]->multiply(
1862  *Interpolation_matrices_storage_pt[i - 1],
1863  *Mg_matrices_storage_pt[i][0]);
1864 
1865  // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1866  // was just calculated. This updates A^2h to give us the true
1867  // Galerkin approximation to the finer grid matrix
1868  Restriction_matrices_storage_pt[i - 1]->multiply(
1869  *Mg_matrices_storage_pt[i][0], *Mg_matrices_storage_pt[i][0]);
1870 
1871  //---------------------------------------------------
1872  // Imaginary component of the Galerkin approximation:
1873  //---------------------------------------------------
1874  // First we need to calculate A^h * I^h_2h which we store as A^2h
1875  Mg_matrices_storage_pt[i - 1][1]->multiply(
1876  *Interpolation_matrices_storage_pt[i - 1],
1877  *Mg_matrices_storage_pt[i][1]);
1878 
1879  // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1880  // was just calculated. This updates A^2h to give us the true
1881  // Galerkin approximation to the finer grid matrix
1882  Restriction_matrices_storage_pt[i - 1]->multiply(
1883  *Mg_matrices_storage_pt[i][1], *Mg_matrices_storage_pt[i][1]);
1884 
1885  // If the user did not choose to suppress everything
1886  if (!Suppress_all_output)
1887  {
1888  // End timer for Galerkin matrix calculation
1889  double t_gal_end = TimingHelpers::timer();
1890 
1891  // Calculate setup time
1892  double galerkin_matrix_calculation_time = t_gal_end - t_gal_start;
1893 
1894  // Document the time taken
1895  oomph_info << " - Time for system block matrix formation using the "
1896  << "Galerkin approximation [sec]: "
1897  << galerkin_matrix_calculation_time << "\n"
1898  << std::endl;
1899  }
1900  } // if (i==0) else
1901 
1902  // We only
1903  if (i < Nlevel - 1)
1904  {
1905  // Find the maximum edge width, h:
1906  //--------------------------------
1907  // Initialise the timer start variable
1908  double t_para_start = 0.0;
1909 
1910  // If we're allowed to output things
1911  if (!Suppress_all_output)
1912  {
1913  // Start timer for parameter calculation on the i-th level
1914  t_para_start = TimingHelpers::timer();
1915  }
1916 
1917  // Calculate the i-th entry of Wavenumber and Max_edge_width
1918  maximum_edge_width(i, Max_edge_width[i]);
1919 
1920  // If the user did not choose to suppress everything
1921  if (!Suppress_all_output)
1922  {
1923  // End timer for Galerkin matrix calculation
1924  double t_para_end = TimingHelpers::timer();
1925 
1926  // Calculate setup time
1927  double parameter_calculation_time = t_para_end - t_para_start;
1928 
1929  // Document the time taken
1930  oomph_info << " - Time for maximum edge width calculation [sec]: "
1931  << parameter_calculation_time << "\n"
1932  << std::endl;
1933  }
1934  } // if (i<Nlevel-1)
1935  } // for (unsigned i=0;i<Nlevel;i++)
1936 
1937  // The last system matrix that needs to be setup is the fully expanded
1938  // version of the system matrix on the coarsest level. This is needed
1939  // for the direct solve on the coarsest level
1940  setup_coarsest_level_structures();
1941 
1942  // If we're allowed to output
1943  if (!Suppress_all_output)
1944  {
1945  // Stop clock
1946  double t_m_end = TimingHelpers::timer();
1947  double total_setup_time = double(t_m_end - t_m_start);
1948  oomph_info << "CPU time for setup of MG structures [sec]: "
1949  << total_setup_time << std::endl;
1950 
1951  // Notify user that the hierarchy of levels is complete
1952  oomph_info << "\n============"
1953  << "Multigrid Structures Setup Complete"
1954  << "===========\n"
1955  << std::endl;
1956  }
1957  } // End of setup_mg_structures
oomph::HelmholtzMGPreconditioner::setup_coarsest_level_structures
void setup_coarsest_level_structures()
Definition: helmholtz_geometric_multigrid.h:1980
oomph::HelmholtzMGPreconditioner::block_preconditioner_self_test
void block_preconditioner_self_test()
Definition: helmholtz_geometric_multigrid.h:873
oomph::HelmholtzMGPreconditioner::maximum_edge_width
void maximum_edge_width(const unsigned &level, double &h)
oomph::HelmholtzMGPreconditioner::Wavenumber
double Wavenumber
The value of the wavenumber, k.
Definition: helmholtz_geometric_multigrid.h:680

References oomph::PMLHelmholtzEquations< DIM >::alpha_pt(), GlobalParameters::Alpha_shift, oomph::CRDoubleMatrix::build(), oomph::LinearAlgebraDistribution::communicator_pt(), i, j, oomph::PMLHelmholtzEquations< DIM >::k_squared(), oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, sqrt(), oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), tmp, and GlobalParameters::Wavenumber.

◆ setup_smoothers()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_smoothers
private

Set up the smoothers on all levels.

Function to set up all of the smoothers once the system matrices have been set up

2542  {
2543  // Initialise the timer start variable
2544  double t_m_start = 0.0;
2545 
2546  // Start the clock (if we're allowed to time things)
2547  if (!Suppress_all_output)
2548  {
2549  // Notify user
2550  oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
2551 
2552  // Start the clock
2553  t_m_start = TimingHelpers::timer();
2554  }
2555 
2556  // Loop over the levels and assign the pre- and post-smoother pointers
2557  for (unsigned i = 0; i < Nlevel - 1; i++)
2558  {
2559  // If the pre-smoother factory function pointer hasn't been assigned
2560  // then we simply create a new instance of the ComplexDampedJacobi
2561  // smoother which is the default pre-smoother
2562  if (0 == Pre_smoother_factory_function_pt)
2563  {
2564  // If the value of kh is strictly less than 0.5 then use damped Jacobi
2565  // as a smoother on this level otherwise use complex GMRES as a smoother
2566  // [see Elman et al."A multigrid method enhanced by Krylov subspace
2567  // iteration for discrete Helmholtz equations", p.1305]
2568  if (Wavenumber * Max_edge_width[i] < 0.5)
2569  {
2570  // Make a new ComplexDampedJacobi object and assign its pointer
2571  ComplexDampedJacobi<CRDoubleMatrix>* damped_jacobi_pt =
2572  new ComplexDampedJacobi<CRDoubleMatrix>;
2573 
2574  // Assign the pre-smoother pointer
2575  Pre_smoothers_storage_pt[i] = damped_jacobi_pt;
2576 
2577  // Make the smoother calculate the value of omega and store it
2578  damped_jacobi_pt->calculate_omega(Wavenumber, Max_edge_width[i]);
2579  }
2580  else
2581  {
2582  // Make a new ComplexGMRES object and assign its pointer
2583  ComplexGMRES<CRDoubleMatrix>* gmres_pt =
2584  new ComplexGMRES<CRDoubleMatrix>;
2585 
2586  // Assign the pre-smoother pointer
2587  Pre_smoothers_storage_pt[i] = gmres_pt;
2588  }
2589  }
2590  // Otherwise we use the pre-smoother factory function pointer to
2591  // generate a new pre-smoother
2592  else
2593  {
2594  // Get a pointer to an object of the same type as the pre-smoother
2595  Pre_smoothers_storage_pt[i] = (*Pre_smoother_factory_function_pt)();
2596  } // if (0==Pre_smoother_factory_function_pt)
2597 
2598  // If the post-smoother factory function pointer hasn't been assigned
2599  // then we simply create a new instance of the ComplexDampedJacobi
2600  // smoother which is the default post-smoother
2601  if (0 == Post_smoother_factory_function_pt)
2602  {
2603  // If the value of kh is strictly less than 0.52 then use damped Jacobi
2604  // as a smoother on this level otherwise use complex GMRES as a smoother
2605  // [see Elman et al."A multigrid method enhanced by Krylov subspace
2606  // iteration for discrete Helmholtz equations", p.1295]
2607  if (Wavenumber * Max_edge_width[i] < 0.5)
2608  {
2609  // Make a new ComplexDampedJacobi object and assign its pointer
2610  ComplexDampedJacobi<CRDoubleMatrix>* damped_jacobi_pt =
2611  new ComplexDampedJacobi<CRDoubleMatrix>;
2612 
2613  // Assign the pre-smoother pointer
2614  Post_smoothers_storage_pt[i] = damped_jacobi_pt;
2615 
2616  // Make the smoother calculate the value of omega and store it
2617  damped_jacobi_pt->calculate_omega(Wavenumber, Max_edge_width[i]);
2618  }
2619  else
2620  {
2621  // Make a new ComplexGMRES object and assign its pointer
2622  ComplexGMRES<CRDoubleMatrix>* gmres_pt =
2623  new ComplexGMRES<CRDoubleMatrix>;
2624 
2625  // Assign the pre-smoother pointer
2626  Post_smoothers_storage_pt[i] = gmres_pt;
2627  }
2628  }
2629  // Otherwise we use the post-smoother factory function pointer to
2630  // generate a new post-smoother
2631  else
2632  {
2633  // Get a pointer to an object of the same type as the post-smoother
2634  Post_smoothers_storage_pt[i] = (*Post_smoother_factory_function_pt)();
2635  }
2636  } // if (0==Post_smoother_factory_function_pt)
2637 
2638  // Set the tolerance for the pre- and post-smoothers. The norm of the
2639  // solution which is compared against the tolerance is calculated
2640  // differently to the multigrid solver. To ensure that the smoother
2641  // continues to solve for the prescribed number of iterations we
2642  // lower the tolerance
2643  for (unsigned i = 0; i < Nlevel - 1; i++)
2644  {
2645  // Set the tolerance of the i-th level pre-smoother
2646  Pre_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
2647 
2648  // Set the tolerance of the i-th level post-smoother
2649  Post_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
2650  }
2651 
2652  // Set the number of pre- and post-smoothing iterations in each
2653  // pre- and post-smoother
2654  for (unsigned i = 0; i < Nlevel - 1; i++)
2655  {
2656  // Set the number of pre-smoothing iterations as the value of Npre_smooth
2657  Pre_smoothers_storage_pt[i]->max_iter() = Npre_smooth;
2658 
2659  // Set the number of pre-smoothing iterations as the value of Npost_smooth
2660  Post_smoothers_storage_pt[i]->max_iter() = Npost_smooth;
2661  }
2662 
2663  // Complete the setup of all of the smoothers
2664  for (unsigned i = 0; i < Nlevel - 1; i++)
2665  {
2666  // Pass a pointer to the system matrix on the i-th level to the i-th
2667  // level pre-smoother
2668  Pre_smoothers_storage_pt[i]->complex_smoother_setup(
2669  Mg_matrices_storage_pt[i]);
2670 
2671  // Pass a pointer to the system matrix on the i-th level to the i-th
2672  // level post-smoother
2673  Post_smoothers_storage_pt[i]->complex_smoother_setup(
2674  Mg_matrices_storage_pt[i]);
2675  }
2676 
2677  // Set up the distributions of each smoother
2678  for (unsigned i = 0; i < Nlevel - 1; i++)
2679  {
2680  // Get the number of dofs on the i-th level of the hierarchy.
2681  // This will be the same as the size of the solution vector
2682  // associated with the i-th level
2683  unsigned n_dof = X_mg_vectors_storage[i][0].nrow();
2684 
2685  // Create a LinearAlgebraDistribution which will be passed to the
2686  // linear solver
2687  LinearAlgebraDistribution dist(
2688  Mg_hierarchy_pt[i]->communicator_pt(), n_dof, false);
2689 
2690 #ifdef PARANOID
2691 #ifdef OOMPH_HAS_MPI
2692  // Set up the warning messages
2693  std::string warning_message =
2694  "Setup of pre- and post-smoother distribution ";
2695  warning_message += "has not been tested with MPI.";
2696 
2697  // If we're not running the code in serial
2698  if (dist.communicator_pt()->nproc() > 1)
2699  {
2700  // Throw a warning
2701  OomphLibWarning(
2702  warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2703  }
2704 #endif
2705 #endif
2706 
2707  // Build the distribution of the pre-smoother
2708  Pre_smoothers_storage_pt[i]->build_distribution(dist);
2709 
2710  // Build the distribution of the post-smoother
2711  Post_smoothers_storage_pt[i]->build_distribution(dist);
2712  }
2713 
2714  // Disable the smoother output
2715  if (!Doc_time)
2716  {
2717  // Disable time documentation from the smoothers and the SuperLU linear
2718  // solver (this latter is only done on the coarsest level where it is
2719  // required)
2720  disable_smoother_and_superlu_doc_time();
2721  }
2722 
2723  // If we're allowed to output
2724  if (!Suppress_all_output)
2725  {
2726  // Stop clock
2727  double t_m_end = TimingHelpers::timer();
2728  double total_setup_time = double(t_m_end - t_m_start);
2729  oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
2730  << total_setup_time << std::endl;
2731 
2732  // Notify user that the extraction is complete
2733  oomph_info << "\n=================="
2734  << "Smoother Setup Complete"
2735  << "=================" << std::endl;
2736  }
2737  } // End of setup_smoothers
oomph::HelmholtzMGPreconditioner::disable_smoother_and_superlu_doc_time
void disable_smoother_and_superlu_doc_time()
Suppress the output of both smoothers and SuperLU.
Definition: helmholtz_geometric_multigrid.h:287

References oomph::ComplexDampedJacobi< MATRIX >::calculate_omega(), oomph::LinearAlgebraDistribution::communicator_pt(), i, oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Global_string_for_annotation::string(), oomph::TimingHelpers::timer(), oomph::IterativeLinearSolver::tolerance(), and GlobalParameters::Wavenumber.

◆ setup_transfer_matrices()

template<unsigned DIM>
void oomph::HelmholtzMGPreconditioner< DIM >::setup_transfer_matrices

Setup the transfer matrices on each level.

Set up the transfer matrices. Both the pure injection and full weighting method have been implemented here but it is highly recommended that full weighting is used in general. In both methods the transpose of the transfer matrix is used to transfer a vector back

1423  {
1424  // Initialise the timer start variable
1425  double t_r_start = 0.0;
1426 
1427  // Notify the user (if we're allowed)
1428  if (!Suppress_all_output)
1429  {
1430  // Notify user of progress
1431  oomph_info << "Creating the transfer matrices." << std::endl;
1432 
1433  // Start the clock!
1434  t_r_start = TimingHelpers::timer();
1435  }
1436 
1437  // Using full weighting so use setup_interpolation_matrices.
1438  // Note: There are two methods to choose from here, the ideal choice is
1439  // setup_interpolation_matrices() but that requires a refineable mesh base
1440  if (dynamic_cast<TreeBasedRefineableMeshBase*>(
1441  Mg_problem_pt->mg_bulk_mesh_pt()))
1442  {
1443  setup_interpolation_matrices();
1444  }
1445  // If the mesh is unstructured we have to use the locate_zeta function
1446  // to set up the interpolation matrices
1447  else
1448  {
1449  setup_interpolation_matrices_unstructured();
1450  }
1451 
1452  // Loop over all levels that will be assigned a restriction matrix
1453  set_restriction_matrices_as_interpolation_transposes();
1454 
1455  // If we're allowed
1456  if (!Suppress_all_output)
1457  {
1458  // Stop the clock
1459  double t_r_end = TimingHelpers::timer();
1460  double total_G_setup_time = double(t_r_end - t_r_start);
1461  oomph_info << "\nCPU time for transfer matrices setup [sec]: "
1462  << total_G_setup_time << std::endl;
1463 
1464  // Notify user that the hierarchy of levels is complete
1465  oomph_info
1466  << "\n============Transfer Matrices Setup Complete=============="
1467  << "\n"
1468  << std::endl;
1469  }
1470  } // End of setup_transfer_matrices function
oomph::HelmholtzMGPreconditioner::setup_interpolation_matrices
void setup_interpolation_matrices()
Setup the interpolation matrix on each level.
Definition: helmholtz_geometric_multigrid.h:2744
oomph::HelmholtzMGPreconditioner::set_restriction_matrices_as_interpolation_transposes
void set_restriction_matrices_as_interpolation_transposes()
Definition: helmholtz_geometric_multigrid.h:461
oomph::HelmholtzMGPreconditioner::setup_interpolation_matrices_unstructured
void setup_interpolation_matrices_unstructured()
Setup the interpolation matrices.
Definition: helmholtz_geometric_multigrid.h:3243

References oomph::oomph_info, and oomph::TimingHelpers::timer().

◆ tolerance()

template<unsigned DIM>
double& oomph::HelmholtzMGPreconditioner< DIM >::tolerance ( )
inline

Access function for the variable Tolerance (lvalue)

205  {
206  // Return the variable, Tolerance
207  return Tolerance;
208  } // End of tolerance

References oomph::HelmholtzMGPreconditioner< DIM >::Tolerance.

Member Data Documentation

◆ Alpha_shift

template<unsigned DIM>
double oomph::HelmholtzMGPreconditioner< DIM >::Alpha_shift
private

Temporary version of the shift – to run bash scripts.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::alpha_shift().

◆ Coarsest_matrix_mg_pt

template<unsigned DIM>
CRDoubleMatrix* oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_matrix_mg_pt
private

Stores the system matrix on the coarest level in the fully expanded format: |--—|---—| | A_r | -A_c | A = |--—|---—|. | A_c | A_r | |--—|---—| Note: this is NOT the same as A = A_r + iA_c

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), oomph::HelmholtzMGPreconditioner< DIM >::direct_solve(), and oomph::HelmholtzMGPreconditioner< DIM >::disable_smoother_and_superlu_doc_time().

◆ Coarsest_rhs_mg

template<unsigned DIM>
DoubleVector oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_rhs_mg
private

Assuming we're solving the system Ax=b, this vector will contain the expanded solution vector on the coarsest level of the heirarchy. This will have the form: |--—| | b_r | b = |--—|. | b_c | |--—|

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::direct_solve().

◆ Coarsest_x_mg

template<unsigned DIM>
DoubleVector oomph::HelmholtzMGPreconditioner< DIM >::Coarsest_x_mg
private

Assuming we're solving the system Ax=b, this vector will contain the expanded solution vector on the coarsest level of the heirarchy. This will have the form: |--—| | x_r | x = |--—|. | x_c | |--—|

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::direct_solve().

◆ Doc_time

template<unsigned DIM>
bool oomph::HelmholtzMGPreconditioner< DIM >::Doc_time
private

Indicates whether or not time documentation should be used.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::disable_doc_time(), oomph::HelmholtzMGPreconditioner< DIM >::disable_output(), oomph::HelmholtzMGPreconditioner< DIM >::disable_v_cycle_output(), oomph::HelmholtzMGPreconditioner< DIM >::enable_doc_time(), oomph::HelmholtzMGPreconditioner< DIM >::enable_output(), and oomph::HelmholtzMGPreconditioner< DIM >::enable_v_cycle_output().

◆ Has_been_setup

template<unsigned DIM>
bool oomph::HelmholtzMGPreconditioner< DIM >::Has_been_setup
private

Boolean variable to indicate whether or not the solver has been setup.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory().

◆ Has_been_solved

template<unsigned DIM>
bool oomph::HelmholtzMGPreconditioner< DIM >::Has_been_solved
private

Boolean variable to indicate whether or not the problem was successfully solved

◆ Interpolation_matrices_storage_pt

template<unsigned DIM>
Vector<CRDoubleMatrix*> oomph::HelmholtzMGPreconditioner< DIM >::Interpolation_matrices_storage_pt
private

Vector to store the interpolation matrices.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), oomph::HelmholtzMGPreconditioner< DIM >::interpolation_matrix_set(), and oomph::HelmholtzMGPreconditioner< DIM >::set_restriction_matrices_as_interpolation_transposes().

◆ Max_edge_width

template<unsigned DIM>
Vector<double> oomph::HelmholtzMGPreconditioner< DIM >::Max_edge_width
private

Vector to storage the maximum edge width of each mesh. We only need the maximum edge width on levels where we use a smoother to determine the value of kh

◆ Mg_hierarchy_pt

template<unsigned DIM>
Vector<HelmholtzMGProblem*> oomph::HelmholtzMGPreconditioner< DIM >::Mg_hierarchy_pt
private

Vector containing pointers to problems in hierarchy.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::HelmholtzMGPreconditioner(), and oomph::HelmholtzMGPreconditioner< DIM >::interpolation_matrix_set().

◆ Mg_matrices_storage_pt

template<unsigned DIM>
Vector<Vector<CRDoubleMatrix*> > oomph::HelmholtzMGPreconditioner< DIM >::Mg_matrices_storage_pt
private

Vector of vectors to store the system matrices. The i-th entry in this vector contains a vector of length two. The first entry of which contains the real part of the system matrix which we refer to as A_r and the second entry contains the imaginary part of the system matrix which we refer to as A_c. That is to say, the true system matrix is given by A = A_r + iA_c

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory().

◆ Mg_problem_pt

template<unsigned DIM>
HelmholtzMGProblem* oomph::HelmholtzMGPreconditioner< DIM >::Mg_problem_pt
private

Pointer to the MG problem (deep copy)

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::HelmholtzMGPreconditioner().

◆ Nlevel

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::Nlevel
private

The number of levels in the multigrid heirachy.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), oomph::HelmholtzMGPreconditioner< DIM >::direct_solve(), oomph::HelmholtzMGPreconditioner< DIM >::disable_smoother_and_superlu_doc_time(), and oomph::HelmholtzMGPreconditioner< DIM >::set_restriction_matrices_as_interpolation_transposes().

◆ Npost_smooth

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::Npost_smooth
private

Number of post-smoothing steps.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::npost_smooth().

◆ Npre_smooth

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::Npre_smooth
private

Number of pre-smoothing steps.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::npre_smooth().

◆ Nvcycle

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::Nvcycle
private

Maximum number of V-cycles.

◆ Post_smoother_factory_function_pt

template<unsigned DIM>
PostSmootherFactoryFctPt oomph::HelmholtzMGPreconditioner< DIM >::Post_smoother_factory_function_pt
private

Function to create post-smoothers.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::set_post_smoother_factory_function().

◆ Post_smoothers_storage_pt

template<unsigned DIM>
Vector<HelmholtzSmoother*> oomph::HelmholtzMGPreconditioner< DIM >::Post_smoothers_storage_pt
private

Vector to store the post-smoothers.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), oomph::HelmholtzMGPreconditioner< DIM >::disable_smoother_and_superlu_doc_time(), and oomph::HelmholtzMGPreconditioner< DIM >::post_smooth().

◆ Pre_smoother_factory_function_pt

template<unsigned DIM>
PreSmootherFactoryFctPt oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoother_factory_function_pt
private

Function to create pre-smoothers.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::set_pre_smoother_factory_function().

◆ Pre_smoothers_storage_pt

template<unsigned DIM>
Vector<HelmholtzSmoother*> oomph::HelmholtzMGPreconditioner< DIM >::Pre_smoothers_storage_pt
private

Vector to store the pre-smoothers.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), oomph::HelmholtzMGPreconditioner< DIM >::disable_smoother_and_superlu_doc_time(), and oomph::HelmholtzMGPreconditioner< DIM >::pre_smooth().

◆ Residual_mg_vectors_storage

template<unsigned DIM>
Vector<Vector<DoubleVector> > oomph::HelmholtzMGPreconditioner< DIM >::Residual_mg_vectors_storage
private

Vector to vectors to store the residual vectors. This uses the same format as the X_mg_vectors_storage vector

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::residual_norm().

◆ Restriction_matrices_storage_pt

template<unsigned DIM>
Vector<CRDoubleMatrix*> oomph::HelmholtzMGPreconditioner< DIM >::Restriction_matrices_storage_pt
private

Vector to store the restriction matrices.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::clean_up_memory(), and oomph::HelmholtzMGPreconditioner< DIM >::set_restriction_matrices_as_interpolation_transposes().

◆ Rhs_mg_vectors_storage

template<unsigned DIM>
Vector<Vector<DoubleVector> > oomph::HelmholtzMGPreconditioner< DIM >::Rhs_mg_vectors_storage
private

Vector of vectors to store the RHS vectors. This uses the same format as the X_mg_vectors_storage vector

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::direct_solve(), oomph::HelmholtzMGPreconditioner< DIM >::post_smooth(), oomph::HelmholtzMGPreconditioner< DIM >::pre_smooth(), and oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve().

◆ Stream_pt

template<unsigned DIM>
std::ostream* oomph::HelmholtzMGPreconditioner< DIM >::Stream_pt
private

Pointer to the output stream – defaults to std::cout.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::disable_output(), oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve(), and oomph::HelmholtzMGPreconditioner< DIM >::setup().

◆ Suppress_all_output

template<unsigned DIM>
bool oomph::HelmholtzMGPreconditioner< DIM >::Suppress_all_output
private

If this is set to true then all output from the solver is suppressed.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::disable_output(), oomph::HelmholtzMGPreconditioner< DIM >::enable_output(), oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve(), and oomph::HelmholtzMGPreconditioner< DIM >::setup().

◆ Suppress_v_cycle_output

template<unsigned DIM>
bool oomph::HelmholtzMGPreconditioner< DIM >::Suppress_v_cycle_output
private

Indicates whether or not the V-cycle output should be suppressed.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::disable_output(), oomph::HelmholtzMGPreconditioner< DIM >::disable_v_cycle_output(), oomph::HelmholtzMGPreconditioner< DIM >::enable_output(), oomph::HelmholtzMGPreconditioner< DIM >::enable_v_cycle_output(), and oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve().

◆ Tolerance

template<unsigned DIM>
double oomph::HelmholtzMGPreconditioner< DIM >::Tolerance
private

The tolerance of the multigrid preconditioner.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::tolerance().

◆ V_cycle_counter

template<unsigned DIM>
unsigned oomph::HelmholtzMGPreconditioner< DIM >::V_cycle_counter
private

Pointer to counter for V-cycles.

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::iterations().

◆ Wavenumber

template<unsigned DIM>
double oomph::HelmholtzMGPreconditioner< DIM >::Wavenumber
private

The value of the wavenumber, k.

◆ X_mg_vectors_storage

template<unsigned DIM>
Vector<Vector<DoubleVector> > oomph::HelmholtzMGPreconditioner< DIM >::X_mg_vectors_storage
private

Vector of vectors to store the solution vectors (X_mg) in two parts; the real and imaginary. To access the real part of the solution vector on the i-th level we need to access X_mg_vectors_storage[i][0] while accessing X_mg_vectors_storage[i][1] will give us the corresponding imaginary part

Referenced by oomph::HelmholtzMGPreconditioner< DIM >::direct_solve(), oomph::HelmholtzMGPreconditioner< DIM >::post_smooth(), oomph::HelmholtzMGPreconditioner< DIM >::pre_smooth(), and oomph::HelmholtzMGPreconditioner< DIM >::preconditioner_solve().

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