#include <geometric_multigrid.h>

Inheritance diagram for oomph::MGSolver< DIM >:

Public Types
typedef Smoother (	PreSmootherFactoryFctPt) ()

typedef Smoother (	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...

	MGSolver (MGProblem *mg_problem_pt)

	~MGSolver ()
	Delete any dynamically allocated data. More...

void	clean_up_memory ()
	Clean up anything that needs to be cleaned up. More...

void	set_self_test_vector ()

void	self_test ()

void	restriction_self_test ()

void	interpolation_self_test ()

void	plot (const unsigned &hierarchy_level, const DoubleVector &input_vector, const std::string &filename)

void	disable_v_cycle_output ()

void	disable_output ()

void	enable_v_cycle_output ()
	Enable the output of the V-cycle timings and other output. More...

void	enable_doc_everything ()
	Enable the output from anything that could have been suppressed. 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)

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	solve (Problem *const &problem_pt, DoubleVector &result)

unsigned	iterations () const
	Number of iterations. More...

unsigned &	max_iter ()
	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)

Public Member Functions inherited from oomph::IterativeLinearSolver
	IterativeLinearSolver ()

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

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

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

Preconditioner *&	preconditioner_pt ()
	Access function to preconditioner. More...

Preconditioner *const &	preconditioner_pt () const
	Access function to preconditioner (const version) More...

double &	tolerance ()
	Access to convergence tolerance. More...

unsigned &	max_iter ()
	Access to max. number of iterations. More...

void	enable_doc_convergence_history ()
	Enable documentation of the convergence history. More...

void	disable_doc_convergence_history ()
	Disable documentation of the convergence history. More...

void	open_convergence_history_file_stream (const std::string &file_name, const std::string &zone_title="")

void	close_convergence_history_file_stream ()
	Close convergence history output stream. More...

double	jacobian_setup_time () const

double	linear_solver_solution_time () const
	return the time taken to solve the linear system More...

virtual double	preconditioner_setup_time () const
	returns the the time taken to setup the preconditioner More...

void	enable_setup_preconditioner_before_solve ()
	Setup the preconditioner before the solve. More...

void	disable_setup_preconditioner_before_solve ()
	Don't set up the preconditioner before the solve. More...

void	enable_error_after_max_iter ()
	Throw an error if we don't converge within max_iter. More...

void	disable_error_after_max_iter ()
	Don't throw an error if we don't converge within max_iter (default). More...

void	enable_iterative_solver_as_preconditioner ()

void	disable_iterative_solver_as_preconditioner ()

Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...

	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...

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

virtual	~LinearSolver ()
	Empty virtual destructor. More...

void	enable_doc_time ()
	Enable documentation of solve times. More...

void	disable_doc_time ()
	Disable documentation of solve times. More...

bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...

bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...

virtual void	enable_resolve ()

virtual void	disable_resolve ()

virtual void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

virtual void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)

virtual void	solve_transpose (Problem *const &problem_pt, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)

virtual void	resolve (const DoubleVector &rhs, DoubleVector &result)

virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)

virtual void	enable_computation_of_gradient ()

void	disable_computation_of_gradient ()

void	reset_gradient ()

void	get_gradient (DoubleVector &gradient)
	function to access the gradient, provided it has been computed More...

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)

Protected Member Functions
void	mg_solve (DoubleVector &result)
	Linear solver. More...

void	modify_restriction_matrices ()
	Modify the restriction matrices. More...

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

Protected Attributes
unsigned	Nvcycle

MGProblem *	Mg_problem_pt

Vector< DoubleVector >	Rhs_mg_vectors_storage

bool	Suppress_v_cycle_output

bool	Suppress_all_output

std::ostream *	Stream_pt

Protected Attributes inherited from oomph::IterativeLinearSolver
bool	Doc_convergence_history

std::ofstream	Output_file_stream
	Output file stream for convergence history. More...

double	Tolerance
	Convergence tolerance. More...

unsigned	Max_iter
	Maximum number of iterations. More...

Preconditioner *	Preconditioner_pt
	Pointer to the preconditioner. More...

double	Jacobian_setup_time
	Jacobian setup time. More...

double	Solution_time
	linear solver solution time More...

double	Preconditioner_setup_time
	Preconditioner setup time. More...

bool	Setup_preconditioner_before_solve

bool	Throw_error_after_max_iter

bool	Use_iterative_solver_as_preconditioner
	Use the iterative solver as preconditioner. More...

bool	First_time_solve_when_used_as_preconditioner

Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve

bool	Doc_time
	Boolean flag that indicates whether the time taken. More...

bool	Compute_gradient

bool	Gradient_has_been_computed
	flag that indicates whether the gradient was computed or not More...

DoubleVector	Gradient_for_glob_conv_newton_solve

Private Member Functions
void	setup_mg_hierarchy ()
	Set up the MG hierarchy. More...

void	setup_mg_structures ()
	Set up the MG hierarchy structures. More...

void	setup_smoothers ()
	Set up the smoothers on all levels. More...

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...

unsigned	Nlevel
	The number of levels in the multigrid heirachy. More...

Vector< MGProblem * >	Mg_hierarchy
	Vector containing pointers to problems in hierarchy. More...

Vector< CRDoubleMatrix * >	Mg_matrices_storage_pt
	Vector to store the system matrices. More...

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< DoubleVector >	X_mg_vectors_storage
	Vector to store the solution vectors (X_mg) More...

Vector< DoubleVector >	Residual_mg_vectors_storage
	Vector to store the residual vectors. More...

Vector< DoubleVector >	Interpolation_self_test_vectors_storage

Vector< DoubleVector >	Restriction_self_test_vectors_storage

Vector< Smoother * >	Pre_smoothers_storage_pt
	Vector to store the pre-smoothers. More...

Vector< Smoother * >	Post_smoothers_storage_pt
	Vector to store the post-smoothers. More...

unsigned	Npre_smooth
	Number of pre-smoothing steps. More...

unsigned	Npost_smooth
	Number of post-smoothing steps. More...

bool	Doc_everything

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

bool	Has_been_solved

unsigned	V_cycle_counter
	Pointer to counter for V-cycles. More...

Additional Inherited Members
Static Protected Attributes inherited from oomph::IterativeLinearSolver
static IdentityPreconditioner	Default_preconditioner

Detailed Description

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

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

Member Typedef Documentation

◆ PostSmootherFactoryFctPt

template<unsigned DIM>

typedef Smoother*(* oomph::MGSolver< DIM >::PostSmootherFactoryFctPt) ()

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

◆ PreSmootherFactoryFctPt

template<unsigned DIM>

typedef Smoother*(* oomph::MGSolver< DIM >::PreSmootherFactoryFctPt) ()

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

Constructor & Destructor Documentation

◆ MGSolver()

template<unsigned DIM>

oomph::MGSolver< DIM >::MGSolver ( MGProblem * mg_problem_pt )

inline

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

       : Nvcycle(200),
         Mg_problem_pt(mg_problem_pt),
         Suppress_v_cycle_output(false),
         Suppress_all_output(false),
         Stream_pt(0),
         Pre_smoother_factory_function_pt(0),
         Post_smoother_factory_function_pt(0),
         Npre_smooth(2),
         Npost_smooth(2),
         Doc_everything(false),
         Has_been_setup(false),
         Has_been_solved(false)
     {
       // Set the tolerance in the base class
       this->Tolerance = 1.0e-09;
  
       // Set the pointer to the finest level as the first
       // entry in Mg_hierarchy
       Mg_hierarchy.push_back(Mg_problem_pt);
     } // End of MGSolver

References oomph::MGSolver< DIM >::Mg_hierarchy, oomph::MGSolver< DIM >::Mg_problem_pt, and oomph::IterativeLinearSolver::Tolerance.

◆ ~MGSolver()

template<unsigned DIM>

oomph::MGSolver< DIM >::~MGSolver ( )

inline

Delete any dynamically allocated data.

     {
       // Run the function written to clean up the memory
       clean_up_memory();
     } // End of ~MGSolver

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

Member Function Documentation

◆ clean_up_memory()

template<unsigned DIM>

void oomph::MGSolver< DIM >::clean_up_memory ( )

inlinevirtual

Clean up anything that needs to be cleaned up.

Reimplemented from oomph::LinearSolver.

Reimplemented in oomph::MGPreconditioner< DIM >.

     {
       // We only need to destroy data if the solver has been set up and
       // the data hasn't already been cleared
       if (Has_been_setup)
       {
         // Loop over all of the levels in the hierarchy
         for (unsigned i = 0; i < Nlevel - 1; i++)
         {
           // Delete the pre-smoother associated with this level
           delete Pre_smoothers_storage_pt[i];
  
           // Make it a null pointer
           Pre_smoothers_storage_pt[i] = 0;
  
           // Delete the post-smoother associated with this level
           delete Post_smoothers_storage_pt[i];
  
           // Make it a null pointer
           Post_smoothers_storage_pt[i] = 0;
  
           // Delete the system matrix associated with the i-th level
           delete Mg_matrices_storage_pt[i];
  
           // Make it a null pointer
           Mg_matrices_storage_pt[i] = 0;
         }
  
         // Loop over all but the coarsest of the levels in the hierarchy
         for (unsigned i = 0; i < Nlevel - 1; i++)
         {
           // Delete the interpolation matrix associated with the i-th level
           delete Interpolation_matrices_storage_pt[i];
  
           // Make it a null pointer
           Interpolation_matrices_storage_pt[i] = 0;
  
           // Delete the restriction matrix associated with the i-th level
           delete Restriction_matrices_storage_pt[i];
  
           // Make it a null pointer
           Restriction_matrices_storage_pt[i] = 0;
         }
  
         // If this solver has been set up then a hierarchy of problems
         // will have been set up. If the user chose to document everything
         // then the coarse-grid multigrid problems will have been kept alive
         // which means we now have to loop over the coarse-grid levels and
         // destroy them
         if (Doc_everything)
         {
           // Loop over the levels
           for (unsigned i = 1; i < Nlevel; i++)
           {
             // Delete the i-th level problem
             delete Mg_hierarchy[i];
  
             // Make the associated pointer a null pointer
             Mg_hierarchy[i] = 0;
           }
         } // if (Doc_everything)
  
         // Everything has been deleted now so we need to indicate that the
         // solver is not set up
         Has_been_setup = false;
       } // if (Has_been_setup)
     } // End of clean_up_memory

References oomph::MGSolver< DIM >::Doc_everything, oomph::MGSolver< DIM >::Has_been_setup, i, oomph::MGSolver< DIM >::Interpolation_matrices_storage_pt, oomph::MGSolver< DIM >::Mg_hierarchy, oomph::MGSolver< DIM >::Mg_matrices_storage_pt, oomph::MGSolver< DIM >::Nlevel, oomph::MGSolver< DIM >::Post_smoothers_storage_pt, oomph::MGSolver< DIM >::Pre_smoothers_storage_pt, and oomph::MGSolver< DIM >::Restriction_matrices_storage_pt.

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

◆ direct_solve()

template<unsigned DIM>

void oomph::MGSolver< DIM >::direct_solve ( )

inline

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

     {
       // Get solution by direct solve:
       Mg_matrices_storage_pt[Nlevel - 1]->solve(
         Rhs_mg_vectors_storage[Nlevel - 1], X_mg_vectors_storage[Nlevel - 1]);
     } // End of direct_solve

References oomph::MGSolver< DIM >::Mg_matrices_storage_pt, oomph::MGSolver< DIM >::Nlevel, oomph::MGSolver< DIM >::Rhs_mg_vectors_storage, and oomph::MGSolver< DIM >::X_mg_vectors_storage.

◆ disable_output()

template<unsigned DIM>

