helmholtz_geometric_multigrid.h
Go to the documentation of this file.
1 // Include guards
2 #ifndef OOMPH_HELMHOLTZ_GEOMETRIC_MULTIGRID_HEADER
3 #define OOMPH_HELMHOLTZ_GEOMETRIC_MULTIGRID_HEADER
4 
5 // Config header generated by autoconfig
6 #ifdef HAVE_CONFIG_H
7 #include <config.h>
8 #endif
9 
10 // Oomph-lib headers
11 #include "generic/problem.h"
12 #include "generic/matrices.h"
13 #include "generic/preconditioner.h"
14 
15 // Include the complex smoother
16 #include "complex_smoother.h"
17 
18 // Namespace extension
19 namespace oomph
20 {
21 
22 //======================================================================
23 /// HelmholtzMGProblem class; subclass of Problem
24 //======================================================================
25  class HelmholtzMGProblem : public virtual Problem
26  {
27 
28  public:
29 
30  /// Constructor. Initialise pointers to coarser and finer levels
32  {}
33 
34  /// Destructor (empty)
36  {}
37 
38  /// \short This function needs to be implemented in the derived problem:
39  /// Returns a pointer to a new object of the same type as the derived
40  /// problem
42 
43  /// \short Function to get a pointer to the mesh we will be working
44  /// with. If there are flux elements present in the mesh this will
45  /// be overloaded to return a pointer to the bulk mesh otherwise
46  /// it can be overloaded to point to the global mesh but it must
47  /// be of type RefineableMeshBase
49 
50  }; // End of HelmholtzMGProblem class
51 
52 
53 //////////////////////////////////////////////////////////
54 //////////////////////////////////////////////////////////
55 //////////////////////////////////////////////////////////
56 
57 
58 //======================================================================
59 // MG solver class
60 //======================================================================
61  template<unsigned DIM>
62  class HelmholtzMGPreconditioner : public BlockPreconditioner<CRDoubleMatrix>
63  {
64 
65  public:
66 
67  /// \short typedef for a function that returns a pointer to an object
68  /// of the class HelmholtzSmoother to be used as the pre-smoother
69  typedef HelmholtzSmoother* (*PreSmootherFactoryFctPt)();
70 
71  /// \short typedef for a function that returns a pointer to an object
72  /// of the class HelmholtzSmoother to be used as the post-smoother
73  typedef HelmholtzSmoother* (*PostSmootherFactoryFctPt)();
74 
75  /// Access function to set the pre-smoother creation function.
77  PreSmootherFactoryFctPt pre_smoother_fn)
78  {
79  // Assign the function pointer
80  Pre_smoother_factory_function_pt=pre_smoother_fn;
81  }
82 
83  /// Access function to set the post-smoother creation function.
85  PostSmootherFactoryFctPt post_smoother_fn)
86  {
87  // Assign the function pointer
88  Post_smoother_factory_function_pt=post_smoother_fn;
89  }
90 
91  /// \short Constructor: Set up default values for number of V-cycles
92  /// and pre- and post-smoothing steps.
97  Mg_problem_pt(mg_problem_pt),
98  Tolerance(1.0e-09),
99  Npre_smooth(2),
100  Npost_smooth(2),
101  Nvcycle(2),
102  Doc_time(true),
104  Suppress_all_output(false),
105  Has_been_setup(false),
106  Has_been_solved(false),
107  Stream_pt(0),
108  Alpha_shift(0.0)
109  {
110  // Set the pointer to the finest level as the first
111  // entry in Mg_hierarchy_pt
112  Mg_hierarchy_pt.push_back(Mg_problem_pt);
113  } // End of HelmholtzMGPreconditioner
114 
115  /// Delete any dynamically allocated data
117  {
118  // Run the function written to clean up the memory
119  clean_up_memory();
120  } // End of ~HelmholtzMGPreconditioner
121 
122  /// Clean up anything that needs to be cleaned up
124  {
125  // We only need to destroy data if the solver has been set up and
126  // the data hasn't already been cleared
127  if (Has_been_setup)
128  {
129  // Loop over all of the levels in the hierarchy
130  for (unsigned i=0;i<Nlevel-1;i++)
131  {
132  // Delete the pre-smoother associated with this level
133  delete Pre_smoothers_storage_pt[i];
134 
135  // Make it a null pointer
137 
138  // Delete the post-smoother associated with this level
140 
141  // Make it a null pointer
143 
144  // Loop over the real and imaginary parts of the system matrix
145  // associated with the i-th level
146  for (unsigned j=0;j<2;j++)
147  {
148  // Delete the real/imaginary part of the system matrix
149  delete Mg_matrices_storage_pt[i][j];
150 
151  // Make it a null pointer
153  }
154  }
155 
156  // Delete the expanded matrix associated with the problem on the
157  // coarsest level
158  delete Coarsest_matrix_mg_pt;
159 
160  // Make it a null pointer
162 
163  // Loop over all but the coarsest of the levels in the hierarchy
164  for (unsigned i=0;i<Nlevel-1;i++)
165  {
166  // Delete the interpolation matrix associated with the i-th level
168 
169  // Make it a null pointer
171 
172  // Delete the restriction matrix associated with the i-th level
174 
175  // Make it a null pointer
177  }
178 
179  // Everything has been deleted now so we need to indicate that the
180  // solver is not set up
181  Has_been_setup=false;
182  }
183  } // End of clean_up_memory
184 
185  /// Access function for the variable Tolerance (lvalue)
186  double& tolerance()
187  {
188  // Return the variable, Tolerance
189  return Tolerance;
190  } // End of tolerance
191 
192  /// Function to change the value of the shift
193  double& alpha_shift()
194  {
195  // Return the alpha shift value
196  return Alpha_shift;
197  } // End of alpha_shift
198 
199  /// Disable time documentation
201  {
202  // Set the value of Doc_time to false
203  Doc_time=false;
204  } // End of disable_doc_time
205 
206  /// \short Disable all output from mg_solve apart from the number of
207  /// V-cycles used to solve the problem
209  {
210  // Set the value of Doc_time to false
211  Doc_time=false;
212 
213  // Enable the suppression of output from the V-cycle
215  } // End of disable_v_cycle_output
216 
217  /// \short Suppress anything that can be suppressed, i.e. any timings.
218  /// Things like mesh adaptation can not however be silenced using this
220  {
221  // Set the value of Doc_time to false
222  Doc_time=false;
223 
224  // Enable the suppression of output from the V-cycle
226 
227  // Enable the suppression of everything
228  Suppress_all_output=true;
229 
230  // Store the output stream pointer
232 
233  // Now set the oomph_info stream pointer to the null stream to
234  // disable all possible output
236  } // End of disable_output
237 
238  /// Enable time documentation
240  {
241  // Set the value of Doc_time to true
242  Doc_time=true;
243  } // End of enable_doc_time
244 
245  /// \short Enable the output of the V-cycle timings and other output
247  {
248  // Enable time documentation
249  Doc_time=true;
250 
251  // Enable output from the MG solver
253  } // End of enable_v_cycle_output
254 
255  /// \short Enable the output from anything that could have been suppressed
257  {
258  // Enable time documentation
259  Doc_time=true;
260 
261  // Enable output from everything during the full setup of the solver
262  Suppress_all_output=false;
263 
264  // Enable output from the MG solver
266  } // End of enable_output
267 
268  /// \short Suppress the output of both smoothers and SuperLU
270  {
271  // Loop over all levels of the hierarchy
272  for (unsigned i=0;i<Nlevel-1;i++)
273  {
274  // Disable time documentation on each level (for each pre-smoother)
275  Pre_smoothers_storage_pt[i]->disable_doc_time();
276 
277  // Disable time documentation on each level (for each post-smoother)
278  Post_smoothers_storage_pt[i]->disable_doc_time();
279  }
280 
281  // We only need a direct solve on the coarsest level so this is the
282  // only place we need to silence SuperLU
284  } // End of disable_smoother_and_superlu_doc_time
285 
286  /// Return the number of post-smoothing iterations (lvalue)
287  unsigned& npost_smooth()
288  {
289  // Return the number of post-smoothing iterations to be done on each
290  // level of the hierarchy
291  return Npost_smooth;
292  } // End of npost_smooth
293 
294  /// Return the number of pre-smoothing iterations (lvalue)
295  unsigned& npre_smooth()
296  {
297  // Return the number of pre-smoothing iterations to be done on each
298  // level of the hierarchy
299  return Npre_smooth;
300  } // End of npre_smooth
301 
302  /// \short Pre-smoother: Perform 'max_iter' smoothing steps on the
303  /// linear system Ax=b with current RHS vector, b, starting with
304  /// current solution vector, x. Return the residual vector r=b-Ax.
305  /// Uses the default smoother (set in the HelmholtzMGProblem constructor)
306  /// which can be overloaded for a specific problem.
307  void pre_smooth(const unsigned& level)
308  {
309  // HIERHER: Delete
310  // Calculate the residual norm
311  //double res_initial=residual_norm(level);
312  //double norm=1.0;
313 
314  // Output
315  //std::cout << "\nOn level: " << level
316  // << "\nResidual before: " << norm << std::endl;
317 
318  // Run pre-smoother 'max_iter' times
319  Pre_smoothers_storage_pt[level]->
320  complex_smoother_solve(Rhs_mg_vectors_storage[level],
321  X_mg_vectors_storage[level]);
322 
323  // Calculate the residual norm
324  //double res_new=residual_norm(level)/res_initial;
325 
326  // Output
327  //std::cout << "Residual after: " << res_new << std::endl;
328 
329  // Calculate the residual vector on this level
330  residual_norm(level);
331  } // End of pre_smooth
332 
333  /// \short Post-smoother: Perform max_iter smoothing steps on the
334  /// linear system Ax=b with current RHS vector, b, starting with
335  /// current solution vector, x. Uses the default smoother (set in
336  /// the HelmholtzMGProblem constructor) which can be overloaded for specific
337  /// problem.
338  void post_smooth(const unsigned& level)
339  {
340  // Run post-smoother 'max_iter' times
342  complex_smoother_solve(Rhs_mg_vectors_storage[level],
343  X_mg_vectors_storage[level]);
344 
345  // Calculate the residual vector on this level
346  residual_norm(level);
347  } // End of post_smooth
348 
349  /// \short Return norm of residual r=b-Ax and the residual vector itself
350  /// on the level-th level
351  double residual_norm(const unsigned& level)
352  {
353  // Loop over the real and imaginary part of the residual vector on
354  // the given level
355  for (unsigned j=0;j<2;j++)
356  {
357  // And zero the entries of residual
358  Residual_mg_vectors_storage[level][j].initialise(0.0);
359  }
360 
361  // Return the norm of the residual
362  return residual_norm(level,Residual_mg_vectors_storage[level]);
363  } // End of residual_norm
364 
365  /// Calculate the norm of the residual vector, r=b-Ax
366  double residual_norm(const unsigned& level,Vector<DoubleVector>& residual);
367 
368  /// \short Function to create the fully expanded system matrix on the
369  /// coarsest level
371 
372  /// \short Call the direct solver (SuperLU) to solve the problem exactly.
373  // The result is placed in X_mg
375  {
376  // Concatenate the vectors in X_mg into one vector, coarsest_x_mg
378  Coarsest_x_mg);
379 
380  // Concatenate the vectors in Rhs_mg into one vector, coarsest_rhs_mg
383 
384  // Get solution by direct solve:
386 
387  // Split the solution vector into a vector of DoubleVectors and store it
389  } // End of direct_solve
390 
391  /// \short Builds a CRDoubleMatrix that is used to interpolate the
392  /// residual between levels. The transpose can be used as the full
393  /// weighting restriction.
394  void interpolation_matrix_set(const unsigned& level,
395  double* value,
396  int* col_index,
397  int* row_st,
398  unsigned& ncol,
399  unsigned& nnz)
400  {
401  // Dynamically allocate the interpolation matrix
403 
404  // Build the matrix
406  build_without_copy(ncol,nnz,value,col_index,row_st);
407 
408  } // End of interpolation_matrix_set
409 
410  /// \short Builds a CRDoubleMatrix that is used to interpolate the
411  /// residual between levels. The transpose can be used as the full
412  /// weighting restriction.
413  void interpolation_matrix_set(const unsigned& level,
414  Vector<double>& value,
415  Vector<int>& col_index,
416  Vector<int>& row_st,
417  unsigned& ncol,
418  unsigned& nrow)
419  {
420  // Dynamically allocate the interpolation matrix
422 
423  // Make the distribution pointer
424  LinearAlgebraDistribution* dist_pt=
426  communicator_pt(),nrow,false);
427 
428 #ifdef PARANOID
429 #ifdef OOMPH_HAS_MPI
430  // Set up the warning messages
431  std::string warning_message="Setup of interpolation matrix distribution ";
432  warning_message+="has not been tested with MPI.";
433 
434  // If we're not running the code in serial
435  if (dist_pt->communicator_pt()->nproc()>1)
436  {
437  // Throw a warning
438  OomphLibWarning(warning_message,
439  OOMPH_CURRENT_FUNCTION,
440  OOMPH_EXCEPTION_LOCATION);
441  }
442 #endif
443 #endif
444 
445  // Build the matrix itself
447  build(dist_pt,ncol,value,col_index,row_st);
448 
449  // Delete the newly created distribution pointer
450  delete dist_pt;
451 
452  // Make it a null pointer
453  dist_pt=0;
454 
455  } // End of interpolation_matrix_set
456 
457  /// \short Builds a CRDoubleMatrix on each level that is used to
458  /// restrict the residual between levels. The transpose can be used
459  /// as the interpolation matrix
461  {
462  for (unsigned i=0;i<Nlevel-1;i++)
463  {
464  // Dynamically allocate the restriction matrix
466 
467  // Set the restriction matrix to be the transpose of the
468  // interpolation matrix
470  get_matrix_transpose(Restriction_matrices_storage_pt[i]);
471  }
472  } // End of set_restriction_matrices_as_interpolation_transposes
473 
474  /// \short Restrict residual (computed on level-th MG level) to the next
475  /// coarser mesh and stick it into the coarse mesh RHS vector.
476  void restrict_residual(const unsigned& level);
477 
478  /// \short Interpolate solution at current level onto next finer mesh
479  /// and correct the solution x at that level
480  void interpolate_and_correct(const unsigned& level);
481 
482  /// \short Given the son_type of an element and a local node number
483  /// j in that element with nnode_1d nodes per coordinate direction,
484  /// return the local coordinate s in its father element. Needed in
485  /// the setup of the interpolation matrices
486  void level_up_local_coord_of_node(const int& son_type,
487  Vector<double>& s);
488 
489  /// \short Setup the interpolation matrix on each level
491 
492  /// \short Setup the interpolation matrix on each level (used for
493  /// unstructured meshes)
495 
496  /// \short Setup the transfer matrices on each level
498 
499  /// \short Do a full setup (assumes everything will be setup around the
500  /// HelmholtzMGProblem pointer given in the constructor)
501  void full_setup();
502 
503  /// \short Function applies MG to the vector r for a full solve
505  {
506  // Split up the RHS vector into DoubleVectors, whose entries are
507  // arranged to match the matrix blocks and assign it
509 
510  // Split up the solution vector into DoubleVectors, whose entries are
511  // arranged to match the matrix blocks and assign it
513 
514  // Run the MG method and assign the solution to z
515  this->mg_solve(X_mg_vectors_storage[0]);
516 
517  // Copy solution in block vectors block_z back to z
519 
520  // Only output if the V-cycle output isn't suppressed
521  if (!(this->Suppress_v_cycle_output))
522  {
523  // Notify user that the hierarchy of levels is complete
524  oomph_info << "\n=========="
525  << "Multigrid Preconditioner Solve Complete"
526  << "=========" << std::endl;
527  }
528 
529  // Only enable and assign the stream pointer again if we originally
530  // suppressed everything otherwise it won't be set yet
531  if (this->Suppress_all_output)
532  {
533  // Now enable the stream pointer again
535  }
536  } // End of preconditioner_solve
537 
538  /// Number of iterations
539  unsigned iterations() const
540  {
541  // Return the number of V-cycles which have been done
542  return V_cycle_counter;
543  } // End of iterations
544 
545  /// \short Use the version in the Preconditioner base class for the
546  /// alternative setup function that takes a matrix pointer as an argument.
547  using Preconditioner::setup;
548 
549  private:
550 
551  /// Function to create pre-smoothers
553 
554  /// Function to create post-smoothers
556 
557  /// \short Do the actual solve -- this is called through the pure virtual
558  /// solve function in the LinearSolver base class. The function is stored
559  /// as protected to allow the HelmholtzMGPreconditioner derived class to use the
560  /// solver
561  void mg_solve(Vector<DoubleVector>& result);
562 
563  /// \short Function to ensure the block form of the Jacobian matches
564  /// the form described, i.e. we should have:
565  /// |-----|------|
566  /// | A_r | -A_c |
567  /// A = |-----|------|.
568  /// | A_c | A_r |
569  /// |-----|------|
571 
572  /// \short Function to set up the hierachy of levels. Creates a vector
573  /// of pointers to each MG level
574  void setup()
575  {
576  // Run the full setup
577  this->full_setup();
578 
579  // Only enable and assign the stream pointer again if we originally
580  // suppressed everything otherwise it won't be set yet
581  if (this->Suppress_all_output)
582  {
583  // Now enable the stream pointer again
585  }
586  } // End of setup
587 
588  /// \short Function to set up the hierachy of levels. Creates a vector
589  /// of pointers to each MG level
590  void setup_mg_hierarchy();
591 
592  /// \short Function to set up the hierachy of levels. Creates a vector
593  /// of pointers to each MG level
594  void setup_mg_structures();
595 
596  /// \short Estimate the value of the parameter h on the level-th problem
597  /// in the hierarchy.
598  void maximum_edge_width(const unsigned& level,double& h);
599 
600  /// \short Function to set up all of the smoothers once the system matrices
601  /// have been set up
602  void setup_smoothers();
603 
604  /// Pointer to the MG problem (deep copy)
606 
607  /// Vector containing pointers to problems in hierarchy
609 
610  /// \short Vector of vectors to store the system matrices. The i-th
611  /// entry in this vector contains a vector of length two. The first
612  /// entry of which contains the real part of the system matrix which
613  /// we refer to as A_r and the second entry contains the imaginary
614  /// part of the system matrix which we refer to as A_c. That is to say,
615  /// the true system matrix is given by A = A_r + iA_c
617 
618  /// \short Stores the system matrix on the coarest level in the fully
619  /// expanded format:
620  /// |-----|------|
621  /// | A_r | -A_c |
622  /// A = |-----|------|.
623  /// | A_c | A_r |
624  /// |-----|------|
625  /// Note: this is NOT the same as A = A_r + iA_c
627 
628  /// \short Assuming we're solving the system Ax=b, this vector will
629  /// contain the expanded solution vector on the coarsest level of the
630  /// heirarchy. This will have the form:
631  /// |-----|
632  /// | x_r |
633  /// x = |-----|.
634  /// | x_c |
635  /// |-----|
637 
638  /// \short Assuming we're solving the system Ax=b, this vector will
639  /// contain the expanded solution vector on the coarsest level of the
640  /// heirarchy. This will have the form:
641  /// |-----|
642  /// | b_r |
643  /// b = |-----|.
644  /// | b_c |
645  /// |-----|
647 
648  /// Vector to store the interpolation matrices
650 
651  /// Vector to store the restriction matrices
653 
654  /// \short Vector of vectors to store the solution vectors (X_mg) in two
655  /// parts; the real and imaginary. To access the real part of the solution
656  /// vector on the i-th level we need to access X_mg_vectors_storage[i][0]
657  /// while accessing X_mg_vectors_storage[i][1] will give us the
658  /// corresponding imaginary part
660 
661  /// \short Vector of vectors to store the RHS vectors. This uses the same
662  /// format as the X_mg_vectors_storage vector
664 
665  /// \short Vector to vectors to store the residual vectors. This uses
666  /// the same format as the X_mg_vectors_storage vector
668 
669  /// Vector to store the pre-smoothers
671 
672  /// Vector to store the post-smoothers
674 
675  /// \short Vector to storage the maximum edge width of each mesh. We only
676  /// need the maximum edge width on levels where we use a smoother to
677  /// determine the value of kh
679 
680  /// The value of the wavenumber, k
681  double Wavenumber;
682 
683  /// The tolerance of the multigrid preconditioner
684  double Tolerance;
685 
686  /// The number of levels in the multigrid heirachy
687  unsigned Nlevel;
688 
689  /// Number of pre-smoothing steps
690  unsigned Npre_smooth;
691 
692  /// Number of post-smoothing steps
693  unsigned Npost_smooth;
694 
695  /// Maximum number of V-cycles
696  unsigned Nvcycle;
697 
698  /// Pointer to counter for V-cycles
699  unsigned V_cycle_counter;
700 
701  /// Indicates whether or not time documentation should be used
702  bool Doc_time;
703 
704  /// Indicates whether or not the V-cycle output should be suppressed.
706 
707  /// If this is set to true then all output from the solver is suppressed
709 
710  /// Boolean variable to indicate whether or not the solver has been setup
712 
713  /// \short Boolean variable to indicate whether or not the problem was
714  /// successfully solved
716 
717  /// Pointer to the output stream -- defaults to std::cout
718  std::ostream* Stream_pt;
719 
720  /// Temporary version of the shift -- to run bash scripts
721  double Alpha_shift;
722  };
723 
724  //========================================================================
725  /// \short Calculating the residual r=b-Ax in the complex case requires
726  /// more care than the real case. To calculate the residual vector we
727  /// split A, x and b into their complex components:
728  /// r = b - A*x,
729  /// = (b_r + i*b_c) - (A_r + i*A_c)*(x_r + i*x_c),
730  /// = [b_r - A_r*x_r + A_c*x_c] + i*[b_c - A_r*x_c - A_c*x_r],
731  /// ==> real(r) = b_r - A_r*x_r + A_c*x_c,
732  /// & imag(r) = b_c - A_r*x_c - A_c*x_r.
733  //========================================================================
734  template<unsigned DIM>
736  residual_norm(const unsigned& level,Vector<DoubleVector>& residual)
737  {
738  // Number of rows in each block vector
739  unsigned n_rows=X_mg_vectors_storage[level][0].nrow();
740 
741 #ifdef PARANOID
742  // PARANOID check - if the residual vector doesn't have length 2 it cannot
743  // be used here since we need two vectors corresponding to the real and
744  // imaginary part of the residual
745  if (residual.size()!=2)
746  {
747  // Throw an error if the residual vector doesn't have the correct length
748  throw OomphLibError("This residual vector must have length 2. ",
749  OOMPH_CURRENT_FUNCTION,
750  OOMPH_EXCEPTION_LOCATION);
751  }
752  if (residual[0].nrow()!=residual[1].nrow())
753  {
754  // Create an output stream
755  std::ostringstream error_message_stream;
756 
757  // Store the error message
758  error_message_stream << "Residual[0] has length: " << residual[0].nrow()
759  << "\nResidual[1] has length: " << residual[1].nrow()
760  << "\nThis method requires that the constituent "
761  << "DoubleVectors in residual have the same length. "
762  << std::endl;
763 
764  // Throw an error
765  throw OomphLibError(error_message_stream.str(),
766  OOMPH_CURRENT_FUNCTION,
767  OOMPH_EXCEPTION_LOCATION);
768  }
769 #endif
770 
771  // Loop over the block vectors
772  for (unsigned i=0;i<2;i++)
773  {
774  // Start by setting the distribution of the residuals vector if
775  // it is not set up
776  if (!residual[i].distribution_built())
777  {
778  // Set up distribution
779  LinearAlgebraDistribution dist(Mg_hierarchy_pt[level]->
780  communicator_pt(),n_rows,false);
781 
782  // Build the distribution
783  residual[i].build(&dist,0.0);
784  }
785  // Otherwise just zero the entries of residual
786  else
787  {
788 #ifdef PARANOID
789  // PARANOID check - if the residuals are distributed then this method
790  // cannot be used, a distributed residuals can only be assembled by
791  // get_jacobian(...) for CRDoubleMatrices
792  if (residual[i].distributed())
793  {
794  throw OomphLibError
795  ("This method can only assemble a non-distributed residuals vector ",
796  OOMPH_CURRENT_FUNCTION,
797  OOMPH_EXCEPTION_LOCATION);
798  }
799 #endif
800 
801  // And zero the entries of residual
802  residual[i].initialise(0.0);
803  }
804  }
805 
806  // Store the pointer to the distribution of Matrix_real_pt (the same as
807  // Matrix_imag_pt presumably so we only need to use one)
808  LinearAlgebraDistribution* dist_pt=
809  Mg_matrices_storage_pt[level][0]->distribution_pt();
810 
811  // Create 4 temporary vectors to store the various matrix-vector products
812  // required. The appropriate combination has been appended to the name of
813  // the vector:
814  // Vector to store A_r*x_r:
815  DoubleVector temp_vec_rr(dist_pt,0.0);
816 
817  // Vector to store A_r*x_c:
818  DoubleVector temp_vec_rc(dist_pt,0.0);
819 
820  // Vector to store A_c*x_r:
821  DoubleVector temp_vec_cr(dist_pt,0.0);
822 
823  // Vector to store A_c*x_c:
824  DoubleVector temp_vec_cc(dist_pt,0.0);
825 
826  // We can't delete the distribution pointer because the Jacobian on the
827  // finest level is using it but we can make it a null pointer
828  dist_pt=0;
829 
830  // Calculate the different combinations of A*x (or A*e depending on the
831  // level in the heirarchy) in the complex case:
832  // A_r*x_r:
833  Mg_matrices_storage_pt[level][0]->
834  multiply(X_mg_vectors_storage[level][0],temp_vec_rr);
835 
836  // A_r*x_c:
837  Mg_matrices_storage_pt[level][0]->
838  multiply(X_mg_vectors_storage[level][1],temp_vec_rc);
839 
840  // A_c*x_r:
841  Mg_matrices_storage_pt[level][1]->
842  multiply(X_mg_vectors_storage[level][0],temp_vec_cr);
843 
844  // A_c*x_c:
845  Mg_matrices_storage_pt[level][1]->
846  multiply(X_mg_vectors_storage[level][1],temp_vec_cc);
847 
848  // Real part of the residual:
849  residual[0]=Rhs_mg_vectors_storage[level][0];
850  residual[0]-=temp_vec_rr;
851  residual[0]+=temp_vec_cc;
852 
853  // Imaginary part of the residual:
854  residual[1]=Rhs_mg_vectors_storage[level][1];
855  residual[1]-=temp_vec_rc;
856  residual[1]-=temp_vec_cr;
857 
858  // Get the residual norm of the real part of the residual vector
859  double norm_r=residual[0].norm();
860 
861  // Get the residual norm of the imaginary part of the residual vector
862  double norm_c=residual[1].norm();
863 
864  // Return the true norm of the residual vector which depends on the
865  // norm of the real part and the norm of the imaginary part
866  return sqrt(norm_r*norm_r+norm_c*norm_c);
867  }
868 
869  //=====================================================================
870  /// \short Check the block preconditioner framework returns the correct
871  /// system matrix
872  //=====================================================================
873  template<unsigned DIM>
875  {
876  // Start clock
877  clock_t t_bl_start=clock();
878 
879 #ifdef PARANOID
880  if (Mg_hierarchy_pt[0]->mesh_pt()==0)
881  {
882  std::stringstream err;
883  err << "Please set pointer to mesh using set_bulk_helmholtz_mesh(...).\n";
884  throw OomphLibError(err.str(),
885  OOMPH_CURRENT_FUNCTION,
886  OOMPH_EXCEPTION_LOCATION);
887  }
888 #endif
889 
890  // The preconditioner works with one mesh; set it! Since we only use the
891  // block preconditioner on the finest level, we use the mesh from that level
892  this->set_nmesh(1);
893 
894  // Elements in actual pml layer are trivially wrapped versions of
895  // their bulk counterparts. Technically they are different
896  // elements so we have to allow different element types
897  bool allow_different_element_types_in_mesh=true;
898  this->set_mesh(0,Mg_problem_pt->mesh_pt(),
899  allow_different_element_types_in_mesh);
900 
901 #ifdef PARANOID
902  // This preconditioner only works for 2 dof types
903  unsigned n_dof_types=this->ndof_types();
904  if (n_dof_types!=2)
905  {
906  std::stringstream tmp;
907  tmp << "This preconditioner only works for problems with 2 dof types\n"
908  << "Yours has " << n_dof_types;
909  throw OomphLibError(tmp.str(),
910  OOMPH_CURRENT_FUNCTION,
911  OOMPH_EXCEPTION_LOCATION);
912  }
913 #endif
914 
915  // Set up the generic block look up scheme
916  this->block_setup();
917 
918  // Extract the number of blocks.
919  unsigned nblock_types=this->nblock_types();
920 #ifdef PARANOID
921  if (nblock_types!=2)
922  {
923  std::stringstream tmp;
924  tmp << "There are supposed to be two block types.\n"
925  << "Yours has " << nblock_types;
926  throw OomphLibError(tmp.str(),
927  OOMPH_CURRENT_FUNCTION,
928  OOMPH_EXCEPTION_LOCATION);
929  }
930 #endif
931 
932  // This is how the 2x2 block matrices are extracted. We retain the sanity
933  // check (i.e. the diagonals are the same and the off-diagonals are negatives
934  // of each other in PARANOID mode. Otherwise we only extract 2 matrices
935  DenseMatrix<CRDoubleMatrix*> block_pt(nblock_types,nblock_types,0);
936  for (unsigned i=0;i<nblock_types;i++)
937  {
938  for (unsigned j=0;j<nblock_types;j++)
939  {
940  // we want...
941  block_pt(i,j)=new CRDoubleMatrix;
942  this->get_block(i,j,*block_pt(i,j));
943  }
944  }
945 
946  // Check that diagonal matrices are the same
947  //------------------------------------------
948  {
949  unsigned nnz1=block_pt(0,0)->nnz();
950  unsigned nnz2=block_pt(1,1)->nnz();
951  if (nnz1!=nnz2)
952  {
953  std::stringstream tmp;
954  tmp << "nnz of diagonal blocks doesn't match: "
955  << nnz1 << " != " << nnz2 << std::endl;
956  throw OomphLibError(tmp.str(),
957  OOMPH_CURRENT_FUNCTION,
958  OOMPH_EXCEPTION_LOCATION);
959  }
960  unsigned nrow1=block_pt(0,0)->nrow();
961  unsigned nrow2=block_pt(1,1)->nrow();
962  if (nrow1!=nrow2)
963  {
964  std::stringstream tmp;
965  tmp << "nrow of diagonal blocks doesn't match: "
966  << nrow1 << " != " << nrow2 << std::endl;
967  throw OomphLibError(tmp.str(),
968  OOMPH_CURRENT_FUNCTION,
969  OOMPH_EXCEPTION_LOCATION);
970  }
971 
972  // Check entries
973  bool fail=false;
974  double tol=1.0e-15;
975  std::stringstream tmp;
976 
977  // Check row starts
978  for (unsigned i=0;i<nrow1+1;i++)
979  {
980  if (block_pt(0,0)->row_start()[i]!=
981  block_pt(1,1)->row_start()[i])
982  {
983  fail=true;
984  tmp << "Row starts of diag matrices don't match for row " << i
985  << " : "
986  << block_pt(0,0)->row_start()[i] << " "
987  << block_pt(1,1)->row_start()[i] << " "
988  << std::endl;
989  }
990  }
991 
992  // Check values and column indices
993  for (unsigned i=0;i<nnz1;i++)
994  {
995  if (block_pt(0,0)->column_index()[i]!=
996  block_pt(1,1)->column_index()[i])
997  {
998  fail=true;
999  tmp << "Column of diag matrices indices don't match for entry " << i
1000  << " : "
1001  << block_pt(0,0)->column_index()[i] << " "
1002  << block_pt(1,1)->column_index()[i] << " "
1003  << std::endl;
1004  }
1005 
1006  if (fabs(block_pt(0,0)->value()[i]-
1007  block_pt(1,1)->value()[i])>tol)
1008  {
1009  fail=true;
1010  tmp << "Values of diag matrices don't match for entry " << i
1011  << " : "
1012  << block_pt(0,0)->value()[i] << " "
1013  << block_pt(1,1)->value()[i] << " "
1014  << std::endl;
1015  }
1016  }
1017  if (fail)
1018  {
1019  throw OomphLibError(tmp.str(),
1020  OOMPH_CURRENT_FUNCTION,
1021  OOMPH_EXCEPTION_LOCATION);
1022  }
1023  }
1024 
1025  // Check that off-diagonal matrices are negatives
1026  //-----------------------------------------------
1027  {
1028  unsigned nnz1=block_pt(0,1)->nnz();
1029  unsigned nnz2=block_pt(1,0)->nnz();
1030  if (nnz1!=nnz2)
1031  {
1032  std::stringstream tmp;
1033  tmp << "nnz of diagonal blocks doesn't match: "
1034  << nnz1 << " != " << nnz2 << std::endl;
1035  throw OomphLibError(tmp.str(),
1036  OOMPH_CURRENT_FUNCTION,
1037  OOMPH_EXCEPTION_LOCATION);
1038  }
1039  unsigned nrow1=block_pt(0,1)->nrow();
1040  unsigned nrow2=block_pt(1,0)->nrow();
1041  if (nrow1!=nrow2)
1042  {
1043  std::stringstream tmp;
1044  tmp << "nrow of off-diagonal blocks doesn't match: "
1045  << nrow1 << " != " << nrow2 << std::endl;
1046  throw OomphLibError(tmp.str(),
1047  OOMPH_CURRENT_FUNCTION,
1048  OOMPH_EXCEPTION_LOCATION);
1049  }
1050 
1051 
1052  // Check entries
1053  bool fail=false;
1054  double tol=1.0e-15;
1055  std::stringstream tmp;
1056 
1057  // Check row starts
1058  for (unsigned i=0;i<nrow1+1;i++)
1059  {
1060  if (block_pt(0,1)->row_start()[i]!=
1061  block_pt(1,0)->row_start()[i])
1062  {
1063  fail=true;
1064  tmp << "Row starts of off-diag matrices don't match for row " << i
1065  << " : "
1066  << block_pt(0,1)->row_start()[i] << " "
1067  << block_pt(1,0)->row_start()[i] << " "
1068  << std::endl;
1069  }
1070  }
1071 
1072  // Check values and column indices
1073  for (unsigned i=0;i<nnz1;i++)
1074  {
1075  if (block_pt(0,1)->column_index()[i]!=
1076  block_pt(1,0)->column_index()[i])
1077  {
1078  fail=true;
1079  tmp << "Column indices of off-diag matrices don't match for entry " << i
1080  << " : "
1081  << block_pt(0,1)->column_index()[i] << " "
1082  << block_pt(1,0)->column_index()[i] << " "
1083  << std::endl;
1084  }
1085 
1086  if (fabs(block_pt(0,1)->value()[i]+
1087  block_pt(1,0)->value()[i])>tol)
1088  {
1089  fail=true;
1090  tmp << "Values of off-diag matrices aren't negatives of "
1091  << "each other for entry " << i
1092  << " : "
1093  << block_pt(0,1)->value()[i] << " "
1094  << block_pt(1,0)->value()[i] << " "
1095  << std::endl;
1096  }
1097  }
1098  if (fail)
1099  {
1100  throw OomphLibError(tmp.str(),
1101  OOMPH_CURRENT_FUNCTION,
1102  OOMPH_EXCEPTION_LOCATION);
1103  }
1104  }
1105 
1106  // Clean up
1107  for (unsigned i=0;i<nblock_types;i++)
1108  {
1109  for (unsigned j=0;j<nblock_types;j++)
1110  {
1111  delete block_pt(i,j);
1112  block_pt(i,j)=0;
1113  }
1114  }
1115 
1116  // Stop clock
1117  clock_t t_bl_end=clock();
1118  double total_setup_time=double(t_bl_end-t_bl_start)/CLOCKS_PER_SEC;
1119  oomph_info << "CPU time for block preconditioner self-test [sec]: "
1120  << total_setup_time << "\n" << std::endl;
1121 
1122  }
1123 
1124 //===================================================================
1125 /// Runs a full setup of the MG solver
1126 //===================================================================
1127  template<unsigned DIM>
1129  {
1130  // Initialise the timer start variable
1131  double t_fs_start=0.0;
1132 
1133  // If we're allowed to output
1134  if (!Suppress_all_output)
1135  {
1136  // Start the timer
1137  t_fs_start=TimingHelpers::timer();
1138 
1139  // Notify user that the hierarchy of levels is complete
1140  oomph_info << "\n========Starting Setup of Multigrid Preconditioner========"
1141  << std::endl;
1142 
1143  // Notify user that the hierarchy of levels is complete
1144  oomph_info << "\nStarting the full setup of the Helmholtz multigrid solver."
1145  << std::endl;
1146  }
1147 
1148 #ifdef PARANOID
1149  // PARANOID check - Make sure the dimension of the solver matches the
1150  // dimension of the elements used in the problems mesh
1151  if (dynamic_cast<FiniteElement*>
1152  (Mg_problem_pt->mesh_pt()->element_pt(0))->dim()!=DIM)
1153  {
1154  // Create an error message
1155  std::string err_strng="The dimension of the elements used in the mesh ";
1156  err_strng+="does not match the dimension of the solver.";
1157 
1158  // Throw the error
1159  throw OomphLibError(err_strng,
1160  OOMPH_CURRENT_FUNCTION,
1161  OOMPH_EXCEPTION_LOCATION);
1162  }
1163 
1164  // PARANOID check - The elements of the bulk mesh must all be refineable
1165  // elements otherwise we cannot deal with this
1166  if (Mg_problem_pt->mg_bulk_mesh_pt()!=0)
1167  {
1168  // Find the number of elements in the bulk mesh
1169  unsigned n_elements=Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
1170 
1171  // Loop over the elements in the mesh and ensure that they are
1172  // all refineable elements
1173  for (unsigned el_counter=0;el_counter<n_elements;el_counter++)
1174  {
1175  // Upcast global mesh element to a refineable element
1176  RefineableElement* el_pt=dynamic_cast<RefineableElement*>
1177  (Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));
1178 
1179  // Check if the upcast worked or not; if el_pt is a null pointer the
1180  // element is not refineable
1181  if (el_pt==0)
1182  {
1183  // Create an output steam
1184  std::ostringstream error_message_stream;
1185 
1186  // Store the error message
1187  error_message_stream << "Element in bulk mesh could not be upcast to "
1188  << "a refineable element." << std::endl;
1189 
1190  // Throw an error
1191  throw OomphLibError(error_message_stream.str(),
1192  OOMPH_CURRENT_FUNCTION,
1193  OOMPH_EXCEPTION_LOCATION);
1194  }
1195  } // for (unsigned el_counter=0;el_counter<n_elements;el_counter++)
1196  }
1197  // If the provided bulk mesh pointer is a null pointer
1198  else
1199  {
1200  // Create an output steam
1201  std::ostringstream error_message_stream;
1202 
1203  // Store the error message
1204  error_message_stream << "The provided bulk mesh pointer is a null pointer. "
1205  << std::endl;
1206 
1207  // Throw an error
1208  throw OomphLibError(error_message_stream.str(),
1209  OOMPH_CURRENT_FUNCTION,
1210  OOMPH_EXCEPTION_LOCATION);
1211  } // if (Mg_problem_pt->mg_bulk_mesh_pt()!=0)
1212 #endif
1213 
1214  // If this is not the first Newton step then we will already have things
1215  // in storage. If this is the case, delete them
1216  clean_up_memory();
1217 
1218  // Resize the Mg_hierarchy_pt vector
1219  Mg_hierarchy_pt.resize(1,0);
1220 
1221  // Set the pointer to the finest level as the first entry in Mg_hierarchy_pt
1222  Mg_hierarchy_pt[0]=Mg_problem_pt;
1223 
1224  // Create the hierarchy of levels
1225  setup_mg_hierarchy();
1226 
1227  // Set up the interpolation and restriction matrices
1228  setup_transfer_matrices();
1229 
1230  // Set up the data structures on each level, i.e. the system matrix,
1231  // LHS and RHS vectors and smoothers
1232  setup_mg_structures();
1233 
1234  // Set up the smoothers on all of the levels
1235  setup_smoothers();
1236 
1237  // Loop over all of the coarser levels
1238  for (unsigned i=1;i<Nlevel;i++)
1239  {
1240  // Delete the i-th coarse-grid HelmholtzMGProblem object
1241  delete Mg_hierarchy_pt[i];
1242 
1243  // Set it to be a null pointer
1244  Mg_hierarchy_pt[i]=0;
1245  }
1246 
1247  // Indicate that the full setup has been completed
1248  Has_been_setup=true;
1249 
1250  // If we're allowed to output to the screen
1251  if (!Suppress_all_output)
1252  {
1253  // Output the time taken to complete the full setup
1254  double t_fs_end=TimingHelpers::timer();
1255  double full_setup_time=t_fs_end-t_fs_start;
1256 
1257  // Output the CPU time
1258  oomph_info << "\nCPU time for full setup [sec]: "
1259  << full_setup_time << std::endl;
1260 
1261  // Notify user that the hierarchy of levels is complete
1262  oomph_info << "\n==============="
1263  << "Multigrid Full Setup Complete"
1264  << "==============\n" << std::endl;
1265  } // if (!Suppress_all_output)
1266  } // End of full_setup
1267 
1268  //===================================================================
1269  /// \short Set up the MG hierarchy. Creates a vector of pointers to
1270  /// each MG level and resizes internal storage for multigrid data
1271  //===================================================================
1272  // Function to set up the hierachy of levels.
1273  template<unsigned DIM>
1275  {
1276  // Initialise the timer start variable
1277  double t_m_start=0.0;
1278 
1279  // Notify the user if it is allowed
1280  if (!Suppress_all_output)
1281  {
1282  // Notify user of progress
1283  oomph_info << "\n==============="
1284  << "Creating Multigrid Hierarchy"
1285  << "===============\n" << std::endl;
1286 
1287  // Start clock
1288  t_m_start=TimingHelpers::timer();
1289  }
1290 
1291  // Create a bool to indicate whether or not we could create an unrefined
1292  // copy. This bool will be assigned the value FALSE when the current copy
1293  // is the last level of the multigrid hierarchy
1294  bool managed_to_create_unrefined_copy=true;
1295 
1296  // Now keep making copies and try to make an unrefined copy of
1297  // the mesh
1298  unsigned level=0;
1299 
1300  // Set up all of the levels by making a completely unrefined copy
1301  // of the problem using the function make_new_problem
1302  while (managed_to_create_unrefined_copy)
1303  {
1304  // Make a new object of the same type as the derived problem
1305  HelmholtzMGProblem* new_problem_pt=Mg_problem_pt->make_new_problem();
1306 
1307  // Do anything that needs to be done before we can refine the mesh
1308  new_problem_pt->actions_before_adapt();
1309 
1310  // To create the next level in the hierarchy we need to create a mesh
1311  // which matches the refinement of the current problem and then unrefine
1312  // the mesh. This can alternatively be done using the function
1313  // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
1314  // reference mesh to do precisely the above with
1315  managed_to_create_unrefined_copy=
1316  new_problem_pt->mg_bulk_mesh_pt()->
1317  refine_base_mesh_as_in_reference_mesh_minus_one(
1318  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt());
1319 
1320  // If we were able to unrefine the problem on the current level
1321  // then add the unrefined problem to a vector of the levels
1322  if (managed_to_create_unrefined_copy)
1323  {
1324  // Another level has been created so increment the level counter
1325  level++;
1326 
1327  // If the documentation of everything has not been suppressed
1328  // then tell the user we managed to create another level
1329  if (!Suppress_all_output)
1330  {
1331  // Notify user that unrefinement was successful
1332  oomph_info << "Success! Level " << level
1333  << " has been created:" << std::endl;
1334  }
1335 
1336  // Do anything that needs to be done after refinement
1337  new_problem_pt->actions_after_adapt();
1338 
1339  // Do the equation numbering for the new problem
1340  oomph_info << "\n - Number of equations: "
1341  << new_problem_pt->assign_eqn_numbers()
1342  << "\n" << std::endl;
1343 
1344  // Add the new problem pointer onto the vector of MG levels
1345  // and increment the value of level by 1
1346  Mg_hierarchy_pt.push_back(new_problem_pt);
1347  }
1348  // If we weren't able to create an unrefined copy
1349  else
1350  {
1351  // Delete the new problem
1352  delete new_problem_pt;
1353 
1354  // Make it a null pointer
1355  new_problem_pt=0;
1356 
1357  // Assign the number of levels to Nlevel
1358  Nlevel=Mg_hierarchy_pt.size();
1359 
1360  // If we're allowed to document then tell the user we've reached
1361  // the coarsest level of the hierarchy
1362  if (!Suppress_all_output)
1363  {
1364  // Notify the user
1365  oomph_info << "Reached the coarsest level! "
1366  << "Number of levels: " << Nlevel << std::endl;
1367  }
1368  } // if (managed_to_unrefine)
1369  } // while (managed_to_unrefine)
1370 
1371  //------------------------------------------------------------------
1372  // Given that we know the number of levels in the hierarchy we can
1373  // resize the vectors which will store all the information required
1374  // for our solver:
1375  //------------------------------------------------------------------
1376  // Resize the vector storing all of the system matrices
1377  Mg_matrices_storage_pt.resize(Nlevel);
1378 
1379  // Resize the vector storing all of the solution vectors (X_mg)
1380  X_mg_vectors_storage.resize(Nlevel);
1381 
1382  // Resize the vector storing all of the RHS vectors (Rhs_mg)
1383  Rhs_mg_vectors_storage.resize(Nlevel);
1384 
1385  // Resize the vector storing all of the residual vectors
1386  Residual_mg_vectors_storage.resize(Nlevel);
1387 
1388  // Allocate space for the Max_edge_width vector, we only need the value
1389  // of h on the levels where we use a smoother
1390  Max_edge_width.resize(Nlevel-1,0.0);
1391 
1392  // Allocate space for the pre-smoother storage vector (remember, we do
1393  // not need a smoother on the coarsest level; we use a direct solve there)
1394  Pre_smoothers_storage_pt.resize(Nlevel-1,0);
1395 
1396  // Allocate space for the post-smoother storage vector (remember, we do
1397  // not need a smoother on the coarsest level; we use a direct solve there)
1398  Post_smoothers_storage_pt.resize(Nlevel-1,0);
1399 
1400  // Resize the vector storing all of the interpolation matrices
1401  Interpolation_matrices_storage_pt.resize(Nlevel-1,0);
1402 
1403  // Resize the vector storing all of the restriction matrices
1404  Restriction_matrices_storage_pt.resize(Nlevel-1,0);
1405 
1406  // If we're allowed to output information to the screen
1407  if (!Suppress_all_output)
1408  {
1409  // Stop clock
1410  double t_m_end=TimingHelpers::timer();
1411  double total_setup_time=double(t_m_end-t_m_start);
1412  oomph_info << "\nCPU time for creation of hierarchy of MG problems [sec]: "
1413  << total_setup_time << std::endl;
1414 
1415  // Notify user that the hierarchy of levels is complete
1416  oomph_info << "\n==============="
1417  << "Hierarchy Creation Complete"
1418  << "================\n" << std::endl;
1419  }
1420  } // End of setup_mg_hierarchy
1421 
1422  //===================================================================
1423  /// \short Set up the transfer matrices. Both the pure injection and
1424  /// full weighting method have been implemented here but it is highly
1425  /// recommended that full weighting is used in general. In both
1426  /// methods the transpose of the transfer matrix is used to transfer
1427  /// a vector back
1428  //===================================================================
1429  template<unsigned DIM>
1431  {
1432  // Initialise the timer start variable
1433  double t_r_start=0.0;
1434 
1435  // Notify the user (if we're allowed)
1436  if (!Suppress_all_output)
1437  {
1438  // Notify user of progress
1439  oomph_info << "Creating the transfer matrices." << std::endl;
1440 
1441  // Start the clock!
1442  t_r_start=TimingHelpers::timer();
1443  }
1444 
1445  // Using full weighting so use setup_interpolation_matrices.
1446  // Note: There are two methods to choose from here, the ideal choice is
1447  // setup_interpolation_matrices() but that requires a refineable mesh base
1448  if (dynamic_cast<TreeBasedRefineableMeshBase*>
1449  (Mg_problem_pt->mg_bulk_mesh_pt()))
1450  {
1451  setup_interpolation_matrices();
1452  }
1453  // If the mesh is unstructured we have to use the locate_zeta function
1454  // to set up the interpolation matrices
1455  else
1456  {
1457  setup_interpolation_matrices_unstructured();
1458  }
1459 
1460  // Loop over all levels that will be assigned a restriction matrix
1461  set_restriction_matrices_as_interpolation_transposes();
1462 
1463  // If we're allowed
1464  if (!Suppress_all_output)
1465  {
1466  // Stop the clock
1467  double t_r_end=TimingHelpers::timer();
1468  double total_G_setup_time=double(t_r_end-t_r_start);
1469  oomph_info << "\nCPU time for transfer matrices setup [sec]: "
1470  << total_G_setup_time << std::endl;
1471 
1472  // Notify user that the hierarchy of levels is complete
1473  oomph_info << "\n============Transfer Matrices Setup Complete=============="
1474  << "\n" << std::endl;
1475  }
1476  } // End of setup_transfer_matrices function
1477 
1478  //===================================================================
1479  /// \short Set up the MG structures on each level
1480  //===================================================================
1481  template<unsigned DIM>
1483  {
1484  // Initialise the timer start variable
1485  double t_m_start=0.0;
1486 
1487  // Start the clock (if we're allowed to time things)
1488  if (!Suppress_all_output)
1489  {
1490  // Start the clock
1491  t_m_start=TimingHelpers::timer();
1492  }
1493 
1494  // Calculate the wavenumber value:
1495  //--------------------------------
1496  // Reset the value of Wavenumber
1497  Wavenumber=0.0;
1498 
1499  // Upcast the first element in the bulk mesh
1500  PMLHelmholtzEquations<DIM>* pml_helmholtz_el_pt=
1501  dynamic_cast<PMLHelmholtzEquations<DIM>*>
1502  (Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(0));
1503 
1504  // Grab the k_squared value from the element pointer and square root.
1505  // Note, we assume the wavenumber is the same everywhere in the mesh
1506  // and it is also the same on every level.
1507  // hierher: This is true unless we wish to look at dispersion analysis
1508  Wavenumber=sqrt(pml_helmholtz_el_pt->k_squared());
1509 
1510  // We don't need the pointer anymore so make it a null pointer but don't
1511  // delete the underlying element data
1512  pml_helmholtz_el_pt=0;
1513 
1514  // Set up the system matrix and accompanying vectors on each level:
1515  //-----------------------------------------------------------------
1516  // Loop over each level and extract the system matrix, solution vector
1517  // right-hand side vector and residual vector (to store the value of r=b-Ax)
1518  for (unsigned i=0;i<Nlevel;i++)
1519  {
1520  // If we're allowed to output
1521  if (!Suppress_all_output)
1522  {
1523  // Output the level we're working on
1524  oomph_info << "Setting up MG structures on level: " << i
1525  << "\n" << std::endl;
1526  }
1527 
1528  // Resize the solution and RHS vector
1529  unsigned n_dof_per_block=Mg_hierarchy_pt[i]->ndof()/2;
1530 
1531  // Create the linear algebra distribution to build the vectors
1532  LinearAlgebraDistribution* dist_pt=
1533  new LinearAlgebraDistribution(Mg_hierarchy_pt[i]->communicator_pt(),
1534  n_dof_per_block,false);
1535 
1536 #ifdef PARANOID
1537 #ifdef OOMPH_HAS_MPI
1538  // Set up the warning messages
1539  std::string warning_message="Setup of distribution has not been ";
1540  warning_message+="tested with MPI.";
1541 
1542  // If we're not running the code in serial
1543  if (dist_pt->communicator_pt()->nproc()>1)
1544  {
1545  // Throw a warning
1546  OomphLibWarning(warning_message,
1547  OOMPH_CURRENT_FUNCTION,
1548  OOMPH_EXCEPTION_LOCATION);
1549  }
1550 #endif
1551 #endif
1552 
1553  // Create storage for the i-th level system matrix:
1554  //-------------------------------------------------
1555  // Resize the i-th entry of the matrix storage vector to contain two
1556  // CRDoubleMatrix pointers
1557  Mg_matrices_storage_pt[i].resize(2,0);
1558 
1559  // Loop over the real and imaginary part
1560  for (unsigned j=0;j<2;j++)
1561  {
1562  // Dynamically allocate a new CRDoubleMatrix
1563  Mg_matrices_storage_pt[i][j]=new CRDoubleMatrix;
1564  }
1565 
1566  // Build the matrix distribution (real part)
1567  Mg_matrices_storage_pt[i][0]->build(dist_pt);
1568 
1569  // Build the matrix distribution (imaginary part)
1570  Mg_matrices_storage_pt[i][1]->build(dist_pt);
1571 
1572  // Create storage for the i-th level solution vector:
1573  //---------------------------------------------------
1574  // Resize the i-th level solution vector to contain two DoubleVectors
1575  X_mg_vectors_storage[i].resize(2);
1576 
1577  // Build the approximate solution (real part)
1578  X_mg_vectors_storage[i][0].build(dist_pt);
1579 
1580  // Build the approximate solution (imaginary part)
1581  X_mg_vectors_storage[i][1].build(dist_pt);
1582 
1583  // Create storage for the i-th level RHS vector:
1584  //----------------------------------------------
1585  // Resize the i-th level RHS vector to contain two DoubleVectors
1586  Rhs_mg_vectors_storage[i].resize(2);
1587 
1588  // Build the point source function (real part)
1589  Rhs_mg_vectors_storage[i][0].build(dist_pt);
1590 
1591  // Build the point source function (imaginary part)
1592  Rhs_mg_vectors_storage[i][1].build(dist_pt);
1593 
1594  // Create storage for the i-th level residual vector:
1595  //---------------------------------------------------
1596  // Resize the i-th level residual vector to contain two DoubleVectors
1597  Residual_mg_vectors_storage[i].resize(2);
1598 
1599  // Build the residual vector, r=b-Ax (real part)
1600  Residual_mg_vectors_storage[i][0].build(dist_pt);
1601 
1602  // Build the residual vector, r=b-Ax (imaginary part)
1603  Residual_mg_vectors_storage[i][1].build(dist_pt);
1604 
1605  // Delete the distribution pointer
1606  delete dist_pt;
1607 
1608  // Make it a null pointer
1609  dist_pt=0;
1610 
1611  // Compute system matrix on the current level. On the finest level of the
1612  // hierarchy the system matrix is given by the complex-shifted Laplacian
1613  // preconditioner. On the subsequent levels the Galerkin approximation
1614  // is used to give us a coarse-grid representation of the problem
1615  if (i==0)
1616  {
1617  // Initialise the timer start variable
1618  double t_jac_start=0.0;
1619 
1620  // If we're allowed to output things
1621  if (!Suppress_all_output)
1622  {
1623  // Start timer for Jacobian setup
1624  t_jac_start=TimingHelpers::timer();
1625  }
1626 
1627 #ifdef PARANOID
1628  // Make sure the block preconditioner returns the correct sort of matrix
1629  block_preconditioner_self_test();
1630 #endif
1631 
1632  // The preconditioner works with one mesh; set it!
1633  this->set_nmesh(1);
1634 
1635  // Elements in actual pml layer are trivially wrapped versions of their
1636  // bulk counterparts. Technically they are different elements so we have
1637  // to allow different element types
1638  bool allow_different_element_types_in_mesh=true;
1639  this->set_mesh(0,Mg_hierarchy_pt[0]->mesh_pt(),
1640  allow_different_element_types_in_mesh);
1641 
1642 #ifdef PARANOID
1643  // Find the number of dof types we're dealing with
1644  unsigned n_dof_types=this->ndof_types();
1645 
1646  // This preconditioner only works for 2 dof types
1647  if (n_dof_types!=2)
1648  {
1649  // Create an error message
1650  std::stringstream tmp;
1651  tmp << "This preconditioner only works for problems with 2 dof types\n"
1652  << "Yours has " << n_dof_types;
1653 
1654  // Throw an error
1655  throw OomphLibError(tmp.str(),
1656  OOMPH_CURRENT_FUNCTION,
1657  OOMPH_EXCEPTION_LOCATION);
1658  }
1659 #endif
1660 
1661  // If we're not using a value of zero for the alpha shift
1662  if (Alpha_shift!=0.0)
1663  {
1664  //------------------------------------------------------------------
1665  // Set the damping in all of the PML elements to create the complex-
1666  // shifted Laplacian preconditioner:
1667  //------------------------------------------------------------------
1668  // Create a pointer which will contain the value of the shift given
1669  // by Alpha_shift
1670  double* alpha_shift_pt=new double(Alpha_shift);
1671 
1672  // Calculate the number of elements in the mesh
1673  unsigned n_element=Mg_hierarchy_pt[0]->mesh_pt()->nelement();
1674 
1675  // To grab the static variable used to set the value of alpha we first
1676  // need to get an element of type PMLHelmholtzEquations (we arbitrarily
1677  // chose the first element in the mesh)
1679  dynamic_cast<PMLHelmholtzEquations<DIM>*>
1680  (Mg_hierarchy_pt[0]->mesh_pt()->element_pt(0));
1681 
1682  // Now grab the pointer from the element
1683  static double* default_physical_constant_value_pt=el_pt->alpha_pt();
1684 
1685  // Loop over all of the elements
1686  for (unsigned i_el=0;i_el<n_element;i_el++)
1687  {
1688  // Upcast from a GeneralisedElement to a PmlHelmholtzElement
1690  dynamic_cast<PMLHelmholtzEquations<DIM>*>
1691  (Mg_hierarchy_pt[0]->mesh_pt()->element_pt(i_el));
1692 
1693  // Make the internal element alpha pointer point to Alpha_shift (the
1694  // chosen value of the shift to be applied to the problem)
1695  el_pt->alpha_pt()=alpha_shift_pt;
1696  }
1697 
1698  //---------------------------------------------------------------------
1699  // We want to solve the problem with the modified Jacobian but the
1700  // Preconditioner will have a handle to pointer which points to the
1701  // system matrix. If we wish to use the block preconditioner to extract
1702  // the appropriate blocks of the matrix we need to assign it. To avoid
1703  // losing the original system matrix we will create a temporary pointer
1704  // to it which will be reassigned after we have use the block
1705  // preconditioner to extract the blocks of the shifted matrix:
1706  //---------------------------------------------------------------------
1707  // Create a temporary pointer to the Jacobian
1708  CRDoubleMatrix* jacobian_pt=this->matrix_pt();
1709 
1710  // Create a new CRDoubleMatrix to hold the shifted Jacobian matrix
1711  CRDoubleMatrix* shifted_jacobian_pt=new CRDoubleMatrix;
1712 
1713  // Allocate space for the residuals
1714  DoubleVector residuals;
1715 
1716  // Get the residuals vector and the shifted Jacobian. Note, we do
1717  // not need to assign the residuals vector; we're simply using
1718  // MG as a preconditioner
1719  Mg_hierarchy_pt[0]->get_jacobian(residuals,*shifted_jacobian_pt);
1720 
1721  // Replace the current matrix used in Preconditioner by the new matrix
1722  this->set_matrix_pt(shifted_jacobian_pt);
1723 
1724  // Set up the generic block look up scheme
1725  this->block_setup();
1726 
1727  // Extract the number of blocks.
1728  unsigned nblock_types=this->nblock_types();
1729 
1730 #ifdef PARANOID
1731  // PARANOID check - there must only be two block types
1732  if (nblock_types!=2)
1733  {
1734  // Create the error message
1735  std::stringstream tmp;
1736  tmp << "There are supposed to be two block types.\nYours has "
1737  << nblock_types << std::endl;
1738 
1739  // Throw an error
1740  throw OomphLibError(tmp.str(),
1741  OOMPH_CURRENT_FUNCTION,
1742  OOMPH_EXCEPTION_LOCATION);
1743  }
1744 #endif
1745 
1746  // Store the level
1747  unsigned level=0;
1748 
1749  // Loop over the rows of the block matrix
1750  for (unsigned i_row=0;i_row<nblock_types;i_row++)
1751  {
1752  // Fix the column index
1753  unsigned j_col=0;
1754 
1755  // Extract the required blocks, i.e. the first column
1756  this->get_block(i_row,j_col,*Mg_matrices_storage_pt[level][i_row]);
1757  }
1758 
1759  // The blocks have been extracted from the shifted Jacobian therefore
1760  // we no longer have any use for it
1761  delete shifted_jacobian_pt;
1762 
1763  // Now the appropriate blocks have been extracted from the shifted
1764  // Jacobian we reset the matrix pointer in Preconditioner to the
1765  // original Jacobian so the linear solver isn't affected
1766  this->set_matrix_pt(jacobian_pt);
1767 
1768  //--------------------------------------------------------
1769  // Reassign the damping factor in all of the PML elements:
1770  //--------------------------------------------------------
1771  // Loop over all of the elements
1772  for (unsigned i_el=0;i_el<n_element;i_el++)
1773  {
1774  // Upcast from a GeneralisedElement to a PmlHelmholtzElement
1776  dynamic_cast<PMLHelmholtzEquations<DIM>*>
1777  (Mg_hierarchy_pt[0]->mesh_pt()->element_pt(i_el));
1778 
1779  // Set the value of alpha
1780  el_pt->alpha_pt()=default_physical_constant_value_pt;
1781  }
1782 
1783  // We've finished using the alpha_shift_pt pointer so delete it
1784  // as it was dynamically allocated
1785  delete alpha_shift_pt;
1786 
1787  // Make it a null pointer
1788  alpha_shift_pt=0;
1789  }
1790  // If the value of the shift is zero then we use the original
1791  // Jacobian matrix
1792  else
1793  {
1794  // The Jacobian has already been provided so now we just need to set
1795  // up the generic block look up scheme
1796  this->block_setup();
1797 
1798  // Extract the number of blocks.
1799  unsigned nblock_types=this->nblock_types();
1800 
1801 #ifdef PARANOID
1802  // If there are not only two block types we have a problem
1803  if (nblock_types!=2)
1804  {
1805  std::stringstream tmp;
1806  tmp << "There are supposed to be two block types.\n"
1807  << "Yours has " << nblock_types;
1808  throw OomphLibError(tmp.str(),
1809  OOMPH_CURRENT_FUNCTION,
1810  OOMPH_EXCEPTION_LOCATION);
1811  }
1812 #endif
1813 
1814  // Store the level
1815  unsigned level=0;
1816 
1817  // Loop over the rows of the block matrix
1818  for (unsigned i_row=0;i_row<nblock_types;i_row++)
1819  {
1820  // Fix the column index (since we only want the first column)
1821  unsigned j_col=0;
1822 
1823  // Extract the required blocks
1824  this->get_block(i_row,j_col,*Mg_matrices_storage_pt[level][i_row]);
1825  }
1826  } // if (Alpha_shift!=0.0) else
1827 
1828  if (!Suppress_all_output)
1829  {
1830  // Document the time taken
1831  double t_jac_end=TimingHelpers::timer();
1832  double jacobian_setup_time=t_jac_end-t_jac_start;
1833  oomph_info << " - Time for setup of Jacobian block matrix [sec]: "
1834  << jacobian_setup_time << "\n" << std::endl;
1835  }
1836  }
1837  // If we're not dealing with the finest level we're dealing with a
1838  // coarse-grid problem
1839  else
1840  {
1841  // Initialise the timer start variable
1842  double t_gal_start=0.0;
1843 
1844  // If we're allowed
1845  if (!Suppress_all_output)
1846  {
1847  // Start timer for Galerkin matrix calculation
1848  t_gal_start=TimingHelpers::timer();
1849  }
1850 
1851  //---------------------------------------------------------------------
1852  // The system matrix on the coarser levels must be formed using the
1853  // Galerkin approximation which we do by calculating the product
1854  // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
1855  // finer grid system matrix is formed by multiplying by the (fine grid)
1856  // restriction matrix from the left and the (fine grid) interpolation
1857  // matrix from the left. Fortunately, since the restriction and
1858  // interpolation acts on the real and imaginary parts separately we
1859  // can calculate the real and imaginary parts of the Galerkin
1860  // approximation separately.
1861  //---------------------------------------------------------------------
1862 
1863  //----------------------------------------------
1864  // Real component of the Galerkin approximation:
1865  //----------------------------------------------
1866  // First we need to calculate A^h * I^h_2h which we store as A^2h
1867  Mg_matrices_storage_pt[i-1][0]->
1868  multiply(*Interpolation_matrices_storage_pt[i-1],
1869  *Mg_matrices_storage_pt[i][0]);
1870 
1871  // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1872  // was just calculated. This updates A^2h to give us the true
1873  // Galerkin approximation to the finer grid matrix
1874  Restriction_matrices_storage_pt[i-1]->
1875  multiply(*Mg_matrices_storage_pt[i][0],
1876  *Mg_matrices_storage_pt[i][0]);
1877 
1878  //---------------------------------------------------
1879  // Imaginary component of the Galerkin approximation:
1880  //---------------------------------------------------
1881  // First we need to calculate A^h * I^h_2h which we store as A^2h
1882  Mg_matrices_storage_pt[i-1][1]->
1883  multiply(*Interpolation_matrices_storage_pt[i-1],
1884  *Mg_matrices_storage_pt[i][1]);
1885 
1886  // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1887  // was just calculated. This updates A^2h to give us the true
1888  // Galerkin approximation to the finer grid matrix
1889  Restriction_matrices_storage_pt[i-1]->
1890  multiply(*Mg_matrices_storage_pt[i][1],
1891  *Mg_matrices_storage_pt[i][1]);
1892 
1893  // If the user did not choose to suppress everything
1894  if (!Suppress_all_output)
1895  {
1896  // End timer for Galerkin matrix calculation
1897  double t_gal_end=TimingHelpers::timer();
1898 
1899  // Calculate setup time
1900  double galerkin_matrix_calculation_time=t_gal_end-t_gal_start;
1901 
1902  // Document the time taken
1903  oomph_info << " - Time for system block matrix formation using the "
1904  << "Galerkin approximation [sec]: "
1905  << galerkin_matrix_calculation_time << "\n" << std::endl;
1906  }
1907  } // if (i==0) else
1908 
1909  // We only
1910  if (i<Nlevel-1)
1911  {
1912  // Find the maximum edge width, h:
1913  //--------------------------------
1914  // Initialise the timer start variable
1915  double t_para_start=0.0;
1916 
1917  // If we're allowed to output things
1918  if (!Suppress_all_output)
1919  {
1920  // Start timer for parameter calculation on the i-th level
1921  t_para_start=TimingHelpers::timer();
1922  }
1923 
1924  // Calculate the i-th entry of Wavenumber and Max_edge_width
1925  maximum_edge_width(i,Max_edge_width[i]);
1926 
1927  // If the user did not choose to suppress everything
1928  if (!Suppress_all_output)
1929  {
1930  // End timer for Galerkin matrix calculation
1931  double t_para_end=TimingHelpers::timer();
1932 
1933  // Calculate setup time
1934  double parameter_calculation_time=t_para_end-t_para_start;
1935 
1936  // Document the time taken
1937  oomph_info << " - Time for maximum edge width calculation [sec]: "
1938  << parameter_calculation_time << "\n" << std::endl;
1939  }
1940  } // if (i<Nlevel-1)
1941  } // for (unsigned i=0;i<Nlevel;i++)
1942 
1943  // The last system matrix that needs to be setup is the fully expanded
1944  // version of the system matrix on the coarsest level. This is needed
1945  // for the direct solve on the coarsest level
1946  setup_coarsest_level_structures();
1947 
1948  // If we're allowed to output
1949  if (!Suppress_all_output)
1950  {
1951  // Stop clock
1952  double t_m_end=TimingHelpers::timer();
1953  double total_setup_time=double(t_m_end-t_m_start);
1954  oomph_info << "CPU time for setup of MG structures [sec]: "
1955  << total_setup_time << std::endl;
1956 
1957  // Notify user that the hierarchy of levels is complete
1958  oomph_info << "\n============"
1959  << "Multigrid Structures Setup Complete"
1960  << "===========\n" << std::endl;
1961  }
1962  } // End of setup_mg_structures
1963 
1964  //=========================================================================
1965  /// \short Function to set up structures on the coarsest level of the MG
1966  /// hierarchy. This includes setting up the CRDoubleMatrix version of the
1967  /// coarsest level system matrix. This would otherwise be stored as a
1968  /// vector of pointers to the constituent CRDoubleMatrix objects which
1969  /// has the form:
1970  /// |-----|
1971  /// | A_r |
1972  /// Matrix_mg_pt = |-----|
1973  /// | A_i |
1974  /// |-----|
1975  /// and we want to construct:
1976  /// |-----|------|
1977  /// | A_r | -A_c |
1978  /// Coarse_matrix_mg_pt = |-----|------|
1979  /// | A_c | A_r |
1980  /// |-----|------|
1981  /// Once this is done we have to set up the distributions of the vectors
1982  /// associated with Coarse_matrix_mg_pt
1983  //=========================================================================
1984  template<unsigned DIM>
1986  {
1987  // Start timer
1988  double t_cm_start=TimingHelpers::timer();
1989 
1990  //---------------------------------------------------
1991  // Extract information from the constituent matrices:
1992  //---------------------------------------------------
1993 
1994  // Grab the real and imaginary matrix parts from storage
1995  CRDoubleMatrix* real_matrix_pt=Mg_matrices_storage_pt[Nlevel-1][0];
1996  CRDoubleMatrix* imag_matrix_pt=Mg_matrices_storage_pt[Nlevel-1][1];
1997 
1998  // Number of nonzero entries in A_r
1999  unsigned nnz_r=real_matrix_pt->nnz();
2000  unsigned nnz_c=imag_matrix_pt->nnz();
2001 
2002  // Calculate the total number of rows (and thus columns) in the real matrix
2003  unsigned n_rows_r=real_matrix_pt->nrow();
2004 
2005  // Acquire access to the value, row_start and column_index arrays from
2006  // the real and imaginary portions of the full matrix.
2007  // Real part:
2008  const double* value_r_pt=real_matrix_pt->value();
2009  const int* row_start_r_pt=real_matrix_pt->row_start();
2010  const int* column_index_r_pt=real_matrix_pt->column_index();
2011 
2012  // Imaginary part:
2013  const double* value_c_pt=imag_matrix_pt->value();
2014  const int* row_start_c_pt=imag_matrix_pt->row_start();
2015  const int* column_index_c_pt=imag_matrix_pt->column_index();
2016 
2017 #ifdef PARANOID
2018  // PARANOID check - make sure the matrices have the same number of rows
2019  // to ensure they are compatible
2020  if (real_matrix_pt->nrow()!=imag_matrix_pt->nrow())
2021  {
2022  std::ostringstream error_message_stream;
2023  error_message_stream << "The real and imaginary matrices do not have "
2024  << "compatible sizes";
2025  throw OomphLibError(error_message_stream.str(),
2026  OOMPH_CURRENT_FUNCTION,
2027  OOMPH_EXCEPTION_LOCATION);
2028  }
2029  // PARANOID check - make sure the matrices have the same number of columns
2030  // to ensure they are compatible
2031  if (real_matrix_pt->ncol()!=imag_matrix_pt->ncol())
2032  {
2033  std::ostringstream error_message_stream;
2034  error_message_stream << "The real and imaginary matrices do not have "
2035  << "compatible sizes";
2036  throw OomphLibError(error_message_stream.str(),
2037  OOMPH_CURRENT_FUNCTION,
2038  OOMPH_EXCEPTION_LOCATION);
2039  }
2040 #endif
2041 
2042  // Calculate the total number of nonzeros in the full matrix
2043  unsigned nnz=2*(nnz_r+nnz_c);
2044 
2045  // Allocate space for the row_start and column_index vectors associated with
2046  // the complete matrix
2047  Vector<double> value(nnz,0.0);
2048  Vector<int> column_index(nnz,0);
2049  Vector<int> row_start(2*n_rows_r+1,0);
2050 
2051  //----------------------------
2052  // Build the row start vector:
2053  //----------------------------
2054 
2055  // Loop over the rows of the row_start vector. This is decomposed into
2056  // two parts. The first (n_rows_r+1) entries are found by computing
2057  // the entry-wise addition of row_start_r+row_start_c. The remaining
2058  // entries are almost the same as the first (n_rows_r+1). The only
2059  // distinction is that we need to shift the values of the entries by
2060  // the number of nonzeros in the top half of A. This is obvious from
2061  // observing the structure of the complete matrix.
2062 
2063  // Loop over the rows in the top half:
2064  for (unsigned i=0;i<n_rows_r;i++)
2065  {
2066  // Set the i-th entry in the row start vector
2067  row_start[i]=*(row_start_r_pt+i)+*(row_start_c_pt+i);
2068  }
2069 
2070  // Loop over the rows in the bottom half:
2071  for (unsigned i=n_rows_r;i<2*n_rows_r;i++)
2072  {
2073  // Set the i-th entry in the row start vector (bottom half)
2074  row_start[i]=row_start[i-n_rows_r]+(nnz_r+nnz_c);
2075  }
2076 
2077  // Set the last entry in the row start vector
2078  row_start[2*n_rows_r]=nnz;
2079 
2080  //-----------------------------------------
2081  // Build the column index and value vector:
2082  //-----------------------------------------
2083 
2084  // Variable to store the index of the nonzero
2085  unsigned i_nnz=0;
2086 
2087  // Loop over the top half of the complete matrix
2088  for (unsigned i=0;i<n_rows_r;i++)
2089  {
2090  // Calculate the number of nonzeros in the i-th row of A_r
2091  unsigned i_row_r_nnz=*(row_start_r_pt+i+1)-*(row_start_r_pt+i);
2092 
2093  // Calculate the number of nonzeros in the i-th row of A_c
2094  unsigned i_row_c_nnz=*(row_start_c_pt+i+1)-*(row_start_c_pt+i);
2095 
2096  // The index of the first entry in the i-th row of A_r
2097  unsigned i_first_entry_r=*(row_start_r_pt+i);
2098 
2099  // The index of the first entry in the i-th row of A_c
2100  unsigned i_first_entry_c=*(row_start_c_pt+i);
2101 
2102  // Loop over the number of nonzeros in the row associated with A_r
2103  for (unsigned j=0;j<i_row_r_nnz;j++)
2104  {
2105  // Assign the column index of the j-th entry in the i-th row of A_r
2106  // to the column_index vector
2107  column_index[i_nnz]=*(column_index_r_pt+i_first_entry_r+j);
2108 
2109  // Assign the corresponding entry in the value vector
2110  value[i_nnz]=*(value_r_pt+i_first_entry_r+j);
2111 
2112  // Increment the value of i_nnz
2113  i_nnz++;
2114  }
2115 
2116  // Loop over the number of nonzeros in the row associated with -A_c
2117  for (unsigned j=0;j<i_row_c_nnz;j++)
2118  {
2119  // Assign the column index of the j-th entry in the i-th row of -A_c
2120  // to the column_index vector
2121  column_index[i_nnz]=*(column_index_c_pt+i_first_entry_c+j)+n_rows_r;
2122 
2123  // Assign the corresponding entry in the value vector
2124  value[i_nnz]=-*(value_c_pt+i_first_entry_c+j);
2125 
2126  // Increment the value of i_nnz
2127  i_nnz++;
2128  }
2129  } // for (unsigned i=0;i<n_rows_r;i++)
2130 
2131  // Loop over the bottom half of the complete matrix
2132  for (unsigned i=n_rows_r;i<2*n_rows_r;i++)
2133  {
2134  // Calculate the number of nonzeros in row i of A_r
2135  unsigned i_row_r_nnz=
2136  *(row_start_r_pt+i-n_rows_r+1)-*(row_start_r_pt+i-n_rows_r);
2137 
2138  // Calculate the number of nonzeros in row i of A_c
2139  unsigned i_row_c_nnz=
2140  *(row_start_c_pt+i-n_rows_r+1)-*(row_start_c_pt+i-n_rows_r);
2141 
2142  // Get the index of the first entry in the i-th row of A_r
2143  unsigned i_first_entry_r=*(row_start_r_pt+i-n_rows_r);
2144 
2145  // Get the index of the first entry in the i-th row of A_c
2146  unsigned i_first_entry_c=*(row_start_c_pt+i-n_rows_r);
2147 
2148  // Loop over the number of nonzeros in the row associated with A_c
2149  for (unsigned j=0;j<i_row_c_nnz;j++)
2150  {
2151  // Assign the column index of the j-th entry in the i-th row of A_c
2152  // to the column_index vector
2153  column_index[i_nnz]=*(column_index_c_pt+i_first_entry_c+j);
2154 
2155  // Assign the corresponding entry in the value vector
2156  value[i_nnz]=*(value_c_pt+i_first_entry_c+j);
2157 
2158  // Increment the value of i_nnz
2159  i_nnz++;
2160  }
2161 
2162  // Loop over the number of nonzeros in the row associated with A_r
2163  for (unsigned j=0;j<i_row_r_nnz;j++)
2164  {
2165  // Assign the column index of the j-th entry in the i-th row of A_r
2166  // to the column_index vector
2167  column_index[i_nnz]=*(column_index_r_pt+i_first_entry_r+j)+n_rows_r;
2168 
2169  // Assign the corresponding entry in the value vector
2170  value[i_nnz]=*(value_r_pt+i_first_entry_r+j);
2171 
2172  // Increment the value of i_nnz
2173  i_nnz++;
2174  }
2175  } // for (unsigned i=n_rows_r;i<2*n_rows_r;i++)
2176 
2177  //-----------------------
2178  // Build the full matrix:
2179  //-----------------------
2180 
2181  // Allocate space for a CRDoubleMatrix
2182  Coarsest_matrix_mg_pt=new CRDoubleMatrix;
2183 
2184  // Make the distribution pointer
2185  LinearAlgebraDistribution* dist_pt=
2186  new LinearAlgebraDistribution(Mg_hierarchy_pt[Nlevel-1]->
2187  communicator_pt(),2*n_rows_r,false);
2188 
2189  // First, we need to build the matrix. Making use of its particular
2190  // structure we know that there are 2*n_rows_r columns in this matrix.
2191  // The remaining information has just been sorted out
2192  Coarsest_matrix_mg_pt->build(dist_pt,
2193  2*n_rows_r,
2194  value,
2195  column_index,
2196  row_start);
2197 
2198  // Build the distribution of associated solution vector
2199  Coarsest_x_mg.build(dist_pt);
2200 
2201  // Build the distribution of associated RHS vector
2202  Coarsest_rhs_mg.build(dist_pt);
2203 
2204  // Delete the associated distribution pointer
2205  delete dist_pt;
2206 
2207  // Summarise setup
2208  double t_cm_end=TimingHelpers::timer();
2209  double total_setup_time=double(t_cm_end-t_cm_start);
2210  oomph_info << " - Time for formation of the full matrix "
2211  << "on the coarsest level [sec]: "
2212  << total_setup_time << "\n" << std::endl;
2213  } // End of setup_coarsest_matrix_mg
2214 
2215 
2216 //==========================================================================
2217 /// \short Find the value of the parameters h on the level-th problem in
2218 /// the hierarchy. The value of h is determined by looping over each element
2219 /// in the mesh and calculating the length of each side and take the maximum
2220 /// value.Note, this is a heuristic calculation; if the mesh is nonuniform
2221 /// or adaptive refinement is used then the value of h, is not the same
2222 /// everywhere so we find the maximum edge width instead. If, however,
2223 /// uniform refinement is used on a uniform mesh (using quad elements) then
2224 /// this will return the correct value of h.
2225 ///
2226 /// This is the explicit template specialisation of the case DIM=2.
2227 //==========================================================================
2228  template<>
2230  maximum_edge_width(const unsigned& level,double& h)
2231  {
2232  // Create a pointer to the "bulk" mesh
2233  TreeBasedRefineableMeshBase* bulk_mesh_pt=
2234  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt();
2235 
2236  // Reset the value of h
2237  h=0.0;
2238 
2239  // Find out how many nodes there are along one edge of the first element.
2240  // We assume here that all elements have the same number of nodes
2241  unsigned nnode_1d=dynamic_cast<FiniteElement*>
2242  (bulk_mesh_pt->element_pt(0))->nnode_1d();
2243 
2244  // Sort out corner node IDs:
2245  //--------------------------
2246  // Initialise a vector to store local node IDs of the corners
2247  Vector<unsigned> corner_node_id(4,0);
2248 
2249  // Identify the local node ID of the first corner
2250  corner_node_id[0]=0;
2251 
2252  // Identify the local node ID of the second corner
2253  corner_node_id[1]=nnode_1d-1;
2254 
2255  // Identify the local node ID of the third corner
2256  corner_node_id[2]=nnode_1d*nnode_1d-1;
2257 
2258  // Identify the local node ID of the fourth corner
2259  corner_node_id[3]=nnode_1d*(nnode_1d-1);
2260 
2261  // Create storage for the nodal information:
2262  //------------------------------------------
2263  // Pointer to the first corner node on the j-th edge
2264  Node* first_node_pt=0;
2265 
2266  // Pointer to the second corner node on the j-th edge
2267  Node* second_node_pt=0;
2268 
2269  // Vector to store the (Eulerian) position of the first corner node
2270  // along a given edge of the element
2271  Vector<double> first_node_x(2,0.0);
2272 
2273  // Vector to store the (Eulerian) position of the second corner node
2274  // along a given edge of the element
2275  Vector<double> second_node_x(2,0.0);
2276 
2277  // Calculate h:
2278  //-------------
2279  // Find out how many elements there are in the bulk mesh
2280  unsigned n_element=bulk_mesh_pt->nelement();
2281 
2282  // Store a pointer which will point to a given element in the bulk mesh
2283  FiniteElement* el_pt=0;
2284 
2285  // Initialise a dummy variable to compare with h
2286  double test_h=0.0;
2287 
2288  // Store the number of edges in a 2D quad element
2289  unsigned n_edge=4;
2290 
2291  // Loop over all of the elements in the bulk mesh
2292  for (unsigned i=0;i<n_element;i++)
2293  {
2294  // Upcast the pointer to the i-th element to a FiniteElement pointer
2295  el_pt=dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(i));
2296 
2297  // Loop over the edges of the element
2298  for (unsigned j=0;j<n_edge;j++)
2299  {
2300  // Get the local node ID of the first corner node on the j-th edge
2301  first_node_pt=el_pt->node_pt(corner_node_id[j]);
2302 
2303  // Get the local node ID of the second corner node on the j-th edge
2304  second_node_pt=el_pt->node_pt(corner_node_id[(j+1)%4]);
2305 
2306  // Obtain the (Eulerian) position of the first corner node
2307  for (unsigned n=0;n<2;n++)
2308  {
2309  // Grab the n-th coordinate of this node
2310  first_node_x[n]=first_node_pt->x(n);
2311  }
2312 
2313  // Obtain the (Eulerian) position of the second corner node
2314  for (unsigned n=0;n<2;n++)
2315  {
2316  // Grab the n-th coordinate of this node
2317  second_node_x[n]=second_node_pt->x(n);
2318  }
2319 
2320  // Reset the value of test_h
2321  test_h=0.0;
2322 
2323  // Calculate the norm of (second_node_x-first_node_x)
2324  for (unsigned n=0;n<2;n++)
2325  {
2326  // Update the value of test_h
2327  test_h+=pow(second_node_x[n]-first_node_x[n],2.0);
2328  }
2329 
2330  // Square root the value of h
2331  test_h=sqrt(test_h);
2332 
2333  // Check if the value of test_h is greater than h
2334  if (test_h>h)
2335  {
2336  // If the value of test_h is greater than h then store it
2337  h=test_h;
2338  }
2339  } // for (unsigned j=0;j<n_edge;j++)
2340  } // for (unsigned i=0;i<n_element;i++)
2341  } // End of maximum_edge_width
2342 
2343 //==========================================================================
2344 /// \short Find the value of the parameters h on the level-th problem in
2345 /// the hierarchy. The value of h is determined by looping over each element
2346 /// in the mesh and calculating the length of each side and take the maximum
2347 /// value.Note, this is a heuristic calculation; if the mesh is nonuniform
2348 /// or adaptive refinement is used then the value of h, is not the same
2349 /// everywhere so we find the maximum edge width instead. If, however,
2350 /// uniform refinement is used on a uniform mesh (using quad elements) then
2351 /// this will return the correct value of h.
2352 ///
2353 /// This is the explicit template specialisation of the case DIM=3. The
2354 /// calculation of h is different here. In 2D we were able to loop over
2355 /// each pair of nodes in an anti-clockwise manner since the only node
2356 /// pairs were {(C0,C1),(C1,C2),(C2,C3),(C3,C0)} where CN denotes the N-th
2357 /// corner in the element. In 3D this method cannot be used since we have
2358 /// 12 edges to consider.
2359 //==========================================================================
2360  template<>
2362  double& h)
2363  {
2364  // Create a pointer to the "bulk" mesh
2365  TreeBasedRefineableMeshBase* bulk_mesh_pt=
2366  Mg_hierarchy_pt[level]->mg_bulk_mesh_pt();
2367 
2368  //--------------------------------
2369  // Find the maximum edge width, h:
2370  //--------------------------------
2371  // Reset the value of h
2372  h=0.0;
2373 
2374  // Find out how many nodes there are along one edge of the first element.
2375  // We assume here that all elements have the same number of nodes
2376  unsigned nnode_1d=dynamic_cast<FiniteElement*>
2377  (bulk_mesh_pt->element_pt(0))->nnode_1d();
2378 
2379  // Sort out corner node IDs:
2380  //--------------------------
2381  // Grab the octree pointer from the first element in the bulk mesh
2382  OcTree* octree_pt=dynamic_cast<RefineableQElement<3>*>
2383  (bulk_mesh_pt->element_pt(0))->octree_pt();
2384 
2385  // Initialise a vector to store local node IDs of the corners
2386  Vector<unsigned> corner_node_id(8,0);
2387 
2388  // Identify the local node ID of the first corner
2389  corner_node_id[0]=octree_pt->
2390  vertex_to_node_number(OcTreeNames::LDB,nnode_1d);
2391 
2392  // Identify the local node ID of the second corner
2393  corner_node_id[1]=octree_pt->
2394  vertex_to_node_number(OcTreeNames::RDB,nnode_1d);
2395 
2396  // Identify the local node ID of the third corner
2397  corner_node_id[2]=octree_pt->
2398  vertex_to_node_number(OcTreeNames::LUB,nnode_1d);
2399 
2400  // Identify the local node ID of the fourth corner
2401  corner_node_id[3]=octree_pt->
2402  vertex_to_node_number(OcTreeNames::RUB,nnode_1d);
2403 
2404  // Identify the local node ID of the fifth corner
2405  corner_node_id[4]=octree_pt->
2406  vertex_to_node_number(OcTreeNames::LDF,nnode_1d);
2407 
2408  // Identify the local node ID of the sixth corner
2409  corner_node_id[5]=octree_pt->
2410  vertex_to_node_number(OcTreeNames::RDF,nnode_1d);
2411 
2412  // Identify the local node ID of the seventh corner
2413  corner_node_id[6]=octree_pt->
2414  vertex_to_node_number(OcTreeNames::LUF,nnode_1d);
2415 
2416  // Identify the local node ID of the eighth corner
2417  corner_node_id[7]=octree_pt->
2418  vertex_to_node_number(OcTreeNames::RUF,nnode_1d);
2419 
2420  // Sort out the local node ID pairs:
2421  //----------------------------------
2422  // Store the number of edges in a 3D cubic element
2423  unsigned n_edge=12;
2424 
2425  // Create a vector to store the pairs of adjacent nodes
2426  Vector<std::pair<unsigned,unsigned> > edge_node_pair(n_edge);
2427 
2428  // First edge
2429  edge_node_pair[0]=std::make_pair(corner_node_id[0],corner_node_id[1]);
2430 
2431  // Second edge
2432  edge_node_pair[1]=std::make_pair(corner_node_id[0],corner_node_id[2]);
2433 
2434  // Third edge
2435  edge_node_pair[2]=std::make_pair(corner_node_id[0],corner_node_id[4]);
2436 
2437  // Fourth edge
2438  edge_node_pair[3]=std::make_pair(corner_node_id[1],corner_node_id[3]);
2439 
2440  // Fifth edge
2441  edge_node_pair[4]=std::make_pair(corner_node_id[1],corner_node_id[5]);
2442 
2443  // Sixth edge
2444  edge_node_pair[5]=std::make_pair(corner_node_id[2],corner_node_id[3]);
2445 
2446  // Seventh edge
2447  edge_node_pair[6]=std::make_pair(corner_node_id[2],corner_node_id[6]);
2448 
2449  // Eighth edge
2450  edge_node_pair[7]=std::make_pair(corner_node_id[3],corner_node_id[7]);
2451 
2452  // Ninth edge
2453  edge_node_pair[8]=std::make_pair(corner_node_id[4],corner_node_id[5]);
2454 
2455  // Tenth edge
2456  edge_node_pair[9]=std::make_pair(corner_node_id[4],corner_node_id[6]);
2457 
2458  // Eleventh edge
2459  edge_node_pair[10]=std::make_pair(corner_node_id[5],corner_node_id[7]);
2460 
2461  // Twelfth edge
2462  edge_node_pair[11]=std::make_pair(corner_node_id[6],corner_node_id[7]);
2463 
2464  // Create storage for the nodal information:
2465  //------------------------------------------
2466  // Pointer to the first corner node on the j-th edge
2467  Node* first_node_pt=0;
2468 
2469  // Pointer to the second corner node on the j-th edge
2470  Node* second_node_pt=0;
2471 
2472  // Vector to store the (Eulerian) position of the first corner node
2473  // along a given edge of the element
2474  Vector<double> first_node_x(3,0.0);
2475 
2476  // Vector to store the (Eulerian) position of the second corner node
2477  // along a given edge of the element
2478  Vector<double> second_node_x(3,0.0);
2479 
2480  // Calculate h:
2481  //-------------
2482  // Find out how many elements there are in the bulk mesh
2483  unsigned n_element=bulk_mesh_pt->nelement();
2484 
2485  // Store a pointer which will point to a given element in the bulk mesh
2486  FiniteElement* el_pt=0;
2487 
2488  // Initialise a dummy variable to compare with h
2489  double test_h=0.0;
2490 
2491  // Loop over all of the elements in the bulk mesh
2492  for (unsigned i=0;i<n_element;i++)
2493  {
2494  // Upcast the pointer to the i-th element to a FiniteElement pointer
2495  el_pt=dynamic_cast<FiniteElement*>(bulk_mesh_pt->element_pt(i));
2496 
2497  // Loop over the edges of the element
2498  for (unsigned j=0;j<n_edge;j++)
2499  {
2500  // Get the local node ID of the first corner node on the j-th edge
2501  first_node_pt=el_pt->node_pt(edge_node_pair[j].first);
2502 
2503  // Get the local node ID of the second corner node on the j-th edge
2504  second_node_pt=el_pt->node_pt(edge_node_pair[j].second);
2505 
2506  // Obtain the (Eulerian) position of the first corner node
2507  for (unsigned n=0;n<3;n++)
2508  {
2509  // Grab the n-th coordinate of this node
2510  first_node_x[n]=first_node_pt->x(n);
2511  }
2512 
2513  // Obtain the (Eulerian) position of the second corner node
2514  for (unsigned n=0;n<3;n++)
2515  {
2516  // Grab the n-th coordinate of this node
2517  second_node_x[n]=second_node_pt->x(n);
2518  }
2519 
2520  // Reset the value of test_h
2521  test_h=0.0;
2522 
2523  // Calculate the norm of (second_node_x-first_node_x)
2524  for (unsigned n=0;n<3;n++)
2525  {
2526  // Update the value of test_h
2527  test_h+=pow(second_node_x[n]-first_node_x[n],2.0);
2528  }
2529 
2530  // Square root the value of h
2531  test_h=sqrt(test_h);
2532 
2533  // Check if the value of test_h is greater than h
2534  if (test_h>h)
2535  {
2536  // If the value of test_h is greater than h then store it
2537  h=test_h;
2538  }
2539  } // for (unsigned j=0;j<n_edge;j++)
2540  } // for (unsigned i=0;i<n_element;i++)
2541  } // End of maximum_edge_width
2542 
2543 //=========================================================================
2544 /// \short Set up the smoothers on all levels
2545 //=========================================================================
2546  template<unsigned DIM>
2548  {
2549  // Initialise the timer start variable
2550  double t_m_start=0.0;
2551 
2552  // Start the clock (if we're allowed to time things)
2553  if (!Suppress_all_output)
2554  {
2555  // Notify user
2556  oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
2557 
2558  // Start the clock
2559  t_m_start=TimingHelpers::timer();
2560  }
2561 
2562 
2563  //===================================================TO_BE_DELETED=========
2564  // Temporary vector
2565  Vector<double> edge_widths(Nlevel,0.0);
2566 
2567  // Create storage for the wavenumber value
2568  double k=Wavenumber;
2569 
2570  // Loop over the levels and assign the pre- and post-smoother pointers
2571  for (unsigned i=0;i<Nlevel;i++)
2572  {
2573  // Calculate the i-th entry of Wavenumber and Max_edge_width
2574  maximum_edge_width(i,edge_widths[i]);
2575 
2576  // Create storage for the interval spacing
2577  double h=edge_widths[i];
2578 
2579  // Estimate the value of kh and output
2580  double kh=k*h;
2581  double omega=(12.0-4.0*pow(kh,2.0))/(18.0-3.0*pow(kh,2.0));
2582  double pi=MathematicalConstants::Pi;
2583  double omega_2=(2.0-pow(kh,2.0))/(2.0*pow(sin(pi*h/2),2.0)-0.5*pow(kh,2.0));
2584  oomph_info << "--------------------------"
2585  << "\nTable of estimated values:"
2586  << "\n--------------------------"
2587  << "\nk: " << k
2588  << "\nh: " << h
2589  << "\nkh: " << kh
2590  << "\nrho: " << 1.0/(1.0-(1.0/3.0)*pow(kh,2.0))
2591  << "\nomega: " << omega;
2592 
2593  // If we're working on the coarsest grids or kh is high enough
2594  if (kh>2*cos(pi*h/2))
2595  {
2596  oomph_info << "\n(0,omega_2): (0," << omega_2 << ")";
2597  }
2598 
2599  // Carry on outputting
2600  oomph_info << "\n--------------------------"
2601  << std::endl;
2602  }
2603  //===================================================TO_BE_DELETED=========
2604 
2605 
2606  // Loop over the levels and assign the pre- and post-smoother pointers
2607  for (unsigned i=0;i<Nlevel-1;i++)
2608  {
2609  // If the pre-smoother factory function pointer hasn't been assigned
2610  // then we simply create a new instance of the ComplexDampedJacobi
2611  // smoother which is the default pre-smoother
2612  if (0==Pre_smoother_factory_function_pt)
2613  {
2614  // If the value of kh is strictly less than 0.5 then use damped Jacobi
2615  // as a smoother on this level otherwise use complex GMRES as a smoother
2616  // [see Elman et al."A multigrid method enhanced by Krylov subspace
2617  // iteration for discrete Helmholtz equations", p.1305]
2618  if (Wavenumber*Max_edge_width[i]<0.5)
2619  {
2620  // Make a new ComplexDampedJacobi object and assign its pointer
2621  ComplexDampedJacobi<CRDoubleMatrix>* damped_jacobi_pt=
2623 
2624  // Assign the pre-smoother pointer
2625  Pre_smoothers_storage_pt[i]=damped_jacobi_pt;
2626 
2627  // Make the smoother calculate the value of omega and store it
2628  damped_jacobi_pt->calculate_omega(Wavenumber,Max_edge_width[i]);
2629  }
2630  else
2631  {
2632  // Make a new ComplexGMRES object and assign its pointer
2634 
2635  // Assign the pre-smoother pointer
2636  Pre_smoothers_storage_pt[i]=gmres_pt;
2637 
2638  // Enable the documentation of the convergence (HIERHER: Temporary!)
2639  //gmres_pt->enable_doc_convergence_history();
2640  }
2641  }
2642  // Otherwise we use the pre-smoother factory function pointer to
2643  // generate a new pre-smoother
2644  else
2645  {
2646  // Get a pointer to an object of the same type as the pre-smoother
2647  Pre_smoothers_storage_pt[i]=
2648  (*Pre_smoother_factory_function_pt)();
2649  } // if (0==Pre_smoother_factory_function_pt)
2650 
2651  // If the post-smoother factory function pointer hasn't been assigned
2652  // then we simply create a new instance of the ComplexDampedJacobi smoother
2653  // which is the default post-smoother
2654  if (0==Post_smoother_factory_function_pt)
2655  {
2656  // If the value of kh is strictly less than 0.52 then use damped Jacobi
2657  // as a smoother on this level otherwise use complex GMRES as a smoother
2658  // [see Elman et al."A multigrid method enhanced by Krylov subspace
2659  // iteration for discrete Helmholtz equations", p.1295]
2660  if (Wavenumber*Max_edge_width[i]<0.5)
2661  {
2662  // Make a new ComplexDampedJacobi object and assign its pointer
2663  ComplexDampedJacobi<CRDoubleMatrix>* damped_jacobi_pt=
2665 
2666  // Assign the pre-smoother pointer
2667  Post_smoothers_storage_pt[i]=damped_jacobi_pt;
2668 
2669  // Make the smoother calculate the value of omega and store it
2670  damped_jacobi_pt->calculate_omega(Wavenumber,Max_edge_width[i]);
2671  }
2672  else
2673  {
2674  // Make a new ComplexGMRES object and assign its pointer
2676 
2677  // Assign the pre-smoother pointer
2678  Post_smoothers_storage_pt[i]=gmres_pt;
2679  }
2680  }
2681  // Otherwise we use the post-smoother factory function pointer to
2682  // generate a new post-smoother
2683  else
2684  {
2685  // Get a pointer to an object of the same type as the post-smoother
2686  Post_smoothers_storage_pt[i]=
2687  (*Post_smoother_factory_function_pt)();
2688  }
2689  } // if (0==Post_smoother_factory_function_pt)
2690 
2691  // Set the tolerance for the pre- and post-smoothers. The norm of the
2692  // solution which is compared against the tolerance is calculated
2693  // differently to the multigrid solver. To ensure that the smoother
2694  // continues to solve for the prescribed number of iterations we
2695  // lower the tolerance
2696  for (unsigned i=0;i<Nlevel-1;i++)
2697  {
2698  // Set the tolerance of the i-th level pre-smoother
2699  Pre_smoothers_storage_pt[i]->tolerance()=1.0e-16;
2700 
2701  // Set the tolerance of the i-th level post-smoother
2702  Post_smoothers_storage_pt[i]->tolerance()=1.0e-16;
2703  }
2704 
2705  // Set the number of pre- and post-smoothing iterations in each
2706  // pre- and post-smoother
2707  for (unsigned i=0;i<Nlevel-1;i++)
2708  {
2709  // Set the number of pre-smoothing iterations as the value of Npre_smooth
2710  Pre_smoothers_storage_pt[i]->max_iter()=Npre_smooth;
2711 
2712  // HIERHER: delete
2713  if (i==1)
2714  {
2715  Pre_smoothers_storage_pt[i]->max_iter()=2;
2716  }
2717 
2718  // Set the number of pre-smoothing iterations as the value of Npost_smooth
2719  Post_smoothers_storage_pt[i]->max_iter()=Npost_smooth;
2720  }
2721 
2722  // Complete the setup of all of the smoothers
2723  for (unsigned i=0;i<Nlevel-1;i++)
2724  {
2725  // Pass a pointer to the system matrix on the i-th level to the i-th
2726  // level pre-smoother
2727  Pre_smoothers_storage_pt[i]->
2728  complex_smoother_setup(Mg_matrices_storage_pt[i]);
2729 
2730  // Pass a pointer to the system matrix on the i-th level to the i-th
2731  // level post-smoother
2732  Post_smoothers_storage_pt[i]->
2733  complex_smoother_setup(Mg_matrices_storage_pt[i]);
2734  }
2735 
2736  // Set up the distributions of each smoother
2737  for (unsigned i=0;i<Nlevel-1;i++)
2738  {
2739  // Get the number of dofs on the i-th level of the hierarchy.
2740  // This will be the same as the size of the solution vector
2741  // associated with the i-th level
2742  unsigned n_dof=X_mg_vectors_storage[i][0].nrow();
2743 
2744  // Create a LinearAlgebraDistribution which will be passed to the
2745  // linear solver
2746  LinearAlgebraDistribution dist(Mg_hierarchy_pt[i]->communicator_pt(),
2747  n_dof,false);
2748 
2749 #ifdef PARANOID
2750 #ifdef OOMPH_HAS_MPI
2751  // Set up the warning messages
2752  std::string warning_message="Setup of pre- and post-smoother distribution ";
2753  warning_message+="has not been tested with MPI.";
2754 
2755  // If we're not running the code in serial
2756  if (dist.communicator_pt()->nproc()>1)
2757  {
2758  // Throw a warning
2759  OomphLibWarning(warning_message,
2760  OOMPH_CURRENT_FUNCTION,
2761  OOMPH_EXCEPTION_LOCATION);
2762  }
2763 #endif
2764 #endif
2765 
2766  // Build the distribution of the pre-smoother
2767  Pre_smoothers_storage_pt[i]->build_distribution(dist);
2768 
2769  // Build the distribution of the post-smoother
2770  Post_smoothers_storage_pt[i]->build_distribution(dist);
2771  }
2772 
2773  // Disable the smoother output
2774  if (!Doc_time)
2775  {
2776  // Disable time documentation from the smoothers and the SuperLU linear
2777  // solver (this latter is only done on the coarsest level where it is
2778  // required)
2779  disable_smoother_and_superlu_doc_time();
2780  }
2781 
2782  // If we're allowed to output
2783  if (!Suppress_all_output)
2784  {
2785  // Stop clock
2786  double t_m_end=TimingHelpers::timer();
2787  double total_setup_time=double(t_m_end-t_m_start);
2788  oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
2789  << total_setup_time << std::endl;
2790 
2791  // Notify user that the extraction is complete
2792  oomph_info << "\n=================="
2793  << "Smoother Setup Complete"
2794  << "=================" << std::endl;
2795  }
2796  } // End of setup_smoothers
2797 
2798 
2799 //===================================================================
2800 /// \short Setup the interpolation matrices
2801 //===================================================================
2802  template<unsigned DIM>
2804  {
2805  // Variable to hold the number of sons in an element
2806  unsigned n_sons;
2807 
2808  // Number of son elements
2809  if (DIM==2)
2810  {
2811  n_sons=4;
2812  }
2813  else if (DIM==3)
2814  {
2815  n_sons=8;
2816  }
2817  else
2818  {
2819  std::ostringstream error_message_stream;
2820  error_message_stream << "DIM should be 2 or 3 not "
2821  << DIM << std::endl;
2822  throw OomphLibError(error_message_stream.str(),
2823  OOMPH_CURRENT_FUNCTION,
2824  OOMPH_EXCEPTION_LOCATION);
2825  }
2826 
2827 #ifdef PARANOID
2828  // PARANOID check - for our implementation we demand that the first
2829  // submesh is the refineable bulk mesh that is provided by the function
2830  // mg_bulk_mesh_pt. Since each mesh in the hierarchy will be set up
2831  // in the same manner through the problem constructor we only needed
2832  // to check the finest level mesh
2833  if (Mg_hierarchy_pt[0]->mesh_pt(0)!=Mg_hierarchy_pt[0]->mg_bulk_mesh_pt())
2834  {
2835  // Create an output stream
2836  std::ostringstream error_message_stream;
2837 
2838  // Create the error message
2839  error_message_stream << "MG solver requires the first submesh be the "
2840  << "refineable bulk mesh provided by the user."
2841  << std::endl;
2842 
2843  // Throw the error message
2844  throw OomphLibError(error_message_stream.str(),
2845  OOMPH_CURRENT_FUNCTION,
2846  OOMPH_EXCEPTION_LOCATION);
2847  }
2848 #endif
2849 
2850  // Vector of local coordinates in the element
2851  Vector<double> s(DIM,0.0);
2852 
2853  // Vector to contain the (Eulerian) spatial location of the fine node.
2854  // This will only be used if we have a PML layer attached to the mesh
2855  // which will require the use of locate_zeta functionality
2856  Vector<double> fine_node_position(DIM,0.0);
2857 
2858  // Find the number of nodes in an element in the mesh. Since this
2859  // quantity will be the same in all meshes we can precompute it here
2860  // using the original problem pointer
2861  unsigned nnod_element=dynamic_cast<FiniteElement*>
2862  (Mg_problem_pt->mesh_pt()->element_pt(0))->nnode();
2863 
2864  // Initialise a null pointer which will be used to point to an object
2865  // of the class MeshAsGeomObject
2866  MeshAsGeomObject* coarse_mesh_from_obj_pt=0;
2867 
2868  // Loop over each level (apart from the coarsest level; an interpolation
2869  // matrix and thus a restriction matrix is not needed here), with 0 being
2870  // the finest level and Nlevel-1 being the coarsest level
2871  for (unsigned level=0;level<Nlevel-1;level++)
2872  {
2873  // Store the ID of the fine level
2874  unsigned fine_level=level;
2875 
2876  // Store the ID of the coarse level
2877  unsigned coarse_level=level+1;
2878 
2879  // Make a pointer to the mesh on the finer level and dynamic_cast
2880  // it as an object of the refineable mesh class
2881  Mesh* ref_fine_mesh_pt=Mg_hierarchy_pt[fine_level]->mesh_pt();
2882 
2883  // Make a pointer to the mesh on the coarse level and dynamic_cast
2884  // it as an object of the refineable mesh class
2885  Mesh* ref_coarse_mesh_pt=Mg_hierarchy_pt[coarse_level]->mesh_pt();
2886 
2887  // Access information about the number of elements in the fine mesh
2888  // from the pointer to the fine mesh (to loop over the rows of the
2889  // interpolation matrix)
2890  unsigned fine_n_element=ref_fine_mesh_pt->nelement();
2891 
2892  // Find the number of elements in the bulk mesh
2893  unsigned n_bulk_mesh_element=
2894  Mg_hierarchy_pt[fine_level]->mg_bulk_mesh_pt()->nelement();
2895 
2896  // If there are elements apart from those in the bulk mesh
2897  // then we require locate_zeta functionality. For this reason
2898  // we have to assign a value to the MeshAsGeomObject pointer
2899  // which we do by creating a MeshAsGeomObject from the coarse
2900  // mesh pointer
2901  if (fine_n_element>n_bulk_mesh_element)
2902  {
2903  // Assign the pointer to coarse_mesh_from_obj_pt
2904  coarse_mesh_from_obj_pt=
2905  new MeshAsGeomObject(Mg_hierarchy_pt[coarse_level]->mesh_pt());
2906  }
2907 
2908  // Find the number of dofs in the fine mesh
2909  unsigned n_fine_dof=Mg_hierarchy_pt[fine_level]->ndof();
2910 
2911  // Find the number of dofs in the coarse mesh
2912  unsigned n_coarse_dof=Mg_hierarchy_pt[coarse_level]->ndof();
2913 
2914  // To get the number of rows in the interpolation matrix we divide
2915  // the number of dofs on each level by 2 since there are 2 dofs at
2916  // each node in the mesh corresponding to the real and imaginary
2917  // part of the solution:
2918 
2919  // Get the numbers of rows in the interpolation matrix.
2920  unsigned n_row=n_fine_dof/2.0;
2921 
2922  // Get the numbers of columns in the interpolation matrix
2923  unsigned n_col=n_coarse_dof/2.0;
2924 
2925  // Mapping relating the pointers to related elements in the coarse
2926  // and fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]. If
2927  // the element in the fine mesh has been unrefined between these two
2928  // levels, this map returns the father element in the coarsened mesh.
2929  // If this element in the fine mesh has not been unrefined, the map
2930  // returns the pointer to the same-sized equivalent element in the
2931  // coarsened mesh.
2932  std::map<RefineableQElement<DIM>*,
2933  RefineableQElement<DIM>*> coarse_mesh_reference_element_pt;
2934 
2935  // Variable to count the number of elements in the fine mesh
2936  unsigned e_fine=0;
2937 
2938  // Variable to count the number of elements in the coarse mesh
2939  unsigned e_coarse=0;
2940 
2941  // While loop over fine elements (while because we're incrementing the
2942  // counter internally if the element was unrefined...). We only bother
2943  // going until we've covered all of the elements in the bulk mesh
2944  // because once we reach the PML layer (if there is one) there will be
2945  // no tree structure to match in between levels as PML elements are
2946  // simply regenerated once the bulk mesh has been refined.
2947  while(e_fine<n_bulk_mesh_element)
2948  {
2949  // Pointer to element in fine mesh
2950  RefineableQElement<DIM>* el_fine_pt=
2951  dynamic_cast<RefineableQElement<DIM>*>
2952  (ref_fine_mesh_pt->finite_element_pt(e_fine));
2953 
2954  // Pointer to element in coarse mesh
2955  RefineableQElement<DIM>* el_coarse_pt=
2956  dynamic_cast<RefineableQElement<DIM>*>
2957  (ref_coarse_mesh_pt->finite_element_pt(e_coarse));
2958 
2959  // If the levels are different then the element in the fine
2960  // mesh has been unrefined between these two levels
2961  if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level())
2962  {
2963  // The element in the fine mesh has been unrefined between these two
2964  // levels. Hence it and its three brothers (ASSUMED to be stored
2965  // consecutively in the Mesh's vector of pointers to its constituent
2966  // elements -- we'll check this!) share the same father element in
2967  // the coarse mesh, currently pointed to by el_coarse_pt.
2968  for(unsigned i=0;i<n_sons;i++)
2969  {
2970  // Set mapping to father element in coarse mesh
2971  coarse_mesh_reference_element_pt[
2972  dynamic_cast<RefineableQElement<DIM>*>
2973  (ref_fine_mesh_pt->finite_element_pt(e_fine))]=el_coarse_pt;
2974 
2975  // Increment counter for elements in fine mesh
2976  e_fine++;
2977  }
2978  }
2979  // The element in the fine mesh has not been unrefined between
2980  // these two levels, so the reference element is the same-sized
2981  // equivalent element in the coarse mesh
2982  else
2983  {
2984  // Set the mapping between the two elements since they are
2985  // the same element
2986  coarse_mesh_reference_element_pt[el_fine_pt]=el_coarse_pt;
2987 
2988  // Increment counter for elements in fine mesh
2989  e_fine++;
2990  } // if (el_fine_pt->tree_pt()->level()!=...)
2991 
2992  // Increment counter for elements in coarse mesh
2993  e_coarse++;
2994  } // while(e_fine<n_bulk_mesh_element)
2995 
2996  // To allow an update of a row only once we use STL vectors for bools
2997  std::vector<bool> contribution_made(n_row,false);
2998 
2999  // Create the row start vector whose i-th entry will contain the index
3000  // in column_start where the entries in the i-th row of the matrix start
3001  Vector<int> row_start(n_row+1);
3002 
3003  // Create the column index vector whose entries will store the column
3004  // index of each contribution, i.e. the global equation of the coarse
3005  // mesh node which supplies a contribution. We don't know how many
3006  // entries this will have so we dynamically allocate entries at run time
3007  // hierher: should we reserve this before run time? We can probably estimate
3008  // the number of contributions from the value of nnode_1d
3009  Vector<int> column_index;
3010 
3011  // Create the value vector which will be used to store the nonzero
3012  // entries in the CRDoubleMatrix. We don't know how many entries this
3013  // will have either so we dynamically allocate its entries at run time
3014  Vector<double> value;
3015 
3016  // The value of index will tell us which row of the interpolation matrix
3017  // we're working on in the following for loop
3018  // hierher: should this worry us? We assume that every node we cover is
3019  // the next node in the mesh (we loop over the elements and the nodes
3020  // inside that). This does work but it may not work for some meshes
3021  unsigned index=0;
3022 
3023  // New loop to go over each element in the fine mesh
3024  for (unsigned k=0;k<fine_n_element;k++)
3025  {
3026  // Set a pointer to the element in the fine mesh
3027  RefineableQElement<DIM>* el_fine_pt=
3028  dynamic_cast<RefineableQElement<DIM>*>
3029  (ref_fine_mesh_pt->finite_element_pt(k));
3030 
3031  // Initialise a null pointer which will point to the corresponding
3032  // coarse mesh element. If we're in the bulk mesh, this will point
3033  // to the parent element of the fine mesh element or itself (if
3034  // there was no refinement between the levels). If, however, we're
3035  // in the PML layer then this will point to element that encloses
3036  // the fine mesh PML element
3037  RefineableQElement<DIM>* el_coarse_pt=0;
3038 
3039  // Variable to store the son type of the fine mesh element with
3040  // respect to the corresponding coarse mesh element (set to OMEGA
3041  // if no unrefinement took place)
3042  int son_type=0;
3043 
3044  // If we're in the bulk mesh we can assign the coarse mesh element
3045  // pointer right now using the map coarse_mesh_reference_element_pt
3046  if (k<n_bulk_mesh_element)
3047  {
3048  // Get the reference element (either the father element or the
3049  // same-sized element) in the coarse mesh
3050  el_coarse_pt=coarse_mesh_reference_element_pt[el_fine_pt];
3051 
3052  // Check if the elements are different on both levels (i.e. to check
3053  // if any unrefinement took place)
3054  if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level())
3055  {
3056  // If the element was refined then we need the son type
3057  son_type=el_fine_pt->tree_pt()->son_type();
3058  }
3059  else
3060  {
3061  // If the element is the same in both meshes
3062  son_type=Tree::OMEGA;
3063  }
3064  } // if (k<n_bulk_mesh_element)
3065 
3066  // Loop through all the nodes in an element in the fine mesh
3067  for(unsigned i=0;i<nnod_element;i++)
3068  {
3069  // Row number in interpolation matrix: Global equation number
3070  // of the d.o.f. stored at this node in the fine element
3071  int ii=el_fine_pt->node_pt(i)->eqn_number(0);
3072 
3073  // Check whether or not the node is a proper d.o.f.
3074  if (ii>=0)
3075  {
3076  // Only assign values to the given row of the interpolation matrix
3077  // if they haven't already been assigned. Notice, we check if the
3078  // (ii/2)-th entry has been set. This is because there are two dofs
3079  // at each node so dividing by two will give us the row number
3080  // associated with this node
3081  if(contribution_made[ii/2]==false)
3082  {
3083  // The value of index was initialised when we allocated space
3084  // for the three vectors to store the CRDoubleMatrix. We use
3085  // index to go through the entries of the row_start vector.
3086  // The next row starts at the value.size() (draw out the entries
3087  // in value if this doesn't make sense noting that the storage
3088  // for the vector 'value' is dynamically allocated)
3089  row_start[index]=value.size();
3090 
3091  // If we're in the PML layer then we need to find the parent element
3092  // associated with a given node using locate_zeta. Unfortunately,
3093  // if this is the case and a given local node in the fine mesh
3094  // element lies on the interface between two elements (E1) and (E2)
3095  // in the coarse mesh then either element will be returned. Suppose,
3096  // for instance, a pointer to (E1) is returned. If the next local
3097  // node lies in the (E2) then the contributions which we pick up
3098  // from the local nodes in (E1) will most likely be incorrect while
3099  // the column index (the global equation number of the coarse mesh
3100  // node) corresponding to the given contribution will definitely be
3101  // incorrect. If we're not using linear quad elements we can get
3102  // around this by using locate_zeta carefully by identifying a node
3103  // which is internal to the fine mesh element. This will allow us to
3104  // correctly obtain the corresponding element in the coarse mesh
3105  if (k>=n_bulk_mesh_element)
3106  {
3107  // Get the (Eulerian) spatial location of the i-th local node
3108  // in the fine mesh element
3109  el_fine_pt->node_pt(i)->position(fine_node_position);
3110 
3111  // Initalise a null pointer to point to an object of the class
3112  // GeomObject
3113  GeomObject* el_pt=0;
3114 
3115  // Get the reference element (either the father element or the
3116  // same-sized element) in the coarse-mesh PML layer using
3117  // the locate_zeta functionality
3118  coarse_mesh_from_obj_pt->locate_zeta(fine_node_position,el_pt,s);
3119 
3120  // Upcast the GeomElement object to a RefineableQElement object
3121  el_coarse_pt=dynamic_cast<RefineableQElement<DIM>*>(el_pt);
3122  }
3123  else
3124  {
3125  // Calculate the local coordinates of the given node
3126  el_fine_pt->local_coordinate_of_node(i,s);
3127 
3128  // Find the local coordinates s, of the present node, in the
3129  // reference element in the coarse mesh, given the element's
3130  // son_type (e.g. SW,SE... )
3131  level_up_local_coord_of_node(son_type,s);
3132  }
3133 
3134  // Allocate space for shape functions in the coarse mesh
3135  Shape psi(el_coarse_pt->nnode());
3136 
3137  // Set the shape function in the reference element
3138  el_coarse_pt->shape(s,psi);
3139 
3140  // Auxiliary storage
3141  std::map<unsigned,double> contribution;
3142  Vector<unsigned> keys;
3143 
3144  // Find the number of nodes in an element in the coarse mesh
3145  unsigned nnod_coarse=el_coarse_pt->nnode();
3146 
3147  // Loop through all the nodes in the reference element
3148  for (unsigned j=0;j<nnod_coarse;j++)
3149  {
3150  // Column number in interpolation matrix: Global equation
3151  // number of the d.o.f. stored at this node in the coarse
3152  // element
3153  int jj=el_coarse_pt->node_pt(j)->eqn_number(0);
3154 
3155  // If the value stored at this node is pinned or hanging
3156  if (jj<0)
3157  {
3158  // Hanging node: In this case we need to accumulate the
3159  // contributions from the master nodes
3160  if (el_coarse_pt->node_pt(j)->is_hanging())
3161  {
3162  // Find the number of master nodes of the hanging
3163  // the node in the reference element
3164  HangInfo* hang_info_pt=el_coarse_pt->node_pt(j)->hanging_pt();
3165  unsigned nmaster=hang_info_pt->nmaster();
3166 
3167  // Loop over the master nodes
3168  for (unsigned i_master=0;i_master<nmaster;i_master++)
3169  {
3170  // Set up a pointer to the master node
3171  Node* master_node_pt=hang_info_pt->master_node_pt(i_master);
3172 
3173  // The column number in the interpolation matrix: the
3174  // global equation number of the d.o.f. stored at this master
3175  // node for the coarse element
3176  int master_jj=master_node_pt->eqn_number(0);
3177 
3178  // Is the master node a proper d.o.f.?
3179  if (master_jj>=0)
3180  {
3181  // If the weight of the master node is non-zero
3182  if (psi(j)*hang_info_pt->master_weight(i_master)!=0.0)
3183  {
3184  contribution[master_jj/2]+=
3185  psi(j)*hang_info_pt->master_weight(i_master);
3186  }
3187  } // if (master_jj>=0)
3188  } // for (unsigned i_master=0;i_master<nmaster;i_master++)
3189  } // if (el_coarse_pt->node_pt(j)->is_hanging())
3190  }
3191  // In the case that the node is not pinned or hanging
3192  else
3193  {
3194  // If we can get a nonzero contribution from the shape function
3195  if (psi(j)!=0.0)
3196  {
3197  contribution[jj/2]+=psi(j);
3198  }
3199  } // if (jj<0) else
3200  } // for (unsigned j=0;j<nnod_coarse;j++)
3201 
3202  // Put the contributions into the value vector
3203  for (std::map<unsigned,double>::iterator it=contribution.begin();
3204  it!=contribution.end(); ++it)
3205  {
3206  if (it->second!=0)
3207  {
3208  // The first entry in contribution is the column index
3209  column_index.push_back(it->first);
3210 
3211  // The second entry in contribution is the value of the
3212  // contribution
3213  value.push_back(it->second);
3214  }
3215  } // for (std::map<unsigned,double>::iterator it=...)
3216 
3217  // Increment the value of index by 1 to indicate that we have
3218  // finished the row, or equivalently, covered another global
3219  // node in the fine mesh
3220  index++;
3221 
3222  // Change the entry in contribution_made to true now to indicate
3223  // that the row has been filled
3224  contribution_made[ii/2]=true;
3225  } // if(contribution_made[ii/2]==false)
3226  } // if (ii>=0)
3227  } // for(unsigned i=0;i<nnod_element;i++)
3228  } // for (unsigned k=0;k<fine_n_element;k++)
3229 
3230  // Set the last entry in the row_start vector
3231  row_start[n_row]=value.size();
3232 
3233  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
3234  // using the vectors value, row_start, column_index and the value
3235  // of fine_n_unknowns and coarse_n_unknowns
3236  interpolation_matrix_set(level,
3237  value,
3238  column_index,
3239  row_start,
3240  n_col,
3241  n_row);
3242  } // for (unsigned level=0;level<Nlevel-1;level++)
3243  } // End of setup_interpolation_matrices
3244 
3245 //===================================================================
3246 /// \short Setup the interpolation matrices
3247 //===================================================================
3248  template<unsigned DIM>
3251  {
3252  // Variable to hold the number of sons in an element
3253  unsigned n_sons;
3254 
3255  // Number of son elements
3256  if (DIM==2)
3257  {
3258  n_sons=4;
3259  }
3260  else if (DIM==3)
3261  {
3262  n_sons=8;
3263  }
3264  else
3265  {
3266  std::ostringstream error_message_stream;
3267  error_message_stream << "DIM should be 2 or 3 not "
3268  << DIM << std::endl;
3269  throw OomphLibError(error_message_stream.str(),
3270  OOMPH_CURRENT_FUNCTION,
3271  OOMPH_EXCEPTION_LOCATION);
3272  }
3273 
3274  // Vector of local coordinates in the element
3275  Vector<double> s(DIM,0.0);
3276 
3277  // Loop over each level (apart from the coarsest level; an interpolation
3278  // matrix and thus a restriction matrix is not needed here), with 0 being
3279  // the finest level and Nlevel-1 being the coarsest level
3280  for (unsigned level=0;level<Nlevel-1;level++)
3281  {
3282  // Assign values to a couple of variables to help out
3283  unsigned coarse_level=level+1;
3284  unsigned fine_level=level;
3285 
3286  // Make a pointer to the mesh on the finer level and dynamic_cast
3287  // it as an object of the refineable mesh class
3288  Mesh* ref_fine_mesh_pt=Mg_hierarchy_pt[fine_level]->mg_bulk_mesh_pt();
3289 
3290  // To use the locate zeta functionality the coarse mesh must be of the
3291  // type MeshAsGeomObject
3292  MeshAsGeomObject* coarse_mesh_from_obj_pt=
3293  new MeshAsGeomObject(Mg_hierarchy_pt[coarse_level]->mg_bulk_mesh_pt());
3294 
3295  // Access information about the number of degrees of freedom
3296  // from the pointers to the problem on each level
3297  unsigned coarse_n_unknowns=Mg_hierarchy_pt[coarse_level]->ndof();
3298  unsigned fine_n_unknowns=Mg_hierarchy_pt[fine_level]->ndof();
3299 
3300  // Make storage vectors to form the interpolation matrix using a
3301  // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
3302  // row_start contains entries which tells us at which entry of the
3303  // vector column_start we start on the i-th row of the actual matrix.
3304  // The entries of column_start indicate the column position of each
3305  // non-zero entry in every row. This runs through the first row
3306  // from left to right then the second row (again, left to right)
3307  // and so on until the end. The entries in value are the entries in
3308  // the interpolation matrix. The vector column_start has the same length
3309  // as the vector value.
3310  Vector<double> value;
3311  Vector<int> column_index;
3312  Vector<int> row_start(fine_n_unknowns+1);
3313 
3314  // Vector to contain the (Eulerian) spatial location of the fine node
3315  Vector<double> fine_node_position(DIM);
3316 
3317  // Find the number of nodes in the mesh
3318  unsigned n_node_fine_mesh=ref_fine_mesh_pt->nnode();
3319 
3320  // Loop over the unknowns in the mesh
3321  for (unsigned i_fine_node=0;i_fine_node<n_node_fine_mesh;i_fine_node++)
3322  {
3323  // Set a pointer to the i_fine_unknown-th node in the fine mesh
3324  Node* fine_node_pt=ref_fine_mesh_pt->node_pt(i_fine_node);
3325 
3326  // Get the global equation number
3327  int i_fine=fine_node_pt->eqn_number(0);
3328 
3329  // If the node is a proper d.o.f.
3330  if (i_fine>=0)
3331  {
3332  // Row number in interpolation matrix: Global equation number
3333  // of the d.o.f. stored at this node in the fine element
3334  row_start[i_fine]=value.size();
3335 
3336  // Get the (Eulerian) spatial location of the fine node
3337  fine_node_pt->position(fine_node_position);
3338 
3339  // Create a null pointer to the GeomObject class
3340  GeomObject* el_pt=0;
3341 
3342  // Get the reference element (either the father element or the
3343  // same-sized element) in the coarse mesh using locate_zeta
3344  coarse_mesh_from_obj_pt->locate_zeta(fine_node_position,el_pt,s);
3345 
3346  // Upcast GeomElement as a FiniteElement
3347  FiniteElement* el_coarse_pt=dynamic_cast<FiniteElement*>(el_pt);
3348 
3349  // Find the number of nodes in the element
3350  unsigned n_node=el_coarse_pt->nnode();
3351 
3352  // Allocate space for shape functions in the coarse mesh
3353  Shape psi(n_node);
3354 
3355  // Calculate the geometric shape functions at local coordinate s
3356  el_coarse_pt->shape(s,psi);
3357 
3358  // Auxiliary storage
3359  std::map<unsigned,double> contribution;
3360  Vector<unsigned> keys;
3361 
3362  // Loop through all the nodes in the (coarse mesh) element containing the
3363  // node pointed to by fine_node_pt (fine mesh)
3364  for(unsigned j_node=0;j_node<n_node;j_node++)
3365  {
3366  // Get the j_coarse_unknown-th node in the coarse element
3367  Node* coarse_node_pt=el_coarse_pt->node_pt(j_node);
3368 
3369  // Column number in interpolation matrix: Global equation number of
3370  // the d.o.f. stored at this node in the coarse element
3371  int j_coarse=coarse_node_pt->eqn_number(0);
3372 
3373  // If the value stored at this node is pinned or hanging
3374  if (j_coarse<0)
3375  {
3376  // Hanging node: In this case we need to accumulate the
3377  // contributions from the master nodes
3378  if (el_coarse_pt->node_pt(j_node)->is_hanging())
3379  {
3380  // Find the number of master nodes of the hanging
3381  // the node in the reference element
3382  HangInfo* hang_info_pt=coarse_node_pt->hanging_pt();
3383  unsigned nmaster=hang_info_pt->nmaster();
3384 
3385  // Loop over the master nodes
3386  for (unsigned i_master=0;i_master<nmaster;i_master++)
3387  {
3388  // Set up a pointer to the master node
3389  Node* master_node_pt=hang_info_pt->master_node_pt(i_master);
3390 
3391  // The column number in the interpolation matrix: the
3392  // global equation number of the d.o.f. stored at this master
3393  // node for the coarse element
3394  int master_jj=master_node_pt->eqn_number(0);
3395 
3396  // Is the master node a proper d.o.f.?
3397  if (master_jj>=0)
3398  {
3399  // If the weight of the master node is non-zero
3400  if (psi(j_node)*hang_info_pt->master_weight(i_master)!=0.0)
3401  {
3402  contribution[master_jj]+=
3403  psi(j_node)*hang_info_pt->master_weight(i_master);
3404  }
3405  } // End of if statement (check it's not a boundary node)
3406  } // End of the loop over the master nodes
3407  } // End of the if statement for only hanging nodes
3408  } // End of the if statement for pinned or hanging nodes
3409  // In the case that the node is not pinned or hanging
3410  else
3411  {
3412  // If we can get a nonzero contribution from the shape function
3413  // at the j_node-th node in the element
3414  if (psi(j_node)!=0.0)
3415  {
3416  contribution[j_coarse]+=psi(j_node);
3417  }
3418  } // End of the if-else statement (check if the node was pinned/hanging)
3419  } // Finished loop over the nodes j in the reference element (coarse)
3420 
3421  // Put the contributions into the value vector
3422  for (std::map<unsigned,double>::iterator it=contribution.begin();
3423  it!=contribution.end(); ++it)
3424  {
3425  if (it->second!=0)
3426  {
3427  value.push_back(it->second);
3428  column_index.push_back(it->first);
3429  }
3430  } // End of putting contributions into the value vector
3431  } // End check (whether or not the fine node was a d.o.f.)
3432  } // End of the for-loop over nodes in the fine mesh
3433 
3434  // Set the last entry of row_start
3435  row_start[fine_n_unknowns]=value.size();
3436 
3437  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
3438  // using the vectors value, row_start, column_index and the value
3439  // of fine_n_unknowns and coarse_n_unknowns
3440  interpolation_matrix_set(level,
3441  value,
3442  column_index,
3443  row_start,
3444  coarse_n_unknowns,
3445  fine_n_unknowns);
3446  } // End of loop over each level
3447  } // End of setup_interpolation_matrices_unstructured
3448 
3449 //=========================================================================
3450 /// \short Given the son type of the element and the local coordinate s of
3451 /// a given node in the son element, return the local coordinate s in its
3452 /// father element. 3D case.
3453 //=========================================================================
3454  template<>
3457  {
3458  // If the element is unrefined between the levels the local coordinate
3459  // of the node in one element is the same as that in the other element
3460  // therefore we only need to perform calculations if the levels are
3461  // different (i.e. son_type is not OMEGA)
3462  if (son_type!=Tree::OMEGA)
3463  {
3464  // Scale the local coordinate from the range [-1,1]x[-1,1]x[-1,1]
3465  // to the range [0,1]x[0,1]x[0,1] to match the width of the local
3466  // coordinate range of the fine element from the perspective of
3467  // the father element. This then simply requires a shift of the
3468  // coordinates to match which type of son element we're dealing with
3469  s[0]=(s[0]+1.0)/2.0;
3470  s[1]=(s[1]+1.0)/2.0;
3471  s[2]=(s[2]+1.0)/2.0;
3472 
3473  // Cases: The son_type determines how the local coordinates should be
3474  // shifted to give the local coordinates in the coarse mesh element
3475  switch(son_type)
3476  {
3477  case OcTreeNames::LDF:
3478  s[0]-=1;
3479  s[1]-=1;
3480  break;
3481 
3482  case OcTreeNames::LDB:
3483  s[0]-=1;
3484  s[1]-=1;
3485  s[2]-=1;
3486  break;
3487 
3488  case OcTreeNames::LUF:
3489  s[0]-=1;
3490  break;
3491 
3492  case OcTreeNames::LUB:
3493  s[0]-=1;
3494  s[2]-=1;
3495  break;
3496 
3497  case OcTreeNames::RDF:
3498  s[1]-=1;
3499  break;
3500 
3501  case OcTreeNames::RDB:
3502  s[1]-=1;
3503  s[2]-=1;
3504  break;
3505 
3506  case OcTreeNames::RUF:
3507  break;
3508 
3509  case OcTreeNames::RUB:
3510  s[2]-=1;
3511  break;
3512  }
3513  } // if son_type!=Tree::OMEGA
3514  } // End of level_up_local_coord_of_node
3515 
3516 //=========================================================================
3517 /// \short Given the son type of the element and the local coordinate s of
3518 /// a given node in the son element, return the local coordinate s in its
3519 /// father element. 2D case.
3520 //=========================================================================
3521  template<>
3524  {
3525  // If the element is unrefined between the levels the local coordinate
3526  // of the node in one element is the same as that in the other element
3527  // therefore we only need to perform calculations if the levels are
3528  // different (i.e. son_type is not OMEGA)
3529  if (son_type!=Tree::OMEGA)
3530  {
3531  // Scale the local coordinate from the range [-1,1]x[-1,1] to the range
3532  // [0,1]x[0,1] to match the width of the local coordinate range of the
3533  // fine element from the perspective of the father element. This
3534  // then simply requires a shift of the coordinates to match which type
3535  // of son element we're dealing with
3536  s[0]=(s[0]+1.0)/2.0;
3537  s[1]=(s[1]+1.0)/2.0;
3538 
3539  // Cases: The son_type determines how the local coordinates should be
3540  // shifted to give the local coordinates in the coarse mesh element
3541  switch(son_type)
3542  {
3543  // If we're dealing with the bottom-left element we need to shift
3544  // the range [0,1]x[0,1] to [-1,0]x[-1,0]
3545  case QuadTreeNames::SW:
3546  s[0]-=1;
3547  s[1]-=1;
3548  break;
3549 
3550  // If we're dealing with the bottom-right element we need to shift
3551  // the range [0,1]x[0,1] to [0,1]x[-1,0]
3552  case QuadTreeNames::SE:
3553  s[1]-=1;
3554  break;
3555 
3556  // If we're dealing with the top-right element we need to shift the
3557  // range [0,1]x[0,1] to [0,1]x[0,1], i.e. nothing needs to be done
3558  case QuadTreeNames::NE:
3559  break;
3560 
3561  // If we're dealing with the top-left element we need to shift
3562  // the range [0,1]x[0,1] to [-1,0]x[0,1]
3563  case QuadTreeNames::NW:
3564  s[0]-=1;
3565  break;
3566  }
3567  } // if son_type!=Tree::OMEGA
3568  } // End of level_up_local_coord_of_node
3569 
3570 //===================================================================
3571 /// \short Restrict residual (computed on current MG level) to
3572 /// next coarser mesh and stick it into the coarse mesh RHS vector
3573 /// using the restriction matrix (if restrict_flag=1) or the transpose
3574 /// of the interpolation matrix (if restrict_flag=2)
3575 //===================================================================
3576  template<unsigned DIM>
3578  {
3579 #ifdef PARANOID
3580  // Check to make sure we can actually restrict the vector
3581  if (!(level<Nlevel-1))
3582  {
3583  // Throw an error
3584  throw OomphLibError("Input level exceeds the possible parameter choice.",
3585  OOMPH_CURRENT_FUNCTION,
3586  OOMPH_EXCEPTION_LOCATION);
3587  }
3588 #endif
3589 
3590  // Multiply the real part of the residual vector by the restriction
3591  // matrix on the level-th level
3592  Restriction_matrices_storage_pt[level]->
3593  multiply(Residual_mg_vectors_storage[level][0],
3594  Rhs_mg_vectors_storage[level+1][0]);
3595 
3596  // Multiply the imaginary part of the residual vector by the restriction
3597  // matrix on the level-th level
3598  Restriction_matrices_storage_pt[level]->
3599  multiply(Residual_mg_vectors_storage[level][1],
3600  Rhs_mg_vectors_storage[level+1][1]);
3601 
3602  } // End of restrict_residual
3603 
3604 //===================================================================
3605 /// \short Interpolate solution at current level onto
3606 /// next finer mesh and correct the solution x at that level
3607 //===================================================================
3608  template<unsigned DIM>
3610  interpolate_and_correct(const unsigned& level)
3611  {
3612 #ifdef PARANOID
3613  // Check to make sure we can actually restrict the vector
3614  if (!(level>0))
3615  {
3616  // Throw an error
3617  throw OomphLibError("Input level exceeds the possible parameter choice.",
3618  OOMPH_CURRENT_FUNCTION,
3619  OOMPH_EXCEPTION_LOCATION);
3620  }
3621 #endif
3622 
3623  // Build distribution of a temporary vector (real part)
3624  DoubleVector temp_soln_r(X_mg_vectors_storage[level-1][0].distribution_pt());
3625 
3626  // Build distribution of a temporary vector (imaginary part)
3627  DoubleVector temp_soln_c(X_mg_vectors_storage[level-1][1].distribution_pt());
3628 
3629  // Interpolate the solution vector
3630  Interpolation_matrices_storage_pt[level-1]->
3631  multiply(X_mg_vectors_storage[level][0],temp_soln_r);
3632 
3633  // Interpolate the solution vector
3634  Interpolation_matrices_storage_pt[level-1]->
3635  multiply(X_mg_vectors_storage[level][1],temp_soln_c);
3636 
3637  // Update the real part of the solution
3638  X_mg_vectors_storage[level-1][0]+=temp_soln_r;
3639 
3640  // Update the imaginary part of the solution
3641  X_mg_vectors_storage[level-1][1]+=temp_soln_c;
3642 
3643  } // End of interpolate_and_correct
3644 
3645 //===================================================================
3646 /// \short Linear solver. This is where the general V-cycle algorithm
3647 /// is implemented
3648 //===================================================================
3649  template<unsigned DIM>
3651  {
3652  // If we're allowed to time things
3653  double t_mg_start=0.0;
3654  if (!Suppress_v_cycle_output)
3655  {
3656  // Start the clock!
3657  t_mg_start=TimingHelpers::timer();
3658  }
3659 
3660  // Current level
3661  unsigned finest_level=0;
3662 
3663  // Initialise the V-cycle counter
3664  V_cycle_counter=0;
3665 
3666  // Calculate the norm of the residual then output
3667  double normalised_residual_norm=residual_norm(finest_level);
3668  if (!Suppress_v_cycle_output)
3669  {
3670  oomph_info << "\nResidual on finest level for V-cycle: "
3671  << normalised_residual_norm << std::endl;
3672  }
3673 
3674  // Outer loop over V-cycles
3675  //-------------------------
3676  // While the tolerance is not met and the maximum number of
3677  // V-cycles has not been completed
3678  while((normalised_residual_norm>Tolerance) &&
3679  (V_cycle_counter!=Nvcycle))
3680  {
3681  // If the user did not wish to suppress the V-cycle output
3682  if (!Suppress_v_cycle_output)
3683  {
3684  // Output the V-cycle counter
3685  oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
3686  }
3687 
3688  //---------------------------------------------------------------------
3689  // Loop downwards over all levels that have coarser levels beneath them
3690  //---------------------------------------------------------------------
3691  for (unsigned i=0;i<Nlevel-1;i++)
3692  {
3693  // Initialise X_mg and Residual_mg to 0.0 except for the finest level
3694  // since X_mg contains the current approximation to the solution and
3695  // Residual_mg contains the RHS vector on the finest level
3696  if(i!=0)
3697  {
3698  // Initialise the real part of the solution vector
3699  X_mg_vectors_storage[i][0].initialise(0.0);
3700 
3701  // Initialise the imaginary part of the solution vector
3702  X_mg_vectors_storage[i][1].initialise(0.0);
3703 
3704  // Initialise the real part of the residual vector
3705  Residual_mg_vectors_storage[i][0].initialise(0.0);
3706 
3707  // Initialise the imaginary part of the residual vector
3708  Residual_mg_vectors_storage[i][1].initialise(0.0);
3709  }
3710 
3711  // Perform a few pre-smoothing steps and return vector that contains
3712  // the residuals of the linear system at this level.
3713  pre_smooth(i);
3714 
3715  // Restrict the residual to the next coarser mesh and
3716  // assign it to the RHS vector at that level
3717  restrict_residual(i);
3718 
3719  } // Moving down the V-cycle
3720 
3721  // // Output the matrix and rhs vector that we want
3722  // // Create an output stream
3723  // std::ofstream outfile;
3724 
3725  // // Create an output name
3726  // std::string file1="matrix_r.dat";
3727  // outfile.open(file1);
3728  // Mg_matrices_storage_pt[Nlevel-2][0]->sparse_indexed_output_helper(outfile);
3729  // outfile.close();
3730 
3731  // std::string file2="matrix_c.dat";
3732  // outfile.open(file2);
3733  // Mg_matrices_storage_pt[Nlevel-2][1]->sparse_indexed_output_helper(outfile);
3734  // outfile.close();
3735 
3736  // std::string rhs_r="rhs_r.dat";
3737  // Rhs_mg_vectors_storage[Nlevel-2][0].output(rhs_r);
3738  // std::string rhs_c="rhs_c.dat";
3739  // Rhs_mg_vectors_storage[Nlevel-2][1].output(rhs_c);
3740 
3741  // // Exit!
3742  //exit(0);
3743 
3744  //-----------------------------------------------------------
3745  // Reached the lowest level: Do a direct solve, using the RHS
3746  // vector obtained by restriction from above.
3747  //-----------------------------------------------------------
3748  direct_solve();
3749 
3750  //---------------------------------------------------------------
3751  // Loop upwards over all levels that have finer levels above them
3752  //---------------------------------------------------------------
3753  for (unsigned i=Nlevel-1;i>0;i--)
3754  {
3755  // Interpolate solution at current level onto
3756  // next finer mesh and correct the solution x at that level
3757  interpolate_and_correct(i);
3758 
3759  // Perform a few post-smoothing steps (ignore
3760  // vector that contains the residuals of the linear system
3761  // at this level)
3762  post_smooth(i-1);
3763  }
3764 
3765  // Set counter for number of cycles (for doc)
3766  V_cycle_counter++;
3767 
3768  // Calculate the new residual norm then output (if allowed)
3769  normalised_residual_norm=residual_norm(finest_level);
3770 
3771  // Print the residual on the finest level
3772  if (!Suppress_v_cycle_output)
3773  {
3774  oomph_info << "Residual on finest level of V-cycle: "
3775  << normalised_residual_norm << std::endl;
3776  }
3777  } // End of the V-cycles
3778 
3779  // Copy the solution into the result vector
3780  result=X_mg_vectors_storage[finest_level];
3781 
3782  // Need an extra line space if V-cycle output is suppressed
3783  if (!Suppress_v_cycle_output)
3784  {
3785  oomph_info << std::endl;
3786  }
3787 
3788  // If all output is to be suppressed
3789  if (!Suppress_all_output)
3790  {
3791  // Output number of V-cycles taken to solve
3792  if (normalised_residual_norm<Tolerance)
3793  {
3794  Has_been_solved=true;
3795  }
3796  } // if (!Suppress_all_output)
3797 
3798  // If the V-cycle output isn't suppressed
3799  if (!Suppress_v_cycle_output)
3800  {
3801  // Stop the clock
3802  double t_mg_end=TimingHelpers::timer();
3803  double total_G_setup_time=double(t_mg_end-t_mg_start);
3804  oomph_info << "CPU time for MG solver [sec]: "
3805  << total_G_setup_time << std::endl;
3806  }
3807  } // end of mg_solve
3808 
3809 } // End of namespace oomph
3810 
3811 #endif
double Alpha_shift
Temporary version of the shift – to run bash scripts.
void setup()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level...
int * column_index()
Access to C-style column index array.
Definition: matrices.h:1056
virtual void actions_before_adapt()
Actions that are to be performed before a mesh adaptation. These might include removing any additiona...
Definition: problem.h:1004
void 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 co...
void setup_mg_hierarchy()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level...
double & tolerance()
Access to convergence tolerance.
Nullstream oomph_nullstream
Single (global) instantiation of the Nullstream.
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
void enable_doc_time()
Enable time documentation.
The GMRES method rewritten for complex matrices.
void pre_smooth(const unsigned &level)
Pre-smoother: Perform 'max_iter' smoothing steps on the linear system Ax=b with current RHS vector...
void setup_smoothers()
Function to set up all of the smoothers once the system matrices have been set up.
unsigned Nvcycle
Maximum number of V-cycles.
void interpolation_matrix_set(const unsigned &level, Vector< double > &value, Vector< int > &col_index, Vector< int > &row_st, unsigned &ncol, unsigned &nrow)
Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be...
cstr elem_len * i
Definition: cfortran.h:607
HelmholtzMGPreconditioner(HelmholtzMGProblem *mg_problem_pt)
Constructor: Set up default values for number of V-cycles and pre- and post-smoothing steps...
Vector< Vector< DoubleVector > > X_mg_vectors_storage
Vector of vectors to store the solution vectors (X_mg) in two parts; the real and imaginary...
void setup_transfer_matrices()
Setup the transfer matrices on each level.
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
virtual TreeBasedRefineableMeshBase * mg_bulk_mesh_pt()=0
Function to get a pointer to the mesh we will be working with. If there are flux elements present in ...
PreSmootherFactoryFctPt Pre_smoother_factory_function_pt
Function to create pre-smoothers.
const double Pi
50 digits from maple
unsigned nrow() const
access function to the number of global rows.
virtual void actions_after_adapt()
Actions that are to be performed after a mesh adaptation.
Definition: problem.h:1007
bool Has_been_solved
Boolean variable to indicate whether or not the problem was successfully solved.
The Problem class.
Definition: problem.h:152
double Wavenumber
The value of the wavenumber, k.
unsigned iterations() const
Number of iterations.
unsigned Npost_smooth
Number of post-smoothing steps.
bool Doc_time
Indicates whether or not time documentation should be used.
A general Finite Element class.
Definition: elements.h:1271
double & tolerance()
Access function for the variable Tolerance (lvalue)
Vector< Vector< CRDoubleMatrix * > > Mg_matrices_storage_pt
Vector of vectors to store the system matrices. The i-th entry in this vector contains a vector of le...
void multiply(const CRDoubleMatrix *matrix, const DoubleVector &x, DoubleVector &soln)
Function to perform a matrix-vector multiplication on a oomph-lib matrix and vector using Trilinos fu...
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
Find the sub geometric object and local coordinate therein that corresponds to the intrinsic coordina...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
bool Suppress_all_output
If this is set to true then all output from the solver is suppressed.
Vector< CRDoubleMatrix * > Restriction_matrices_storage_pt
Vector to store the restriction matrices.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
OomphInfo oomph_info
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
double k_squared()
Get the square of wavenumber.
void interpolate_and_correct(const unsigned &level)
Interpolate solution at current level onto next finer mesh and correct the solution x at that level...
e
Definition: cfortran.h:575
double * value()
Access to C-style value array.
Definition: matrices.h:1062
void setup_mg_structures()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level...
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
void block_preconditioner_self_test()
Function to ensure the block form of the Jacobian matches the form described, i.e. we should have: |--—|---—| | A_r | -A_c | A = |--—|---—|. | A_c | A_r | |--—|---—|.
Vector< HelmholtzSmoother * > Post_smoothers_storage_pt
Vector to store the post-smoothers.
double residual_norm(const unsigned &level)
Return norm of residual r=b-Ax and the residual vector itself on the level-th level.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)...
Definition: matrices.cc:55
void mg_solve(Vector< DoubleVector > &result)
Do the actual solve – this is called through the pure virtual solve function in the LinearSolver base...
void set_restriction_matrices_as_interpolation_transposes()
Builds a CRDoubleMatrix on each level that is used to restrict the residual between levels...
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
void disable_smoother_and_superlu_doc_time()
Suppress the output of both smoothers and SuperLU.
Vector< HelmholtzMGProblem * > Mg_hierarchy_pt
Vector containing pointers to problems in hierarchy.
DoubleVector Coarsest_rhs_mg
Assuming we're solving the system Ax=b, this vector will contain the expanded solution vector on the ...
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1068
LinearSolver *& linear_solver_pt()
Return a pointer to the linear solver object.
Definition: matrices.h:319
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void setup_coarsest_level_structures()
Function to create the fully expanded system matrix on the coarsest level.
void disable_v_cycle_output()
Disable all output from mg_solve apart from the number of V-cycles used to solve the problem...
virtual HelmholtzMGProblem * make_new_problem()=0
This function needs to be implemented in the derived problem: Returns a pointer to a new object of th...
void enable_v_cycle_output()
Enable the output of the V-cycle timings and other output.
void disable_doc_time()
Disable time documentation.
void maximum_edge_width(const unsigned &level, double &h)
Estimate the value of the parameter h on the level-th problem in the hierarchy.
std::ostream *& stream_pt()
Access function for the stream pointer.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:998
Vector< Vector< DoubleVector > > Rhs_mg_vectors_storage
Vector of vectors to store the RHS vectors. This uses the same format as the X_mg_vectors_storage vec...
Base class for tree-based refineable meshes.
void setup_interpolation_matrices()
Setup the interpolation matrix on each level.
HelmholtzMGProblem class; subclass of Problem.
static char t char * s
Definition: cfortran.h:572
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:992
void post_smooth(const unsigned &level)
Post-smoother: Perform max_iter smoothing steps on the linear system Ax=b with current RHS vector...
void interpolation_matrix_set(const unsigned &level, double *value, int *col_index, int *row_st, unsigned &ncol, unsigned &nnz)
Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be...
unsigned & npre_smooth()
Return the number of pre-smoothing iterations (lvalue)
DoubleVector Coarsest_x_mg
Assuming we're solving the system Ax=b, this vector will contain the expanded solution vector on the ...
Vector< HelmholtzSmoother * > Pre_smoothers_storage_pt
Vector to store the pre-smoothers.
PostSmootherFactoryFctPt Post_smoother_factory_function_pt
Function to create post-smoothers.
void disable_doc_time()
Disable documentation of solve times.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void setup_interpolation_matrices_unstructured()
Setup the interpolation matrix on each level (used for unstructured meshes)
double timer()
returns the time in seconds after some point in past
void initialise(const _Tp &__value)
Iterate over all values and set to the desired value.
Definition: Vector.h:171
Class that contains data for hanging nodes.
Definition: nodes.h:684
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
double *& alpha_pt()
Pointer to Alpha, wavenumber complex shift.
HelmholtzMGProblem()
Constructor. Initialise pointers to coarser and finer levels.
Vector< Vector< DoubleVector > > Residual_mg_vectors_storage
Vector to vectors to store the residual vectors. This uses the same format as the X_mg_vectors_storag...
void disable_output()
Suppress anything that can be suppressed, i.e. any timings. Things like mesh adaptation can not howev...
void restrict_residual(const unsigned &level)
Restrict residual (computed on level-th MG level) to the next coarser mesh and stick it into the coar...
void full_setup()
Do a full setup (assumes everything will be setup around the HelmholtzMGProblem pointer given in the ...
HelmholtzMGProblem * Mg_problem_pt
Pointer to the MG problem (deep copy)
HelmholtzSmoother *(* PreSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class HelmholtzSmoother to be used ...
void calculate_omega(const double &k, const double &h)
Function to calculate the value of Omega by passing in the value of k and h [see Elman et al...
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:241
void enable_output()
Enable the output from anything that could have been suppressed.
double Tolerance
The tolerance of the multigrid preconditioner.
void set_pre_smoother_factory_function(PreSmootherFactoryFctPt pre_smoother_fn)
Access function to set the pre-smoother creation function.
double & alpha_shift()
Function to change the value of the shift.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned V_cycle_counter
Pointer to counter for V-cycles.
bool Has_been_setup
Boolean variable to indicate whether or not the solver has been setup.
void return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
Takes the vector of block vectors, s, and copies its entries into the naturally ordered vector...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
HelmholtzSmoother *(* PostSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class HelmholtzSmoother to be used ...
void split(const DoubleVector &in_vector, Vector< DoubleVector * > &out_vector_pt)
Split a DoubleVector into the out DoubleVectors. Let vec_A be the in Vector, and let vec_B and vec_C ...
virtual void setup()=0
Setup the preconditioner. Pure virtual generic interface function.
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
bool Suppress_v_cycle_output
Indicates whether or not the V-cycle output should be suppressed.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
unsigned Nlevel
The number of levels in the multigrid heirachy.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout.
Vector< CRDoubleMatrix * > Interpolation_matrices_storage_pt
Vector to store the interpolation matrices.
unsigned & npost_smooth()
Return the number of post-smoothing iterations (lvalue)
CRDoubleMatrix * Coarsest_matrix_mg_pt
Stores the system matrix on the coarest level in the fully expanded format: |--—|---—| | A_r | -A_c |...
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1050
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Function applies MG to the vector r for a full solve.
Vector< double > Max_edge_width
Vector to storage the maximum edge width of each mesh. We only need the maximum edge width on levels ...
void concatenate(const Vector< DoubleVector * > &in_vector_pt, DoubleVector &out_vector)
Concatenate DoubleVectors. Takes a Vector of DoubleVectors. If the out vector is built, we will not build a new distribution. Otherwise we build a uniform distribution.
~HelmholtzMGPreconditioner()
Delete any dynamically allocated data.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
void set_post_smoother_factory_function(PostSmootherFactoryFctPt post_smoother_fn)
Access function to set the post-smoother creation function.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
A general mesh class.
Definition: mesh.h:74
virtual ~HelmholtzMGProblem()
Destructor (empty)
void build(const LinearAlgebraDistribution *distribution_pt, const unsigned &ncol, const Vector< double > &value, const Vector< int > &column_index, const Vector< int > &row_start)
build method: vector of values, vector of column indices, vector of row starts and number of rows and...
Definition: matrices.cc:1692
void clean_up_memory()
Clean up anything that needs to be cleaned up.
unsigned Npre_smooth
Number of pre-smoothing steps.
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2399
void direct_solve()
Call the direct solver (SuperLU) to solve the problem exactly.