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