void oomph::MGSolver< 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

     {
       // Set the value of Doc_time (inherited from LinearSolver) to false
       Doc_time = false;
  
       // Enable the suppression of output from the V-cycle
       Suppress_v_cycle_output = true;
  
       // Enable the suppression of everything
       Suppress_all_output = true;
  
       // Store the output stream pointer
       Stream_pt = oomph_info.stream_pt();
  
       // Now set the oomph_info stream pointer to the null stream to
       // disable all possible output
       oomph_info.stream_pt() = &oomph_nullstream;
     } // End of disable_output

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

◆ disable_smoother_and_superlu_doc_time()

template<unsigned DIM>

void oomph::MGSolver< DIM >::disable_smoother_and_superlu_doc_time ( )

inline

Suppress the output of both smoothers and SuperLU.

     {
       // Loop over all levels of the hierarchy
       for (unsigned i = 0; i < Nlevel - 1; i++)
       {
         // Disable time documentation on each level (for each pre-smoother)
         Pre_smoothers_storage_pt[i]->disable_doc_time();
  
         // Disable time documentation on each level (for each post-smoother)
         Post_smoothers_storage_pt[i]->disable_doc_time();
       }
  
       // We only do a direct solve on the coarsest level so this is the only
       // place we need to silence SuperLU
       Mg_matrices_storage_pt[Nlevel - 1]
         ->linear_solver_pt()
         ->disable_doc_time();
     } // End of disable_smoother_and_superlu_doc_time

References i, oomph::MGSolver< DIM >::Mg_matrices_storage_pt, oomph::MGSolver< DIM >::Nlevel, oomph::MGSolver< DIM >::Post_smoothers_storage_pt, and oomph::MGSolver< DIM >::Pre_smoothers_storage_pt.

◆ disable_v_cycle_output()

template<unsigned DIM>

void oomph::MGSolver< DIM >::disable_v_cycle_output ( )

inline

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

     {
       // Set the value of Doc_time (inherited from LinearSolver) to false
       Doc_time = false;
  
       // Enable the suppression of output from the V-cycle
       Suppress_v_cycle_output = true;
     } // End of disable_v_cycle_output

References oomph::LinearSolver::Doc_time, and oomph::MGSolver< DIM >::Suppress_v_cycle_output.

◆ enable_doc_everything()

template<unsigned DIM>

void oomph::MGSolver< DIM >::enable_doc_everything ( )

inline

Enable the output from anything that could have been suppressed.

     {
       // Enable the documentation of everything (if this is set to TRUE then
       // the function self_test() will be run which outputs a solution
       // represented on each level of the hierarchy
       Doc_everything = true;
     } // End of enable_doc_everything

References oomph::MGSolver< DIM >::Doc_everything.

◆ enable_output()

template<unsigned DIM>

void oomph::MGSolver< DIM >::enable_output ( )

inline

Enable the output from anything that could have been suppressed.

     {
       // Enable time documentation
       Doc_time = true;
  
       // Enable output from everything during the full setup of the solver
       Suppress_all_output = false;
  
       // Enable output from the MG solver
       Suppress_v_cycle_output = false;
     } // End of enable_output

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

◆ enable_v_cycle_output()

template<unsigned DIM>

void oomph::MGSolver< DIM >::enable_v_cycle_output ( )

inline

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

     {
       // Enable time documentation
       Doc_time = true;
  
       // Enable output from the MG solver
       Suppress_v_cycle_output = false;
     } // End of enable_v_cycle_output

References oomph::LinearSolver::Doc_time, and oomph::MGSolver< DIM >::Suppress_v_cycle_output.

◆ full_setup()

template<unsigned DIM>

void oomph::MGSolver< DIM >::full_setup

Runs a full setup of the MG solver.

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

   {
     // Initialise the timer start variable
     double t_fs_start = 0.0;
  
     // If we're allowed to output
     if (!Suppress_all_output)
     {
       // Start the timer
       t_fs_start = TimingHelpers::timer();
  
       // Notify user that the hierarchy of levels is complete
       oomph_info
         << "\n===============Starting Multigrid Full Setup=============="
         << std::endl;
  
       // Notify user that the hierarchy of levels is complete
       oomph_info << "\nStarting the full setup of the multigrid solver."
                  << std::endl;
     }
  
 #ifdef PARANOID
     // PARANOID check - Make sure the dimension of the solver matches the
     // dimension of the elements used in the problems mesh
     if (dynamic_cast<FiniteElement*>(Mg_problem_pt->mesh_pt()->element_pt(0))
           ->dim() != DIM)
     {
       std::string err_strng = "The dimension of the elements used in the mesh ";
       err_strng += "does not match the dimension of the solver.";
       throw OomphLibError(
         err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
  
     // PARANOID check - The elements of the bulk mesh must all be refineable
     // elements otherwise we cannot deal with this
     if (Mg_problem_pt->mg_bulk_mesh_pt() != 0)
     {
       // Find the number of elements in the bulk mesh
       unsigned n_elements = Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
  
       // Loop over the elements in the mesh and ensure that they are
       // all refineable elements
       for (unsigned el_counter = 0; el_counter < n_elements; el_counter++)
       {
         // Upcast global mesh element to a refineable element
         RefineableElement* el_pt = dynamic_cast<RefineableElement*>(
           Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));
  
         // Check if the upcast worked or not; if el_pt is a null pointer the
         // element is not refineable
         if (el_pt == 0)
         {
           throw OomphLibError(
             "Element in global mesh could not be upcast to a refineable "
             "element. We cannot deal with elements that are not refineable.",
             OOMPH_CURRENT_FUNCTION,
             OOMPH_EXCEPTION_LOCATION);
         }
       }
     }
     else
     {
       throw OomphLibError(
         "The provided bulk mesh pointer is set to be a null pointer. "
         "The multigrid solver operates on the bulk mesh thus a pointer "
         "to the correct mesh must be given.",
         OOMPH_CURRENT_FUNCTION,
         OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // If this is not the first Newton step then we will already have things
     // in storage. If this is the case, delete them
     clean_up_memory();
  
     // Resize the Mg_hierarchy vector
     Mg_hierarchy.resize(1, 0);
  
     // Set the pointer to the finest level as the first entry in Mg_hierarchy
     Mg_hierarchy[0] = Mg_problem_pt;
  
     // Create the hierarchy of levels
     setup_mg_hierarchy();
  
     // Set up the interpolation and restriction matrices
     setup_transfer_matrices();
  
     // Set up the data structures on each level, i.e. the system matrix,
     // LHS and RHS vectors
     setup_mg_structures();
  
     // Set up the smoothers on all of the levels
     setup_smoothers();
  
     // If we do not want to document everything we want to delete all the
     // coarse-grid problems
     if (!Doc_everything)
     {
       // Loop over all of the coarser levels
       for (unsigned i = 1; i < Nlevel; i++)
       {
         // Delete the i-th coarse-grid MGProblem
         delete Mg_hierarchy[i];
  
         // Set it to be a null pointer
         Mg_hierarchy[i] = 0;
       }
     }
     // Otherwise, document everything!
     else
     {
       // If the user wishes to document everything we run the self-test
       self_test();
     } // if (!Doc_everything)
  
     // Indicate that the full setup has been completed
     Has_been_setup = true;
  
     // If we're allowed to output to the screen
     if (!Suppress_all_output)
     {
       // Output the time taken to complete the full setup
       double t_fs_end = TimingHelpers::timer();
       double full_setup_time = t_fs_end - t_fs_start;
  
       // Output the CPU time
       oomph_info << "\nCPU time for full setup [sec]: " << full_setup_time
                  << std::endl;
  
       // Notify user that the hierarchy of levels is complete
       oomph_info
         << "\n===============Multigrid Full Setup Complete=============="
         << std::endl;
     } // if (!Suppress_all_output)
   } // End of full_setup

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

Referenced by oomph::MGPreconditioner< DIM >::setup(), and oomph::MGSolver< DIM >::solve().

◆ interpolate_and_correct()

template<unsigned DIM>

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

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

   {
 #ifdef PARANOID
     // Check to make sure we can actually restrict the vector
     if (!(level > 0))
     {
       throw OomphLibError("Input level exceeds the possible parameter choice.",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Build distribution of a temporary vector
     DoubleVector temp_soln(X_mg_vectors_storage[level - 1].distribution_pt());
  
     // Interpolate the solution vector
     Interpolation_matrices_storage_pt[level - 1]->multiply(
       X_mg_vectors_storage[level], temp_soln);
  
     // Update
     X_mg_vectors_storage[level - 1] += temp_soln;
   } // End of interpolate_and_correct

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ interpolation_matrix_set() [1/2]

template<unsigned DIM>

void oomph::MGSolver< 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.

     {
       // Dynamically allocate the interpolation matrix
       Interpolation_matrices_storage_pt[level] = new CRDoubleMatrix;
  
       // Build the matrix
       Interpolation_matrices_storage_pt[level]->build_without_copy(
         ncol, nnz, value, col_index, row_st);
     } // End of interpolation_matrix_set

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

◆ interpolation_matrix_set() [2/2]

template<unsigned DIM>

void oomph::MGSolver< 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.

     {
       // Dynamically allocate the interpolation matrix
       Interpolation_matrices_storage_pt[level] = new CRDoubleMatrix;
  
       // Make the distribution pointer
       LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
         Mg_hierarchy[level]->communicator_pt(), nrow, false);
  
 #ifdef PARANOID
 #ifdef OOMPH_HAS_MPI
       // Set up the warning messages
       std::string warning_message =
         "Setup of interpolation matrix distribution ";
       warning_message += "has not been tested with MPI.";
  
       // If we're not running the code in serial
       if (dist_pt->communicator_pt()->nproc() > 1)
       {
         // Throw a warning
         OomphLibWarning(
           warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
 #endif
 #endif
  
       // Build the matrix itself
       Interpolation_matrices_storage_pt[level]->build(
         dist_pt, ncol, value, col_index, row_st);
  
       // Delete the newly created distribution pointer
       delete dist_pt;
  
       // Make it a null pointer
       dist_pt = 0;
     } // End of interpolation_matrix_set

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

◆ interpolation_self_test()

template<unsigned DIM>

void oomph::MGSolver< DIM >::interpolation_self_test

Make a self-test to make sure that the interpolation matrices are doing the same thing to interpolate the vectors up.

Function which implements a self-test to interpolate a vector up all of levels of the MG hierarchy and outputs the restricted vectors to file

   {
     // Set the start of the output file name. The full names will
     // set after we calculate the interpolated vectors
     std::string outputfile = "RESLT/interpolation_self_test";
  
     // Loop over the levels of the hierarchy in reverse order
     for (unsigned level = Nlevel - 1; level > 0; level--)
     {
       // Interpolate the vector up a level
       Interpolation_matrices_storage_pt[level - 1]->multiply(
         Interpolation_self_test_vectors_storage[level],
         Interpolation_self_test_vectors_storage[level - 1]);
     } // End of the for loop over the hierarchy levels
  
     for (unsigned level = 0; level < Nlevel; level++)
     {
       // Set output file name
       std::string filename = outputfile;
       std::stringstream string;
       string << level;
       filename += string.str();
       filename += ".dat";
  
       // Plot the RHS vector on the current level
       plot(level, Interpolation_self_test_vectors_storage[level], filename);
  
     } // End of for-loop to output the RHS vector on each level
   } // End of interpolation_self_test

References MergeRestartFiles::filename, PlanarWave::plot(), and oomph::Global_string_for_annotation::string().

◆ iterations()

template<unsigned DIM>

unsigned oomph::MGSolver< DIM >::iterations ( ) const

inlinevirtual

Number of iterations.

Implements oomph::IterativeLinearSolver.

     {
       // Return the number of V-cycles which have been done
       return V_cycle_counter;
     } // End of iterations

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

◆ level_up_local_coord_of_node() [1/3]

void oomph::MGSolver< 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.

   {
     // If the element is unrefined between the levels the local coordinate
     // of the node in one element is the same as that in the other element
     // therefore we only need to perform calculations if the levels are
     // different (i.e. son_type is not OMEGA)
     if (son_type != Tree::OMEGA)
     {
       // Scale the local coordinate from the range [-1,1]x[-1,1] to the range
       // [0,1]x[0,1] to match the width of the local coordinate range of the
       // fine element from the perspective of the father element. This
       // then simply requires a shift of the coordinates to match which type
       // of son element we're dealing with
       s[0] = (s[0] + 1.0) / 2.0;
       s[1] = (s[1] + 1.0) / 2.0;
  
       // Cases: The son_type determines how the local coordinates should be
       // shifted to give the local coordinates in the coarse mesh element
       switch (son_type)
       {
           // If we're dealing with the bottom-left element we need to shift
           // the range [0,1]x[0,1] to [-1,0]x[-1,0]
         case QuadTreeNames::SW:
           s[0] -= 1;
           s[1] -= 1;
           break;
  
           // If we're dealing with the bottom-right element we need to shift
           // the range [0,1]x[0,1] to [0,1]x[-1,0]
         case QuadTreeNames::SE:
           s[1] -= 1;
           break;
  
           // If we're dealing with the top-right element we need to shift the
           // range [0,1]x[0,1] to [0,1]x[0,1], i.e. nothing needs to be done
         case QuadTreeNames::NE:
           break;
  
           // If we're dealing with the top-left element we need to shift
           // the range [0,1]x[0,1] to [-1,0]x[0,1]
         case QuadTreeNames::NW:
           s[0] -= 1;
           break;
       }
     } // if son_type!=Tree::OMEGA
   } // End of level_up_local_coord_of_node

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

◆ level_up_local_coord_of_node() [2/3]

void oomph::MGSolver< 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.

   {
     // If the element is unrefined between the levels the local coordinate
     // of the node in one element is the same as that in the other element
     // therefore we only need to perform calculations if the levels are
     // different (i.e. son_type is not OMEGA)
     if (son_type != Tree::OMEGA)
     {
       // Scale the local coordinate from the range [-1,1]x[-1,1]x[-1,1]
       // to the range [0,1]x[0,1]x[0,1] to match the width of the local
       // coordinate range of the fine element from the perspective of
       // the father element. This then simply requires a shift of the
       // coordinates to match which type of son element we're dealing with
       s[0] = (s[0] + 1.0) / 2.0;
       s[1] = (s[1] + 1.0) / 2.0;
       s[2] = (s[2] + 1.0) / 2.0;
  
       // Cases: The son_type determines how the local coordinates should be
       // shifted to give the local coordinates in the coarse mesh element
       switch (son_type)
       {
         case OcTreeNames::LDF:
           s[0] -= 1;
           s[1] -= 1;
           break;
  
         case OcTreeNames::LDB:
           s[0] -= 1;
           s[1] -= 1;
           s[2] -= 1;
           break;
  
         case OcTreeNames::LUF:
           s[0] -= 1;
           break;
  
         case OcTreeNames::LUB:
           s[0] -= 1;
           s[2] -= 1;
           break;
  
         case OcTreeNames::RDF:
           s[1] -= 1;
           break;
  
         case OcTreeNames::RDB:
           s[1] -= 1;
           s[2] -= 1;
           break;
  
         case OcTreeNames::RUF:
           break;
  
         case OcTreeNames::RUB:
           s[2] -= 1;
           break;
       }
     } // if son_type!=Tree::OMEGA
   } // End of level_up_local_coord_of_node

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]

template<unsigned DIM>

void oomph::MGSolver< 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

◆ max_iter()

template<unsigned DIM>

unsigned& oomph::MGSolver< DIM >::max_iter ( )

inline

Number of iterations.

     {
       // Return the number of V-cycles which have been done
       return Nvcycle;
     } // End of iterations

References oomph::MGSolver< DIM >::Nvcycle.

◆ mg_solve()

template<unsigned DIM>

void oomph::MGSolver< DIM >::mg_solve ( DoubleVector & result )

protected

Linear solver.

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 MGPreconditioner derived class to use the solver

   {
     // If we're allowed to time things
     double t_mg_start = 0.0;
     if (!Suppress_v_cycle_output)
     {
       // Start the clock!
       t_mg_start = TimingHelpers::timer();
     }
  
     // Current level
     unsigned finest_level = 0;
  
     // Initialise the V-cycle counter
     V_cycle_counter = 0;
  
     // Calculate the norm of the residual then output
     double normalised_residual_norm = residual_norm(finest_level);
     if (!Suppress_v_cycle_output)
     {
       oomph_info << "\nResidual on finest level for V-cycle: "
                  << normalised_residual_norm << std::endl;
     }
  
     // Outer loop over V-cycles
     //-------------------------
     // While the tolerance is not met and the maximum number of
     // V-cycles has not been completed
     while ((normalised_residual_norm > (this->Tolerance)) &&
            (V_cycle_counter != Nvcycle))
     {
       if (!Suppress_v_cycle_output)
       {
         // Output the V-cycle counter
         oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
       }
  
       // Loop downwards over all levels that have coarser levels beneath them
       //---------------------------------------------------------------------
       for (unsigned i = 0; i < Nlevel - 1; i++)
       {
         // Initialise x_mg and residual_mg to 0.0 except for the finest level
         // since x_mg contains the current approximation to the solution and
         // residual_mg contains the RHS vector on the finest level
         if (i != 0)
         {
           // Loop over the entries in x_mg and residual_mg
           X_mg_vectors_storage[i].initialise(0.0);
           Residual_mg_vectors_storage[i].initialise(0.0);
         }
  
         // Perform a few pre-smoothing steps and return vector that contains
         // the residuals of the linear system at this level.
         pre_smooth(i);
  
         // Restrict the residual to the next coarser mesh and
         // assign it to the RHS vector at that level
         restrict_residual(i);
       } // Moving down the V-cycle
  
       // Reached the lowest level: Do a direct solve, using the RHS vector
       //------------------------------------------------------------------
       // obtained by restriction from above.
       //------------------------------------
       direct_solve();
  
       // Loop upwards over all levels that have finer levels above them
       //---------------------------------------------------------------
       for (unsigned i = Nlevel - 1; i > 0; i--)
       {
         // Interpolate solution at current level onto
         // next finer mesh and correct the solution x at that level
         interpolate_and_correct(i);
  
         // Perform a few post-smoothing steps (ignore
         // vector that contains the residuals of the linear system
         // at this level)
         post_smooth(i - 1);
       }
  
       // Calculate the new residual norm then output (if allowed)
       normalised_residual_norm = residual_norm(finest_level);
  
       // If required, we will document the convergence history to screen or file
       // (if stream open)
       if (Doc_convergence_history)
       {
         if (!Output_file_stream.is_open())
         {
           oomph_info << V_cycle_counter << " " << normalised_residual_norm
                      << std::endl;
         }
         else
         {
           Output_file_stream << V_cycle_counter << " "
                              << normalised_residual_norm << std::endl;
         }
       } // if (Doc_convergence_history)
  
       // Set counter for number of cycles (for doc)
       V_cycle_counter++;
  
       // Print the residual on the finest level
       if (!Suppress_v_cycle_output)
       {
         oomph_info << "Residual on finest level of V-cycle: "
                    << normalised_residual_norm << std::endl;
       }
     } // End of the V-cycles
  
     // Copy the solution into the result vector
     result = X_mg_vectors_storage[finest_level];
  
     // Need an extra line space if V-cycle output is suppressed
     if (!Suppress_v_cycle_output)
     {
       oomph_info << std::endl;
     }
  
     // If all output is to be suppressed
     if (!Suppress_all_output)
     {
       // Output number of V-cycles taken to solve
       if (normalised_residual_norm < (this->Tolerance))
       {
         // Indicate that the problem has been solved
         Has_been_solved = true;
       }
     } // if (!Suppress_all_output)
  
     // If the V-cycle output isn't suppressed
     if (!Suppress_v_cycle_output)
     {
       // Stop the clock
       double t_mg_end = TimingHelpers::timer();
       double total_G_setup_time = double(t_mg_end - t_mg_start);
       oomph_info << "CPU time for MG solver [sec]: " << total_G_setup_time
                  << std::endl;
     }
   } // End of mg_solve

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

Referenced by oomph::MGPreconditioner< DIM >::preconditioner_solve(), and oomph::MGSolver< DIM >::solve().

◆ modify_restriction_matrices()

template<unsigned DIM>

void oomph::MGSolver< DIM >::modify_restriction_matrices

protected

Modify the restriction matrices.

Normalise the rows of the restriction matrices to avoid amplifications when projecting to the coarser level

   {
     // Create a null pointer
     CRDoubleMatrix* restriction_matrix_pt = 0;
  
     // Loop over the levels
     for (unsigned level = 0; level < Nlevel - 1; level++)
     {
       // Store a pointer to the (level)-th restriction matrix
       restriction_matrix_pt = Restriction_matrices_storage_pt[level];
  
       // Get access to the row start data
       const int* row_start_pt = restriction_matrix_pt->row_start();
  
       // Get access to the matrix entries
       double* value_pt = restriction_matrix_pt->value();
  
       // Initialise an auxiliary variable to store the index of the start
       // of the i-th row
       unsigned start_index = 0;
  
       // Initialise an auxiliary variable to store the index of the start
       // of the (i+1)-th row
       unsigned end_index = 0;
  
       // Store the number of rows in the matrix
       unsigned n_row = restriction_matrix_pt->nrow();
  
       // Loop over the rows of the matrix
       for (unsigned i = 0; i < n_row; i++)
       {
         // Index associated with the start of the i-th row
         start_index = row_start_pt[i];
  
         // Index associated with the start of the (i+1)-th row
         end_index = row_start_pt[i + 1];
  
         // Variable to store the sum of the absolute values of the off-diagonal
         // entries in the i-th row
         double row_sum = 0.0;
  
         // Loop over the entries in the row
         for (unsigned j = start_index; j < end_index; j++)
         {
           // Add the absolute value of this entry to the off-diagonal sum
           row_sum += value_pt[j];
         }
  
         // Loop over the entries in the row
         for (unsigned j = start_index; j < end_index; j++)
         {
           // Normalise the row entry
           value_pt[j] /= row_sum;
         }
       } // for (unsigned i=0;i<n_row;i++)
     } // for (unsigned level=0;level<Nlevel-1;level++)
   } // End of modify_restriction_matrices

References i, j, oomph::CRDoubleMatrix::nrow(), oomph::CRDoubleMatrix::row_start(), and oomph::CRDoubleMatrix::value().

◆ npost_smooth()

template<unsigned DIM>

unsigned& oomph::MGSolver< DIM >::npost_smooth ( )

inline

Return the number of post-smoothing iterations (lvalue)

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

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

◆ npre_smooth()

template<unsigned DIM>

unsigned& oomph::MGSolver< DIM >::npre_smooth ( )

inline

Return the number of pre-smoothing iterations (lvalue)

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

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

◆ plot()

template<unsigned DIM>

void oomph::MGSolver< DIM >::plot	(	const unsigned &	hierarchy_level,
		const DoubleVector &	input_vector,
		const std::string &	filename
	)

Given a level in the hierarchy, an input vector and a filename this function will document the given vector according to the structure of the mesh on the given level

Plots the input vector (assuming we're dealing with scalar nodal data, otherwise I don't know how to implement this...)

   {
     // Setup an output file
     std::ofstream some_file;
     some_file.open(filename.c_str());
  
     // Set up a pointer to the refineable mesh
     TreeBasedRefineableMeshBase* bulk_mesh_pt =
       Mg_hierarchy[hierarchy_level]->mg_bulk_mesh_pt();
  
     // Find the number of elements in the bulk mesh
     unsigned n_el = bulk_mesh_pt->nelement();
  
     // Get the dimension of the problem (assumed to be the dimension of any
     // node in the mesh)
     unsigned n_dim = bulk_mesh_pt->node_pt(0)->ndim();
  
     // Find the number of nodes in an element
     unsigned nnod = bulk_mesh_pt->finite_element_pt(0)->nnode();
  
     // Find the number of nodes in an element in one direction
     unsigned nnod_1d = bulk_mesh_pt->finite_element_pt(0)->nnode_1d();
  
     // Loop over all elements
     for (unsigned e = 0; e < n_el; e++)
     {
       // Get pointer to element
       FiniteElement* el_pt = bulk_mesh_pt->finite_element_pt(e);
  
       // Input string for tecplot zone header (when plotting nnode_1d
       // points in each coordinate direction)
       some_file << el_pt->tecplot_zone_string(nnod_1d) << std::endl;
  
       // Loop over nodes
       for (unsigned j = 0; j < nnod; j++)
       {
         // Pointer to node
         Node* nod_pt = el_pt->node_pt(j);
  
         // Sanity check
         if (nod_pt->nvalue() != 1)
         {
           std::ostringstream error_stream;
  
           error_stream << "Sorry, not sure what to do here!" << nod_pt->nvalue()
                        << std::endl;
  
           // Output the value of dimension to screen
           error_stream << "The dimension is: " << n_dim << std::endl;
  
           // Output the values of x_i for all i
           for (unsigned i = 0; i < n_dim; i++)
           {
             error_stream << nod_pt->x(i) << " ";
           }
  
           // End line
           error_stream << std::endl;
  
           // Throw an error to indicate that it was
           // not possible to plot the data
           throw OomphLibError(error_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
         }
  
         // Output the coordinates of the node
         for (unsigned i = 0; i < n_dim; i++)
         {
           some_file << nod_pt->x(i) << " ";
         }
  
         // Free or pinned?
         int e = nod_pt->eqn_number(0);
         if (e >= 0)
         {
           // Output the e-th entry if the node is free
           some_file << input_vector[e] << std::endl;
         }
         else
         {
           // Boundary node
           if (nod_pt->is_on_boundary())
           {
             // On the finest level the boundary data is pinned so set to 0
             if (hierarchy_level == 0)
             {
               some_file << 0.0 << std::endl;
             }
             // On any other level we output this data since it corresponds
             // to the correction on the boundary nodes
             else
             {
               some_file << input_vector[e] << std::endl;
             }
           }
           // Hanging node
           else if (nod_pt->is_hanging())
           {
             // Allocate space for a double to store the weighted sum
             // of the surrounding master nodes of the hanging node
             double hang_sum = 0.0;
  
             // Find the number of master nodes of the hanging
             // the node in the reference element
             HangInfo* hang_info_pt = nod_pt->hanging_pt();
             unsigned nmaster = hang_info_pt->nmaster();
  
             // Loop over the master nodes
             for (unsigned i_master = 0; i_master < nmaster; i_master++)
             {
               // Set up a pointer to the master node
               Node* master_node_pt = hang_info_pt->master_node_pt(i_master);
  
               // Get the global equation number of this node
               int master_jj = master_node_pt->eqn_number(0);
  
               // Is the master node a proper d.o.f.?
               if (master_jj >= 0)
               {
                 // If the weight of the master node is non-zero
                 if (hang_info_pt->master_weight(i_master) != 0.0)
                 {
                   hang_sum += hang_info_pt->master_weight(i_master) *
                               input_vector[master_jj];
                 }
               } // if (master_jj>=0)
             } // for (unsigned i_master=0;i_master<nmaster;i_master++)
  
             // Output the weighted sum of the surrounding master nodes
             some_file << hang_sum << std::endl;
           }
         } // if (e>=0)
       } // for (unsigned j=0;j<nnod;j++)
     } // for (unsigned e=0;e<n_el;e++)
  
     // Close output file
     some_file.close();
   } // End of plot

References e(), oomph::Data::eqn_number(), MergeRestartFiles::filename, oomph::Mesh::finite_element_pt(), oomph::Node::hanging_pt(), i, oomph::Node::is_hanging(), oomph::Node::is_on_boundary(), j, oomph::HangInfo::master_node_pt(), oomph::HangInfo::master_weight(), oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::HangInfo::nmaster(), oomph::FiniteElement::nnode(), oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::FiniteElement::tecplot_zone_string(), and oomph::Node::x().

◆ post_smooth()

template<unsigned DIM>

void oomph::MGSolver< 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 MGProblem constructor) which can be overloaded for specific problem.

     {
       // Run post-smoother 'max_iter' times
       Post_smoothers_storage_pt[level]->smoother_solve(
         Rhs_mg_vectors_storage[level], X_mg_vectors_storage[level]);
     } // End of post_smooth

References oomph::MGSolver< DIM >::Post_smoothers_storage_pt, oomph::MGSolver< DIM >::Rhs_mg_vectors_storage, and oomph::MGSolver< DIM >::X_mg_vectors_storage.

◆ pre_smooth()

template<unsigned DIM>

void oomph::MGSolver< 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 MGProblem constructor) which can be overloaded for a specific problem.

     {
       // Run pre-smoother 'max_iter' times
       Pre_smoothers_storage_pt[level]->smoother_solve(
         Rhs_mg_vectors_storage[level], X_mg_vectors_storage[level]);
  
       // Calculate the residual r=b-Ax and assign it
       Mg_matrices_storage_pt[level]->residual(
         X_mg_vectors_storage[level],
         Rhs_mg_vectors_storage[level],
         Residual_mg_vectors_storage[level]);
     } // End of pre_smooth

References oomph::MGSolver< DIM >::Mg_matrices_storage_pt, oomph::MGSolver< DIM >::Pre_smoothers_storage_pt, oomph::MGSolver< DIM >::Residual_mg_vectors_storage, oomph::MGSolver< DIM >::Rhs_mg_vectors_storage, and oomph::MGSolver< DIM >::X_mg_vectors_storage.

◆ residual_norm()

template<unsigned DIM>

double oomph::MGSolver< DIM >::residual_norm ( const unsigned & level )

inline

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

     {
       // And zero the entries of residual
       Residual_mg_vectors_storage[level].initialise(0.0);
  
       // Get the residual
       Mg_matrices_storage_pt[level]->residual(
         X_mg_vectors_storage[level],
         Rhs_mg_vectors_storage[level],
         Residual_mg_vectors_storage[level]);
  
       // Return the norm of the residual
       return Residual_mg_vectors_storage[level].norm();
     } // End of residual_norm

References oomph::MGSolver< DIM >::Mg_matrices_storage_pt, oomph::MGSolver< DIM >::Residual_mg_vectors_storage, oomph::MGSolver< DIM >::Rhs_mg_vectors_storage, and oomph::MGSolver< DIM >::X_mg_vectors_storage.

◆ restrict_residual()

template<unsigned DIM>

void oomph::MGSolver< 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)

   {
 #ifdef PARANOID
     // Check to make sure we can actually restrict the vector
     if (!(level < Nlevel - 1))
     {
       throw OomphLibError("Input exceeds the possible parameter choice.",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // Multiply the residual vector by the restriction matrix on the level-th
     // level (to restrict the vector down to the next coarser level)
     Restriction_matrices_storage_pt[level]->multiply(
       Residual_mg_vectors_storage[level], Rhs_mg_vectors_storage[level + 1]);
   } // End of restrict_residual

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ restriction_self_test()

template<unsigned DIM>

void oomph::MGSolver< DIM >::restriction_self_test

Make a self-test to make sure that the interpolation matrices are doing the same thing to restrict the vectors down through the heirachy.

Function which implements a self-test to restrict the residual vector down all of the coarser grids and output the restricted vectors to file

   {
     // Set the start of the output file name. The full names will
     // set after we calculate the restricted vectors
     std::string outputfile = "RESLT/restriction_self_test";
  
     // Loop over the levels of the hierarchy
     for (unsigned level = 0; level < Nlevel - 1; level++)
     {
       // Restrict the vector down to the next level
       Restriction_matrices_storage_pt[level]->multiply(
         Restriction_self_test_vectors_storage[level],
         Restriction_self_test_vectors_storage[level + 1]);
     } // End of the for loop over the hierarchy levels
  
     // Loop over the levels of hierarchy to plot the restricted vectors
     for (unsigned level = 0; level < Nlevel; level++)
     {
       // Set output file name
       std::string filename = outputfile;
       std::stringstream string;
       string << level;
       filename += string.str();
       filename += ".dat";
  
       // Plot the RHS vector on the current level
       plot(level, Restriction_self_test_vectors_storage[level], filename);
  
     } // End of for-loop to output the RHS vector on each level
   } // End of restriction_self_test

References MergeRestartFiles::filename, PlanarWave::plot(), and oomph::Global_string_for_annotation::string().

◆ self_test()

template<unsigned DIM>

void oomph::MGSolver< DIM >::self_test

Makes a vector, restricts it down the levels of the hierarchy and documents it at each level. After this is done the vector is interpolated up the levels of the hierarchy with the output being documented at each level

   {
     // Start the timer
     double t_st_start = TimingHelpers::timer();
  
     // Notify user that the hierarchy of levels is complete
     oomph_info << "\nStarting the multigrid solver self-test." << std::endl;
  
     // Resize the vector storing all of the restricted vectors
     Restriction_self_test_vectors_storage.resize(Nlevel);
  
     // Resize the vector storing all of the interpolated vectors
     Interpolation_self_test_vectors_storage.resize(Nlevel);
  
     // Find the number of dofs on the finest level
     unsigned n_dof = X_mg_vectors_storage[0].nrow();
  
     // Need to set up the distribution of the finest level vector
     LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
       Mg_problem_pt->communicator_pt(), n_dof, false);
  
 #ifdef PARANOID
 #ifdef OOMPH_HAS_MPI
     // Set up the warning messages
     std::string warning_message = "Setup of distribution has not been ";
     warning_message += "tested with MPI.";
  
     // If we're not running the code in serial
     if (dist_pt->communicator_pt()->nproc() > 1)
     {
       // Throw a warning
       OomphLibWarning(
         warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
 #endif
  
     // Build the vector using the distribution of the restriction matrix
     Restriction_self_test_vectors_storage[0].build(dist_pt);
  
     // Delete the distribution pointer
     delete dist_pt;
  
     // Make it a null pointer
     dist_pt = 0;
  
     // Now assign the values to the first vector. This will be restricted down
     // the levels of the hierarchy then back up
     set_self_test_vector();
  
     // Call the restriction self-test
     restriction_self_test();
  
     // Set the coarest level vector in Restriction_self_test_vectors_storage
     // to be the last entry in Interpolation_self_test_vectors_storage. This
     // will be interpolated up to the finest level in interpolation_self_test()
     Interpolation_self_test_vectors_storage[Nlevel - 1] =
       Restriction_self_test_vectors_storage[Nlevel - 1];
  
     // Call the interpolation self-test
     interpolation_self_test();
  
     // Destroy all of the created restriction self-test vectors
     Restriction_self_test_vectors_storage.resize(0);
  
     // Destroy all of the created interpolation self-test vectors
     Interpolation_self_test_vectors_storage.resize(0);
  
     // Stop the clock
     double t_st_end = TimingHelpers::timer();
     double total_st_time = double(t_st_end - t_st_start);
     oomph_info << "\nCPU time for self-test [sec]: " << total_st_time
                << std::endl;
  
     // Notify user that the hierarchy of levels is complete
     oomph_info << "\n====================Self-Test Complete===================="
                << std::endl;
   } // End of self_test

References oomph::LinearAlgebraDistribution::communicator_pt(), oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Global_string_for_annotation::string(), and oomph::TimingHelpers::timer().

◆ set_post_smoother_factory_function()

template<unsigned DIM>

void oomph::MGSolver< DIM >::set_post_smoother_factory_function ( PostSmootherFactoryFctPt post_smoother_fn )

inline

Access function to set the post-smoother creation function.

     {
       // Assign the function pointer
       Post_smoother_factory_function_pt = post_smoother_fn;
     }

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

Referenced by FluxPoissonMGProblem< ELEMENT, MESH >::set_multigrid_solver(), and UnitCubePoissonMGProblem< ELEMENT, MESH >::set_multigrid_solver().

◆ set_pre_smoother_factory_function()

template<unsigned DIM>

void oomph::MGSolver< DIM >::set_pre_smoother_factory_function ( PreSmootherFactoryFctPt pre_smoother_fn )

inline

Access function to set the pre-smoother creation function.

     {
       // Assign the function pointer
       Pre_smoother_factory_function_pt = pre_smoother_fn;
     }

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

Referenced by FluxPoissonMGProblem< ELEMENT, MESH >::set_multigrid_solver(), and UnitCubePoissonMGProblem< ELEMENT, MESH >::set_multigrid_solver().

◆ set_restriction_matrices_as_interpolation_transposes()

template<unsigned DIM>

void oomph::MGSolver< 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

     {
       for (unsigned i = 0; i < Nlevel - 1; i++)
       {
         // Dynamically allocate the restriction matrix
         Restriction_matrices_storage_pt[i] = new CRDoubleMatrix;
  
         // Set the restriction matrix to be the transpose of the
         // interpolation matrix
         Interpolation_matrices_storage_pt[i]->get_matrix_transpose(
           Restriction_matrices_storage_pt[i]);
       }
     } // End of set_restriction_matrices_as_interpolation_transposes

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

◆ set_self_test_vector()

template<unsigned DIM>

void oomph::MGSolver< DIM >::set_self_test_vector

Makes a vector which will be used in the self-test. Is currently set to make the entries of the vector represent a plane wave propagating at an angle of 45 degrees

Sets the initial vector to be used in the restriction and interpolation self-tests

   {
     // Set up a pointer to the refineable mesh
     TreeBasedRefineableMeshBase* bulk_mesh_pt =
       Mg_problem_pt->mg_bulk_mesh_pt();
  
     // Find the number of elements in the bulk mesh
     unsigned n_el = bulk_mesh_pt->nelement();
  
     // Get the dimension of the problem (assumed to be the dimension of any
     // node in the mesh)
     unsigned n_dim = bulk_mesh_pt->node_pt(0)->ndim();
  
     // Find the number of nodes in an element
     unsigned nnod = bulk_mesh_pt->finite_element_pt(0)->nnode();
  
     // Loop over all elements
     for (unsigned e = 0; e < n_el; e++)
     {
       // Get pointer to element
       FiniteElement* el_pt = bulk_mesh_pt->finite_element_pt(e);
  
       // Loop over nodes
       for (unsigned j = 0; j < nnod; j++)
       {
         // Pointer to node
         Node* nod_pt = el_pt->node_pt(j);
  
         // Sanity check
         if (nod_pt->nvalue() != 1)
         {
           // Set up an output stream
           std::ostringstream error_stream;
  
           // Output the error message
           error_stream << "Sorry, not sure what to do here! I can't deal with "
                        << nod_pt->nvalue() << " values!" << std::endl;
  
           // Throw an error to indicate that it was not possible to plot the
           // data
           throw OomphLibError(error_stream.str(),
                               OOMPH_CURRENT_FUNCTION,
                               OOMPH_EXCEPTION_LOCATION);
         }
  
         // Free or pinned?
         int eqn_num = nod_pt->eqn_number(0);
  
         // If it's a free node
         if (eqn_num >= 0)
         {
           // Initialise the variable coordinate_sum
           double coordinate_sum = 0.0;
  
           // Loop over the coordinates of the node
           for (unsigned i = 0; i < n_dim; i++)
           {
             // Increment coordinate_sum by the value of x(i)
             coordinate_sum += nod_pt->x(i);
           }
  
           // Store the value of pi in cache
           double pi = MathematicalConstants::Pi;
  
           // Make the vector represent a constant function
           Restriction_self_test_vectors_storage[0][eqn_num] =
             sin(pi * (nod_pt->x(0))) * sin(pi * (nod_pt->x(1)));
         }
       } // for (unsigned j=0;j<nnod;j++)
     } // for (unsigned e=0;e<n_el;e++)
  
     // // Get the size of the vector
     // unsigned n_row=Restriction_self_test_vectors_storage[0].nrow();
  
     // // Loop over the entries of the vector
     // for (unsigned i=0;i<n_row;i++)
     // {
     //  // Make the vector represent a constant function
     //  Restriction_self_test_vectors_storage[0][i]=1.0;
     // }
   } // End of set_self_test_vector

References e(), oomph::Data::eqn_number(), oomph::Mesh::finite_element_pt(), i, j, oomph::Node::ndim(), oomph::Mesh::nelement(), oomph::FiniteElement::nnode(), oomph::FiniteElement::node_pt(), oomph::Mesh::node_pt(), oomph::Data::nvalue(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, constants::pi, oomph::MathematicalConstants::Pi, sin(), and oomph::Node::x().

◆ setup_interpolation_matrices()

template<unsigned DIM>

void oomph::MGSolver< DIM >::setup_interpolation_matrices

Setup the interpolation matrix on each level.

Setup the interpolation matrices.

   {
     // Variable to hold the number of sons in an element
     unsigned n_sons;
  
     // Number of son elements
     if (DIM == 2)
     {
       n_sons = 4;
     }
     else if (DIM == 3)
     {
       n_sons = 8;
     }
     else
     {
       std::ostringstream error_message_stream;
       error_message_stream << "DIM should be 2 or 3 not " << DIM << std::endl;
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // Vector of local coordinates in the element
     Vector<double> s(DIM, 0.0);
  
     // Loop over each level (apart from the coarsest level; an interpolation
     // matrix and thus a restriction matrix is not needed here), with 0 being
     // the finest level and Nlevel-1 being the coarsest level
     for (unsigned level = 0; level < Nlevel - 1; level++)
     {
       // Assign values to a couple of variables to help out
       unsigned coarse_level = level + 1;
       unsigned fine_level = level;
  
       // Make a pointer to the mesh on the finer level and dynamic_cast
       // it as an object of the refineable mesh class
       RefineableMeshBase* ref_fine_mesh_pt = dynamic_cast<RefineableMeshBase*>(
         Mg_hierarchy[fine_level]->mg_bulk_mesh_pt());
  
       // Make a pointer to the mesh on the coarse level and dynamic_cast
       // it as an object of the refineable mesh class
       RefineableMeshBase* ref_coarse_mesh_pt =
         dynamic_cast<RefineableMeshBase*>(
           Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
  
       // Access information about the number of elements in the fine mesh
       // from the pointer to the fine mesh (to loop over the rows of the
       // interpolation matrix)
       unsigned fine_n_element = ref_fine_mesh_pt->nelement();
  
       // The numbers of rows and columns in the interpolation matrix. The
       // number of unknowns has been divided by 2 since there are 2 dofs at
       // each node in the mesh corresponding to the real and imaginary part
       // of the solution
       unsigned n_rows = Mg_hierarchy[fine_level]->ndof();
       unsigned n_cols = Mg_hierarchy[coarse_level]->ndof();
  
       // Mapping relating the pointers to related elements in the coarse and
       // fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]
       // If the element in the fine mesh has been unrefined between these two
       // levels, this map returns the father element in the coarsened mesh.
       // If this element in the fine mesh has not been unrefined,
       // the map returns the pointer to the same-sized equivalent
       // element in the coarsened mesh.
       std::map<RefineableQElement<DIM>*, RefineableQElement<DIM>*>
         coarse_mesh_reference_element_pt;
  
       // Counter of elements in coarse and fine meshes: Start with element
       // zero in both meshes.
       unsigned e_coarse = 0;
       unsigned e_fine = 0;
  
       // While loop over fine elements (while because we're
       // incrementing the counter internally if the element was
       // unrefined...)
       while (e_fine < fine_n_element)
       {
         // Pointer to element in fine mesh
         RefineableQElement<DIM>* el_fine_pt =
           dynamic_cast<RefineableQElement<DIM>*>(
             ref_fine_mesh_pt->finite_element_pt(e_fine));
  
         // Pointer to element in coarse mesh
         RefineableQElement<DIM>* el_coarse_pt =
           dynamic_cast<RefineableQElement<DIM>*>(
             ref_coarse_mesh_pt->finite_element_pt(e_coarse));
  
         // If the levels are different then the element in the fine
         // mesh has been unrefined between these two levels
         if (el_fine_pt->tree_pt()->level() != el_coarse_pt->tree_pt()->level())
         {
           // The element in the fine mesh has been unrefined between these two
           // levels. Hence it and its three brothers (ASSUMED to be stored
           // consecutively in the Mesh's vector of pointers to its constituent
           // elements -- we'll check this!) share the same father element in
           // the coarse mesh, currently pointed to by el_coarse_pt.
           for (unsigned i = 0; i < n_sons; i++)
           {
             // Set mapping to father element in coarse mesh
             coarse_mesh_reference_element_pt
               [dynamic_cast<RefineableQElement<DIM>*>(
                 ref_fine_mesh_pt->finite_element_pt(e_fine))] = el_coarse_pt;
  
             // Increment counter for elements in fine mesh
             e_fine++;
           }
         }
         // The element in the fine mesh has not been unrefined between
         // these two levels, so the reference element is the same-sized
         // equivalent element in the coarse mesh
         else
         {
           // Set the mapping between the two elements since they are
           // the same element
           coarse_mesh_reference_element_pt[el_fine_pt] = el_coarse_pt;
  
           // Increment counter for elements in fine mesh
           e_fine++;
         }
         // Increment counter for elements in coarse mesh
         e_coarse++;
       } // End of while loop for setting up element-coincidences
  
       // To allow update of a row only once we use stl vectors for bools
       std::vector<bool> contribution_made(n_rows, false);
  
       // Make storage vectors to form the interpolation matrix using a
       // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
       // row_start contains entries which tells us at which entry of the
       // vector column_start we start on the i-th row of the actual matrix.
       // The entries of column_start indicate the column position of each
       // non-zero entry in every row. This runs through the first row
       // from left to right then the second row (again, left to right)
       // and so on until the end. The entries in value are the entries in
       // the interpolation matrix. The vector column_start has the same length
       // as the vector value.
       Vector<double> value;
       Vector<int> column_index;
       Vector<int> row_start(n_rows + 1);
  
       // The value of index will tell us which row of the interpolation matrix
       // we're working on in the following for loop
       unsigned index = 0;
  
       // New loop to go over each element in the fine mesh
       for (unsigned k = 0; k < fine_n_element; k++)
       {
         // Set a pointer to the element in the fine mesh
         RefineableQElement<DIM>* el_fine_pt =
           dynamic_cast<RefineableQElement<DIM>*>(
             ref_fine_mesh_pt->finite_element_pt(k));
  
         // Get the reference element (either the father element or the
         // same-sized element) in the coarse mesh
         RefineableQElement<DIM>* el_coarse_pt =
           coarse_mesh_reference_element_pt[el_fine_pt];
  
         // Find out what type of son it is (set to OMEGA if no unrefinement
         // took place)
         int son_type = 0;
  
         // Check if the elements are different on both levels (i.e. to check
         // if any unrefinement took place)
         if (el_fine_pt->tree_pt()->level() != el_coarse_pt->tree_pt()->level())
         {
           // If there was refinement we need to find the son type
           son_type = el_fine_pt->tree_pt()->son_type();
         }
         else
         {
           // If there was no refinement then the son_type is given by the
           // value of Tree::OMEGA
           son_type = Tree::OMEGA;
         }
  
         // Find the number of nodes in the fine mesh element
         unsigned nnod_fine = el_fine_pt->nnode();
  
         // Loop through all the nodes in an element in the fine mesh
         for (unsigned i = 0; i < nnod_fine; i++)
         {
           // Row number in interpolation matrix: Global equation number
           // of the d.o.f. stored at this node in the fine element
           int ii = el_fine_pt->node_pt(i)->eqn_number(0);
  
           // Check whether or not the node is a proper d.o.f.
           if (ii >= 0)
           {
             // Only assign values to the given row of the interpolation
             // matrix if they haven't already been assigned
             if (contribution_made[ii] == false)
             {
               // The value of index was initialised when we allocated space
               // for the three vectors to store the CRDoubleMatrix. We use
               // index to go through the entries of the row_start vector.
               // The next row starts at the value.size() (draw out the entries
               // in value if this doesn't make sense noting that the storage
               // for the vector 'value' is dynamically allocated)
               row_start[index] = value.size();
  
               // Calculate the local coordinates of the given node
               el_fine_pt->local_coordinate_of_node(i, s);
  
               // Find the local coordinates s, of the present node, in the
               // reference element in the coarse mesh, given the element's
               // son_type (e.g. SW,SE... )
               level_up_local_coord_of_node(son_type, s);
  
               // Allocate space for shape functions in the coarse mesh
               Shape psi(el_coarse_pt->nnode());
  
               // Set the shape function in the reference element
               el_coarse_pt->shape(s, psi);
  
               // Auxiliary storage
               std::map<unsigned, double> contribution;
               Vector<unsigned> keys;
  
               // Find the number of nodes in an element in the coarse mesh
               unsigned nnod_coarse = el_coarse_pt->nnode();
  
               // Loop through all the nodes in the reference element
               for (unsigned j = 0; j < nnod_coarse; j++)
               {
                 // Column number in interpolation matrix: Global equation
                 // number of the d.o.f. stored at this node in the coarse
                 // element
                 int jj = el_coarse_pt->node_pt(j)->eqn_number(0);
  
                 // If the value stored at this node is pinned or hanging
                 if (jj < 0)
                 {
                   // Hanging node: In this case we need to accumulate the
                   // contributions from the master nodes
                   if (el_coarse_pt->node_pt(j)->is_hanging())
                   {
                     // Find the number of master nodes of the hanging
                     // the node in the reference element
                     HangInfo* hang_info_pt =
                       el_coarse_pt->node_pt(j)->hanging_pt();
                     unsigned nmaster = hang_info_pt->nmaster();
  
                     // Loop over the master nodes
                     for (unsigned i_master = 0; i_master < nmaster; i_master++)
                     {
                       // Set up a pointer to the master node
                       Node* master_node_pt =
                         hang_info_pt->master_node_pt(i_master);
  
                       // The column number in the interpolation matrix: the
                       // global equation number of the d.o.f. stored at this
                       // master node for the coarse element
                       int master_jj = master_node_pt->eqn_number(0);
  
                       // Is the master node a proper d.o.f.?
                       if (master_jj >= 0)
                       {
                         // If the weight of the master node is non-zero
                         if (psi(j) * hang_info_pt->master_weight(i_master) !=
                             0.0)
                         {
                           contribution[master_jj] +=
                             psi(j) * hang_info_pt->master_weight(i_master);
                         }
                       } // if (master_jj>=0)
                     } // for (unsigned i_master=0;i_master<nmaster;i_master++)
                   } // if (el_coarse_pt->node_pt(j)->is_hanging())
                 }
                 // In the case that the node is not pinned or hanging
                 else
                 {
                   // If we can get a nonzero contribution from the shape
                   // function
                   if (psi(j) != 0.0)
                   {
                     contribution[jj] += psi(j);
                   }
                 } // if (jj<0) else
               } // for (unsigned j=0;j<nnod_coarse;j++)
  
               // Put the contributions into the value vector
               for (std::map<unsigned, double>::iterator it =
                      contribution.begin();
                    it != contribution.end();
                    ++it)
               {
                 if (it->second != 0)
                 {
                   column_index.push_back(it->first);
                   value.push_back(it->second);
                 }
               } // for (std::map<unsigned,double>::iterator it=...)
  
               // Increment the value of index by 1 to indicate that we have
               // finished the row, or equivalently, covered another global
               // node in the fine mesh
               index++;
  
               // Change the entry in contribution_made to true now to indicate
               // that the row has been filled
               contribution_made[ii] = true;
             } // if(contribution_made[ii]==false)
           } // if (ii>=0)
         } // for(unsigned i=0;i<nnod_element;i++)
       } // for (unsigned k=0;k<fine_n_element;k++)
  
       // Set the last entry in the row_start vector
       row_start[n_rows] = value.size();
  
       // Set the interpolation matrix to be that formed as the CRDoubleMatrix
       // using the vectors value, row_start, column_index and the value
       // of fine_n_unknowns and coarse_n_unknowns
       interpolation_matrix_set(
         level, value, column_index, row_start, n_cols, n_rows);
     } // for (unsigned level=0;level<Nlevel-1;level++)
   } // End of setup_interpolation_matrices

References DIM, oomph::Data::eqn_number(), oomph::Mesh::finite_element_pt(), i, j, k, 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::MGSolver< DIM >::setup_interpolation_matrices_unstructured

Setup the interpolation matrices.

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

   {
     // Vector of local coordinates in the element
     Vector<double> s(DIM, 0.0);
  
     // Loop over each level (apart from the coarsest level; an interpolation
     // matrix and thus a restriction matrix is not needed here), with 0 being
     // the finest level and Nlevel-1 being the coarsest level
     for (unsigned level = 0; level < Nlevel - 1; level++)
     {
       // Assign values to a couple of variables to help out
       unsigned coarse_level = level + 1;
       unsigned fine_level = level;
  
       // Make a pointer to the mesh on the finer level and dynamic_cast
       // it as an object of the refineable mesh class
       Mesh* ref_fine_mesh_pt = Mg_hierarchy[fine_level]->mg_bulk_mesh_pt();
  
       // To use the locate zeta functionality the coarse mesh must be of the
       // type MeshAsGeomObject
       MeshAsGeomObject* coarse_mesh_from_obj_pt =
         new MeshAsGeomObject(Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
  
       // Access information about the number of degrees of freedom
       // from the pointers to the problem on each level
       unsigned coarse_n_unknowns = Mg_hierarchy[coarse_level]->ndof();
       unsigned fine_n_unknowns = Mg_hierarchy[fine_level]->ndof();
  
       // Make storage vectors to form the interpolation matrix using a
       // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
       // row_start contains entries which tells us at which entry of the
       // vector column_start we start on the i-th row of the actual matrix.
       // The entries of column_start indicate the column position of each
       // non-zero entry in every row. This runs through the first row
       // from left to right then the second row (again, left to right)
       // and so on until the end. The entries in value are the entries in
       // the interpolation matrix. The vector column_start has the same length
       // as the vector value.
       Vector<double> value;
       Vector<int> column_index;
       Vector<int> row_start(fine_n_unknowns + 1);
  
       // Vector to contain the (Eulerian) spatial location of the fine node
       Vector<double> fine_node_position(DIM);
  
       // Find the number of nodes in the mesh
       unsigned n_node_fine_mesh = ref_fine_mesh_pt->nnode();
  
       // Loop over the unknowns in the mesh
       for (unsigned i_fine_node = 0; i_fine_node < n_node_fine_mesh;
            i_fine_node++)
       {
         // Set a pointer to the i_fine_unknown-th node in the fine mesh
         Node* fine_node_pt = ref_fine_mesh_pt->node_pt(i_fine_node);
  
         // Get the global equation number
         int i_fine = fine_node_pt->eqn_number(0);
  
         // If the node is a proper d.o.f.
         if (i_fine >= 0)
         {
           // Row number in interpolation matrix: Global equation number
           // of the d.o.f. stored at this node in the fine element
           row_start[i_fine] = value.size();
  
           // Get the (Eulerian) spatial location of the fine node
           fine_node_pt->position(fine_node_position);
  
           // Create a null pointer to the GeomObject class
           GeomObject* el_pt = 0;
  
           // Get the reference element (either the father element or the
           // same-sized element) in the coarse mesh using locate_zeta
           coarse_mesh_from_obj_pt->locate_zeta(fine_node_position, el_pt, s);
  
           // Upcast GeomElement as a FiniteElement
           FiniteElement* el_coarse_pt = dynamic_cast<FiniteElement*>(el_pt);
  
           // Find the number of nodes in the element
           unsigned n_node = el_coarse_pt->nnode();
  
           // Allocate space for shape functions in the coarse mesh
           Shape psi(n_node);
  
           // Calculate the geometric shape functions at local coordinate s
           el_coarse_pt->shape(s, psi);
  
           // Auxiliary storage
           std::map<unsigned, double> contribution;
           Vector<unsigned> keys;
  
           // Loop through all the nodes in the (coarse mesh) element containing
           // the node pointed to by fine_node_pt (fine mesh)
           for (unsigned j_node = 0; j_node < n_node; j_node++)
           {
             // Get the j_coarse_unknown-th node in the coarse element
             Node* coarse_node_pt = el_coarse_pt->node_pt(j_node);
  
             // Column number in interpolation matrix: Global equation number of
             // the d.o.f. stored at this node in the coarse element
             int j_coarse = coarse_node_pt->eqn_number(0);
  
             // If the value stored at this node is pinned or hanging
             if (j_coarse < 0)
             {
               // Hanging node: In this case we need to accumulate the
               // contributions from the master nodes
               if (el_coarse_pt->node_pt(j_node)->is_hanging())
               {
                 // Find the number of master nodes of the hanging
                 // the node in the reference element
                 HangInfo* hang_info_pt = coarse_node_pt->hanging_pt();
                 unsigned nmaster = hang_info_pt->nmaster();
  
                 // Loop over the master nodes
                 for (unsigned i_master = 0; i_master < nmaster; i_master++)
                 {
                   // Set up a pointer to the master node
                   Node* master_node_pt = hang_info_pt->master_node_pt(i_master);
  
                   // The column number in the interpolation matrix: the
                   // global equation number of the d.o.f. stored at this master
                   // node for the coarse element
                   int master_jj = master_node_pt->eqn_number(0);
  
                   // Is the master node a proper d.o.f.?
                   if (master_jj >= 0)
                   {
                     // If the weight of the master node is non-zero
                     if (psi(j_node) * hang_info_pt->master_weight(i_master) !=
                         0.0)
                     {
                       contribution[master_jj] +=
                         psi(j_node) * hang_info_pt->master_weight(i_master);
                     }
                   } // End of if statement (check it's not a boundary node)
                 } // End of the loop over the master nodes
               } // End of the if statement for only hanging nodes
             } // End of the if statement for pinned or hanging nodes
             // In the case that the node is not pinned or hanging
             else
             {
               // If we can get a nonzero contribution from the shape function
               // at the j_node-th node in the element
               if (psi(j_node) != 0.0)
               {
                 contribution[j_coarse] += psi(j_node);
               }
             } // End of the if-else statement (check if the node was
               // pinned/hanging)
           } // Finished loop over the nodes j in the reference element (coarse)
  
           // Put the contributions into the value vector
           for (std::map<unsigned, double>::iterator it = contribution.begin();
                it != contribution.end();
                ++it)
           {
             if (it->second != 0)
             {
               value.push_back(it->second);
               column_index.push_back(it->first);
             }
           } // End of putting contributions into the value vector
         } // End check (whether or not the fine node was a d.o.f.)
       } // End of the for-loop over nodes in the fine mesh
  
       // Set the last entry of row_start
       row_start[fine_n_unknowns] = value.size();
  
       // Set the interpolation matrix to be that formed as the CRDoubleMatrix
       // using the vectors value, row_start, column_index and the value
       // of fine_n_unknowns and coarse_n_unknowns
       interpolation_matrix_set(level,
                                value,
                                column_index,
                                row_start,
                                coarse_n_unknowns,
                                fine_n_unknowns);
     } // End of loop over each level
   } // 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::Node::position(), s, oomph::FiniteElement::shape(), and Eigen::value.

◆ setup_mg_hierarchy()

template<unsigned DIM>

void oomph::MGSolver< DIM >::setup_mg_hierarchy

private

Set up the MG hierarchy.

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

   {
     // Initialise the timer start variable
     double t_m_start = 0.0;
  
     // Notify the user if it is allowed
     if (!Suppress_all_output)
     {
       // Notify user of progress
       oomph_info
         << "\n===============Creating Multigrid Hierarchy==============="
         << std::endl;
  
       // Start clock
       t_m_start = TimingHelpers::timer();
     }
  
     // Create a bool to indicate whether or not we could create an unrefined
     // copy. This bool will be assigned the value FALSE when the current copy
     // is the last level of the multigrid hierarchy
     bool managed_to_create_unrefined_copy = true;
  
     // Now keep making copies and try to make an unrefined copy of
     // the mesh
     unsigned level = 0;
  
     // Set up all of the levels by making a completely unrefined copy
     // of the problem using the function make_new_problem
     while (managed_to_create_unrefined_copy)
     {
       // Make a new object of the same type as the derived problem
       MGProblem* new_problem_pt = Mg_problem_pt->make_new_problem();
  
       // Do anything that needs to be done before we can refine the mesh
       new_problem_pt->actions_before_adapt();
  
       // To create the next level in the hierarchy we need to create a mesh
       // which matches the refinement of the current problem and then unrefine
       // the mesh. This can alternatively be done using the function
       // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
       // reference mesh to do precisely the above with
       managed_to_create_unrefined_copy =
         new_problem_pt->mg_bulk_mesh_pt()
           ->refine_base_mesh_as_in_reference_mesh_minus_one(
             Mg_hierarchy[level]->mg_bulk_mesh_pt());
  
       // If we were able to unrefine the problem on the current level
       // then add the unrefined problem to a vector of the levels
       if (managed_to_create_unrefined_copy)
       {
         // Another level has been created so increment the level counter
         level++;
  
         // If the documentation of everything has not been suppressed
         // then tell the user we managed to create another level
         if (!Suppress_all_output)
         {
           // Notify user that unrefinement was successful
           oomph_info << "\nSuccess! Level " << level << " has been created."
                      << std::endl;
         }
  
         // Do anything that needs to be done after refinement
         new_problem_pt->actions_after_adapt();
  
         // Do the equation numbering for the new problem
         oomph_info << "\n - Number of equations: "
                    << new_problem_pt->assign_eqn_numbers() << std::endl;
  
         // Add the new problem pointer onto the vector of MG levels
         // and increment the value of level by 1
         Mg_hierarchy.push_back(new_problem_pt);
       }
       // If we weren't able to create an unrefined copy
       else
       {
         // Delete the new problem
         delete new_problem_pt;
  
         // Make it a null pointer
         new_problem_pt = 0;
  
         // Assign the number of levels to Nlevel
         Nlevel = Mg_hierarchy.size();
  
         // If we're allowed to document then tell the user we've reached
         // the coarsest level of the hierarchy
         if (!Suppress_all_output)
         {
           // Notify the user
           oomph_info << "\n Reached the coarsest level! "
                      << "Number of levels: " << Nlevel << std::endl;
         }
       } // if (managed_to_unrefine)
     } // while (managed_to_unrefine)
  
     // Given that we know the number of levels in the hierarchy we can resize
     // the vectors which will store all the information required for our solver:
     // Resize the vector storing all of the system matrices
     Mg_matrices_storage_pt.resize(Nlevel, 0);
  
     // Resize the vector storing all of the solution vectors (X_mg)
     X_mg_vectors_storage.resize(Nlevel);
  
     // Resize the vector storing all of the RHS vectors (Rhs_mg)
     Rhs_mg_vectors_storage.resize(Nlevel);
  
     // Resize the vector storing all of the residual vectors
     Residual_mg_vectors_storage.resize(Nlevel);
  
     // Allocate space for the pre-smoother storage vector (remember, we do
     // not need a smoother on the coarsest level we use a direct solve there)
     Pre_smoothers_storage_pt.resize(Nlevel - 1, 0);
  
     // Allocate space for the post-smoother storage vector (remember, we do
     // not need a smoother on the coarsest level we use a direct solve there)
     Post_smoothers_storage_pt.resize(Nlevel - 1, 0);
  
     // Resize the vector storing all of the interpolation matrices
     Interpolation_matrices_storage_pt.resize(Nlevel - 1, 0);
  
     // Resize the vector storing all of the restriction matrices
     Restriction_matrices_storage_pt.resize(Nlevel - 1, 0);
  
     if (!Suppress_all_output)
     {
       // Stop clock
       double t_m_end = TimingHelpers::timer();
       double total_setup_time = double(t_m_end - t_m_start);
       oomph_info
         << "\nCPU time for creation of hierarchy of MG problems [sec]: "
         << total_setup_time << std::endl;
  
       // Notify user that the hierarchy of levels is complete
       oomph_info
         << "\n===============Hierarchy Creation Complete================"
         << "\n"
         << std::endl;
     }
   } // End of setup_mg_hierarchy

References oomph::Problem::actions_after_adapt(), oomph::Problem::actions_before_adapt(), oomph::Problem::assign_eqn_numbers(), oomph::MGProblem::make_new_problem(), oomph::MGProblem::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::MGSolver< DIM >::setup_mg_structures

private

Set up the MG hierarchy structures.

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

   {
     // Initialise the timer start variable
     double t_m_start = 0.0;
  
     // Start the clock (if we're allowed to time things)
     if (!Suppress_all_output)
     {
       // Start the clock
       t_m_start = TimingHelpers::timer();
     }
  
     // Allocate space for the system matrix on each level
     for (unsigned i = 0; i < Nlevel; i++)
     {
       // Dynamically allocate a new CRDoubleMatrix
       Mg_matrices_storage_pt[i] = new CRDoubleMatrix;
     }
  
     // Loop over each level and extract the system matrix, solution vector
     // right-hand side vector and residual vector (to store the value of r=b-Ax)
     for (unsigned i = 0; i < Nlevel; i++)
     {
       // If we're allowed to output
       if (!Suppress_all_output)
       {
         // Output the level we're working on
         oomph_info << "Setting up MG structures on level: " << i << "\n"
                    << std::endl;
       }
  
       // Resize the solution and RHS vector
       unsigned n_dof = Mg_hierarchy[i]->ndof();
       LinearAlgebraDistribution* dist_pt = new LinearAlgebraDistribution(
         Mg_hierarchy[i]->communicator_pt(), n_dof, false);
  
 #ifdef PARANOID
 #ifdef OOMPH_HAS_MPI
       // Set up the warning messages
       std::string warning_message = "Setup of distribution has not been ";
       warning_message += "tested with MPI.";
  
       // If we're not running the code in serial
       if (dist_pt->communicator_pt()->nproc() > 1)
       {
         // Throw a warning
         OomphLibWarning(
           warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
 #endif
 #endif
  
       // Build the approximate solution
       X_mg_vectors_storage[i].clear();
       X_mg_vectors_storage[i].build(dist_pt);
  
       // Build the point source function
       Rhs_mg_vectors_storage[i].clear();
       Rhs_mg_vectors_storage[i].build(dist_pt);
  
       // Build the residual vector (r=b-Ax)
       Residual_mg_vectors_storage[i].clear();
       Residual_mg_vectors_storage[i].build(dist_pt);
  
       // Delete the distribution pointer
       delete dist_pt;
  
       // Make it a null pointer
       dist_pt = 0;
  
       // Build the matrix distribution
       Mg_matrices_storage_pt[i]->clear();
       Mg_matrices_storage_pt[i]->distribution_pt()->build(
         Mg_hierarchy[i]->communicator_pt(), n_dof, false);
  
       // Compute system matrix on the current level. On the finest level of the
       // hierarchy the system matrix and RHS vector is given by the Jacobian and
       // vector of residuals which define the original problem which the
       // Galerkin approximation to the system matrix is used on the subsequent
       // levels so that the correct contributions are taken from each dof on the
       // level above (that is to say, it should match the contribution taken
       // from the solution vector and RHS vector on the level above)
       if (i == 0)
       {
         // Initialise the timer start variable
         double t_jac_start = 0.0;
  
         // If we're allowed to output things
         if (!Suppress_all_output)
         {
           // Start timer for Jacobian setup
           t_jac_start = TimingHelpers::timer();
         }
  
         // The system matrix on the finest level is the Jacobian and the RHS
         // vector is given by the residual vector which accompanies the Jacobian
         Mg_hierarchy[0]->get_jacobian(Rhs_mg_vectors_storage[0],
                                       *Mg_matrices_storage_pt[0]);
  
         if (!Suppress_all_output)
         {
           // Document the time taken
           double t_jac_end = TimingHelpers::timer();
           double jacobian_setup_time = t_jac_end - t_jac_start;
           oomph_info << " - Time for setup of Jacobian [sec]: "
                      << jacobian_setup_time << "\n"
                      << std::endl;
         }
       }
       else
       {
         // Initialise the timer start variable
         double t_gal_start = 0.0;
  
         // If we're allowed
         if (!Suppress_all_output)
         {
           // Start timer for Galerkin matrix calculation
           t_gal_start = TimingHelpers::timer();
         }
  
         // The system matrix on the coarser levels must be formed using the
         // Galerkin approximation which we do by calculating the product
         // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
         // finer grid system matrix is formed by multiplying by the (fine grid)
         // restriction matrix from the left and the (fine grid) interpolation
         // matrix from the left
  
         // First we need to calculate A^h * I^h_2h which we store as A^2h
         Mg_matrices_storage_pt[i - 1]->multiply(
           *Interpolation_matrices_storage_pt[i - 1],
           *Mg_matrices_storage_pt[i]);
  
         // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
         // was just calculated. This updates A^2h to give us the true
         // Galerkin approximation to the finer grid matrix
         Restriction_matrices_storage_pt[i - 1]->multiply(
           *Mg_matrices_storage_pt[i], *Mg_matrices_storage_pt[i]);
  
         // If the user did not choose to suppress everything
         if (!Suppress_all_output)
         {
           // End timer for Galerkin matrix calculation
           double t_gal_end = TimingHelpers::timer();
  
           // Calculate setup time
           double galerkin_matrix_calculation_time = t_gal_end - t_gal_start;
  
           // Document the time taken
           oomph_info
             << " - Time for system matrix formation using the Galerkin "
             << "approximation [sec]: " << galerkin_matrix_calculation_time
             << "\n"
             << std::endl;
         }
       } // if (i==0) else
     } // for (unsigned i=0;i<Nlevel;i++)
  
     // If we're allowed to output
     if (!Suppress_all_output)
     {
       // Stop clock
       double t_m_end = TimingHelpers::timer();
       double total_setup_time = double(t_m_end - t_m_start);
       oomph_info << "Total CPU time for setup of MG structures [sec]: "
                  << total_setup_time << std::endl;
  
       // Notify user that the hierarchy of levels is complete
       oomph_info << "\n============"
                  << "Multigrid Structures Setup Complete"
                  << "===========\n"
                  << std::endl;
     }
   } // End of setup_mg_structures

References oomph::LinearAlgebraDistribution::clear(), oomph::LinearAlgebraDistribution::communicator_pt(), i, oomph::OomphCommunicator::nproc(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::Global_string_for_annotation::string(), and oomph::TimingHelpers::timer().

◆ setup_smoothers()

template<unsigned DIM>

void oomph::MGSolver< 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

   {
     // Initialise the timer start variable
     double t_m_start = 0.0;
  
     // Start the clock (if we're allowed to time things)
     if (!Suppress_all_output)
     {
       // Notify user
       oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
  
       // Start the clock
       t_m_start = TimingHelpers::timer();
     }
  
     // Loop over the levels and assign the pre- and post-smoother points
     for (unsigned i = 0; i < Nlevel - 1; i++)
     {
       // If the pre-smoother factory function pointer hasn't been assigned
       // then we simply create a new instance of the DampedJacobi smoother
       // which is the default pre-smoother
       if (0 == Pre_smoother_factory_function_pt)
       {
         Pre_smoothers_storage_pt[i] = new DampedJacobi<CRDoubleMatrix>;
       }
       // Otherwise we use the pre-smoother factory function pointer to
       // generate a new pre-smoother
       else
       {
         // Get a pointer to an object of the same type as the pre-smoother
         Pre_smoothers_storage_pt[i] = (*Pre_smoother_factory_function_pt)();
       }
  
       // If the post-smoother factory function pointer hasn't been assigned
       // then we simply create a new instance of the DampedJacobi smoother
       // which is the default post-smoother
       if (0 == Post_smoother_factory_function_pt)
       {
         Post_smoothers_storage_pt[i] = new DampedJacobi<CRDoubleMatrix>;
       }
       // Otherwise we use the post-smoother factory function pointer to
       // generate a new post-smoother
       else
       {
         // Get a pointer to an object of the same type as the post-smoother
         Post_smoothers_storage_pt[i] = (*Post_smoother_factory_function_pt)();
       }
     }
  
     // Set the tolerance for the pre- and post-smoothers. The norm of the
     // solution which is compared against the tolerance is calculated
     // differently to the multigrid solver. To ensure that the smoother
     // continues to solve for the prescribed number of iterations we
     // lower the tolerance
     for (unsigned i = 0; i < Nlevel - 1; i++)
     {
       // Set the tolerance of the i-th level pre-smoother
       Pre_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
  
       // Set the tolerance of the i-th level post-smoother
       Post_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
     }
  
     // Set the number of pre- and post-smoothing iterations in each
     // pre- and post-smoother
     for (unsigned i = 0; i < Nlevel - 1; i++)
     {
       // Set the number of pre-smoothing iterations as the value of Npre_smooth
       Pre_smoothers_storage_pt[i]->max_iter() = Npre_smooth;
  
       // Set the number of pre-smoothing iterations as the value of Npost_smooth
       Post_smoothers_storage_pt[i]->max_iter() = Npost_smooth;
     }
  
     // Complete the setup of all of the smoothers
     for (unsigned i = 0; i < Nlevel - 1; i++)
     {
       // Pass a pointer to the system matrix on the i-th level to the i-th
       // level pre-smoother
       Pre_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
  
       // Pass a pointer to the system matrix on the i-th level to the i-th
       // level post-smoother
       Post_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
     }
  
     // Set up the distributions of each smoother
     for (unsigned i = 0; i < Nlevel - 1; i++)
     {
       // Get the number of dofs on the i-th level of the hierarchy.
       // This will be the same as the size of the solution vector
       // associated with the i-th level
       unsigned n_dof = X_mg_vectors_storage[i].nrow();
  
       // Create a LinearAlgebraDistribution which will be passed to the
       // linear solver
       LinearAlgebraDistribution dist(
         Mg_hierarchy[i]->communicator_pt(), n_dof, false);
  
 #ifdef PARANOID
 #ifdef OOMPH_HAS_MPI
       // Set up the warning messages
       std::string warning_message =
         "Setup of pre- and post-smoother distribution ";
       warning_message += "has not been tested with MPI.";
  
       // If we're not running the code in serial
       if (dist.communicator_pt()->nproc() > 1)
       {
         // Throw a warning
         OomphLibWarning(
           warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
       }
 #endif
 #endif
  
       // Build the distribution of the pre-smoother
       Pre_smoothers_storage_pt[i]->build_distribution(dist);
  
       // Build the distribution of the post-smoother
       Post_smoothers_storage_pt[i]->build_distribution(dist);
     }
  
     // Disable the smoother output on this level
     if (!Doc_time)
     {
       disable_smoother_and_superlu_doc_time();
     }
  
     // If we're allowed to output
     if (!Suppress_all_output)
     {
       // Stop clock
       double t_m_end = TimingHelpers::timer();
       double total_setup_time = double(t_m_end - t_m_start);
       oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
                  << total_setup_time << std::endl;
  
       // Notify user that the extraction is complete
       oomph_info
         << "\n==================Smoother Setup Complete================="
         << std::endl;
     }
   } // End of setup_smoothers

References 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(), and oomph::IterativeLinearSolver::tolerance().

◆ setup_transfer_matrices()

template<unsigned DIM>

void oomph::MGSolver< 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

   {
     // Initialise the timer start variable
     double t_r_start = 0.0;
  
     // Notify the user (if we're allowed)
     if (!Suppress_all_output)
     {
       // Notify user of progress
       oomph_info << "Creating the transfer matrices ";
  
       // Start the clock!
       t_r_start = TimingHelpers::timer();
     }
  
     // If we're allowed to output information
     if (!Suppress_all_output)
     {
       // Say what method we're using
       oomph_info << "using full weighting (recommended).\n" << std::endl;
     }
  
     // Using full weighting so use setup_interpolation_matrices.
     // Note: There are two methods to choose from here, the ideal choice is
     // setup_interpolation_matrices() but that requires a refineable mesh base
     if (dynamic_cast<TreeBasedRefineableMeshBase*>(
           Mg_problem_pt->mg_bulk_mesh_pt()))
     {
       setup_interpolation_matrices();
     }
     // If the mesh is unstructured we have to use the locate_zeta function
     // to set up the interpolation matrices
     else
     {
       setup_interpolation_matrices_unstructured();
     }
  
     // Loop over all levels that will be assigned a restriction matrix
     set_restriction_matrices_as_interpolation_transposes();
  
     // If we're allowed
     if (!Suppress_all_output)
     {
       // Stop the clock
       double t_r_end = TimingHelpers::timer();
       double total_G_setup_time = double(t_r_end - t_r_start);
       oomph_info << "CPU time for transfer matrices setup [sec]: "
                  << total_G_setup_time << std::endl;
  
       // Notify user that the hierarchy of levels is complete
       oomph_info
         << "\n============Transfer Matrices Setup Complete=============="
         << "\n"
         << std::endl;
     }
   } // End of setup_transfer_matrices function

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

◆ solve()

template<unsigned DIM>

void oomph::MGSolver< DIM >::solve	(	Problem *const &	problem_pt,
		DoubleVector &	result
	)

inlinevirtual

Virtual function in the base class that needs to be implemented later but for now just leave it empty

Implements oomph::LinearSolver.

     {
       // Dynamically cast problem_pt of type Problem to a MGProblem pointer
       MGProblem* mg_problem_pt = dynamic_cast<MGProblem*>(problem_pt);
  
 #ifdef PARANOID
       // PARANOID check - If the dynamic_cast produces a null pointer the
       // input was not a MGProblem
       if (0 == mg_problem_pt)
       {
         throw OomphLibError("Input problem must be of type MGProblem.",
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
       // PARANOID check - If a node in the input problem has more than
       // one value we cannot deal with it so arbitarily check the first one
       if (problem_pt->mesh_pt()->node_pt(0)->nvalue() != 1)
       {
         // How many dofs are there in the first node
         unsigned n_value = problem_pt->mesh_pt()->node_pt(0)->nvalue();
  
         // Make the error message
         std::ostringstream error_message_stream;
         error_message_stream << "Cannot currently deal with more than 1 dof"
                              << " per node. This problem has " << n_value
                              << " dofs per node." << std::endl;
  
         // Throw the error message
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // Assign the new MGProblem pointer to Mg_problem_pt
       Mg_problem_pt = mg_problem_pt;
  
       // Set up all of the required MG structures
       full_setup();
  
       // Run the MG method and assign the solution to result
       mg_solve(result);
  
       // Only output if the V-cycle output isn't suppressed
       if (!Suppress_v_cycle_output)
       {
         // Notify user that the hierarchy of levels is complete
         oomph_info << "\n================="
                    << "Multigrid Solve Complete"
                    << "=================\n"
                    << std::endl;
       }
  
       // If the user did not request all output be suppressed
       if (!Suppress_all_output)
       {
         // If the user requested all V-cycle output be suppressed
         if (Suppress_v_cycle_output)
         {
           // Create an extra line spacing
           oomph_info << std::endl;
         }
  
         // Output number of V-cycles taken to solve
         if (Has_been_solved)
         {
           oomph_info << "Total number of V-cycles required for solve: "
                      << V_cycle_counter << std::endl;
         }
         else
         {
           oomph_info << "Total number of V-cycles used: " << V_cycle_counter
                      << std::endl;
         }
       } // if (!Suppress_all_output)
  
       // Only enable and assign the stream pointer again if we originally
       // suppressed everything otherwise it won't be set yet
       if (Suppress_all_output)
       {
         // Now enable the stream pointer again
         oomph_info.stream_pt() = Stream_pt;
       }
     } // End of solve

Vector<DoubleVector> oomph::MGSolver< DIM >::X_mg_vectors_storage

private

Vector to store the solution vectors (X_mg)

Referenced by oomph::MGSolver< DIM >::direct_solve(), oomph::MGSolver< DIM >::post_smooth(), oomph::MGSolver< DIM >::pre_smooth(), and oomph::MGSolver< DIM >::residual_norm().

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

geometric_multigrid.h

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

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

Member Typedef Documentation

◆ PostSmootherFactoryFctPt

◆ PreSmootherFactoryFctPt

Constructor & Destructor Documentation

◆ MGSolver()

◆ ~MGSolver()

Member Function Documentation

◆ clean_up_memory()

◆ direct_solve()

◆ disable_output()

◆ disable_smoother_and_superlu_doc_time()

◆ disable_v_cycle_output()

◆ enable_doc_everything()

◆ enable_output()

◆ enable_v_cycle_output()

◆ full_setup()

◆ interpolate_and_correct()

◆ interpolation_matrix_set() [1/2]

◆ interpolation_matrix_set() [2/2]

◆ interpolation_self_test()

◆ iterations()

◆ level_up_local_coord_of_node() [1/3]

◆ level_up_local_coord_of_node() [2/3]

◆ level_up_local_coord_of_node() [3/3]

◆ max_iter()

◆ mg_solve()

◆ modify_restriction_matrices()

◆ npost_smooth()

◆ npre_smooth()

◆ plot()

◆ post_smooth()

◆ pre_smooth()

◆ residual_norm()

◆ restrict_residual()

◆ restriction_self_test()

◆ self_test()

◆ set_post_smoother_factory_function()

◆ set_pre_smoother_factory_function()

◆ set_restriction_matrices_as_interpolation_transposes()

◆ set_self_test_vector()

◆ setup_interpolation_matrices()

◆ setup_interpolation_matrices_unstructured()

◆ setup_mg_hierarchy()

◆ setup_mg_structures()

◆ setup_smoothers()

◆ setup_transfer_matrices()

◆ solve()

Member Data Documentation

◆ Doc_everything

◆ Has_been_setup

◆ Has_been_solved

◆ Interpolation_matrices_storage_pt

◆ Interpolation_self_test_vectors_storage

◆ Mg_hierarchy

◆ Mg_matrices_storage_pt

◆ Mg_problem_pt

◆ Nlevel

◆ Npost_smooth

◆ Npre_smooth

◆ Nvcycle

◆ Post_smoother_factory_function_pt

◆ Post_smoothers_storage_pt

◆ Pre_smoother_factory_function_pt

◆ Pre_smoothers_storage_pt

◆ Residual_mg_vectors_storage

◆ Restriction_matrices_storage_pt

◆ Restriction_self_test_vectors_storage

◆ Rhs_mg_vectors_storage

◆ Stream_pt

◆ Suppress_all_output

◆ Suppress_v_cycle_output

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