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_GEOMETRIC_MULTIGRID_HEADER
28#define OOMPH_GEOMETRIC_MULTIGRID_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35// For the Problem class
36#include "problem.h"
37
38// For the TreeBasedRefineableMeshBase and MeshAsGeomObject class
39#include "refineable_mesh.h"
41
42// For the RefineableQElement class
43#include "Qelements.h"
44
45// Other obvious stuff
46#include "matrices.h"
48#include "preconditioner.h"
49
50// Namespace extension
51namespace oomph
52{
53 //======================================================================
54 /// MGProblem class; subclass of Problem
55 //======================================================================
56 class MGProblem : public virtual Problem
57 {
58 public:
59 /// Constructor. Initialise pointers to coarser and finer levels
61
62 /// Destructor (empty)
63 virtual ~MGProblem() {}
64
65 /// This function needs to be implemented in the derived problem:
66 /// Returns a pointer to a new object of the same type as the derived
67 /// problem
69
70 /// Function to get a pointer to the mesh we will be working
71 /// with. If there are flux elements present in the mesh this will
72 /// be overloaded to return a pointer to the bulk mesh otherwise
73 /// it can be overloaded to point to the global mesh but it must
74 /// be of type RefineableMeshBase
76
77 }; // End of MGProblem class
78
79
80 /// ///////////////////////////////////////////////////////
81 /// ///////////////////////////////////////////////////////
82 /// ///////////////////////////////////////////////////////
83
84
85 //======================================================================
86 // MG solver class
87 //======================================================================
88 template<unsigned DIM>
90 {
91 public:
92 /// typedef for a function that returns a pointer to an object
93 /// of the class Smoother to be used as the pre-smoother
94 typedef Smoother* (*PreSmootherFactoryFctPt)();
95
96 /// typedef for a function that returns a pointer to an object
97 /// of the class Smoother to be used as the post-smoother
98 typedef Smoother* (*PostSmootherFactoryFctPt)();
99
100 /// Access function to set the pre-smoother creation function.
102 PreSmootherFactoryFctPt pre_smoother_fn)
103 {
104 // Assign the function pointer
105 Pre_smoother_factory_function_pt = pre_smoother_fn;
106 }
107
108 /// Access function to set the post-smoother creation function.
110 PostSmootherFactoryFctPt post_smoother_fn)
111 {
112 // Assign the function pointer
113 Post_smoother_factory_function_pt = post_smoother_fn;
114 }
115
116 /// Constructor: Set up default values for number of V-cycles
117 /// and pre- and post-smoothing steps.
118 MGSolver(MGProblem* mg_problem_pt)
119 : Nvcycle(200),
120 Mg_problem_pt(mg_problem_pt),
122 Suppress_all_output(false),
123 Stream_pt(0),
126 Npre_smooth(2),
127 Npost_smooth(2),
128 Doc_everything(false),
129 Has_been_setup(false),
130 Has_been_solved(false)
131 {
132 // Set the tolerance in the base class
133 this->Tolerance = 1.0e-09;
134
135 // Set the pointer to the finest level as the first
136 // entry in Mg_hierarchy
137 Mg_hierarchy.push_back(Mg_problem_pt);
138 } // End of MGSolver
139
140 /// Delete any dynamically allocated data
142 {
143 // Run the function written to clean up the memory
145 } // End of ~MGSolver
146
147 /// Clean up anything that needs to be cleaned up
149 {
150 // We only need to destroy data if the solver has been set up and
151 // the data hasn't already been cleared
152 if (Has_been_setup)
153 {
154 // Loop over all of the levels in the hierarchy
155 for (unsigned i = 0; i < Nlevel - 1; i++)
156 {
157 // Delete the pre-smoother associated with this level
159
160 // Make it a null pointer
162
163 // Delete the post-smoother associated with this level
165
166 // Make it a null pointer
168
169 // Delete the system matrix associated with the i-th level
171
172 // Make it a null pointer
174 }
175
176 // Loop over all but the coarsest of the levels in the hierarchy
177 for (unsigned i = 0; i < Nlevel - 1; i++)
178 {
179 // Delete the interpolation matrix associated with the i-th level
181
182 // Make it a null pointer
184
185 // Delete the restriction matrix associated with the i-th level
187
188 // Make it a null pointer
190 }
191
192 // If this solver has been set up then a hierarchy of problems
193 // will have been set up. If the user chose to document everything
194 // then the coarse-grid multigrid problems will have been kept alive
195 // which means we now have to loop over the coarse-grid levels and
196 // destroy them
197 if (Doc_everything)
198 {
199 // Loop over the levels
200 for (unsigned i = 1; i < Nlevel; i++)
201 {
202 // Delete the i-th level problem
203 delete Mg_hierarchy[i];
204
205 // Make the associated pointer a null pointer
206 Mg_hierarchy[i] = 0;
207 }
208 } // if (Doc_everything)
209
210 // Everything has been deleted now so we need to indicate that the
211 // solver is not set up
212 Has_been_setup = false;
213 } // if (Has_been_setup)
214 } // End of clean_up_memory
215
216 /// Makes a vector which will be used in the self-test. Is currently
217 /// set to make the entries of the vector represent a plane wave propagating
218 /// at an angle of 45 degrees
220
221 /// Makes a vector, restricts it down the levels of the hierarchy
222 /// and documents it at each level. After this is done the vector is
223 /// interpolated up the levels of the hierarchy with the output
224 /// being documented at each level
225 void self_test();
226
227 /// Make a self-test to make sure that the interpolation matrices
228 /// are doing the same thing to restrict the vectors down through the
229 /// heirachy.
231
232 /// Make a self-test to make sure that the interpolation matrices
233 /// are doing the same thing to interpolate the vectors up.
235
236 /// Given a level in the hierarchy, an input vector and a filename
237 /// this function will document the given vector according to the structure
238 /// of the mesh on the given level
239 void plot(const unsigned& hierarchy_level,
240 const DoubleVector& input_vector,
241 const std::string& filename);
242
243 /// Disable all output from mg_solve apart from the number of
244 /// V-cycles used to solve the problem
246 {
247 // Set the value of Doc_time (inherited from LinearSolver) to false
248 Doc_time = false;
249
250 // Enable the suppression of output from the V-cycle
252 } // End of disable_v_cycle_output
253
254 /// Suppress anything that can be suppressed, i.e. any timings.
255 /// Things like mesh adaptation can not however be silenced using this
257 {
258 // Set the value of Doc_time (inherited from LinearSolver) to false
259 Doc_time = false;
260
261 // Enable the suppression of output from the V-cycle
263
264 // Enable the suppression of everything
265 Suppress_all_output = true;
266
267 // Store the output stream pointer
269
270 // Now set the oomph_info stream pointer to the null stream to
271 // disable all possible output
273 } // End of disable_output
274
275 /// Enable the output of the V-cycle timings and other output
277 {
278 // Enable time documentation
279 Doc_time = true;
280
281 // Enable output from the MG solver
283 } // End of enable_v_cycle_output
284
285 /// Enable the output from anything that could have been suppressed
287 {
288 // Enable the documentation of everything (if this is set to TRUE then
289 // the function self_test() will be run which outputs a solution
290 // represented on each level of the hierarchy
291 Doc_everything = true;
292 } // End of enable_doc_everything
293
294 /// Enable the output from anything that could have been suppressed
296 {
297 // Enable time documentation
298 Doc_time = true;
299
300 // Enable output from everything during the full setup of the solver
301 Suppress_all_output = false;
302
303 // Enable output from the MG solver
305 } // End of enable_output
306
307 /// Suppress the output of both smoothers and SuperLU
309 {
310 // Loop over all levels of the hierarchy
311 for (unsigned i = 0; i < Nlevel - 1; i++)
312 {
313 // Disable time documentation on each level (for each pre-smoother)
314 Pre_smoothers_storage_pt[i]->disable_doc_time();
315
316 // Disable time documentation on each level (for each post-smoother)
317 Post_smoothers_storage_pt[i]->disable_doc_time();
318 }
319
320 // We only do a direct solve on the coarsest level so this is the only
321 // place we need to silence SuperLU
323 ->linear_solver_pt()
324 ->disable_doc_time();
325 } // End of disable_smoother_and_superlu_doc_time
326
327 /// Return the number of post-smoothing iterations (lvalue)
328 unsigned& npost_smooth()
329 {
330 // Return the number of post-smoothing iterations to be done on each
331 // level of the hierarchy
332 return Npost_smooth;
333 } // End of npost_smooth
334
335 /// Return the number of pre-smoothing iterations (lvalue)
336 unsigned& npre_smooth()
337 {
338 // Return the number of pre-smoothing iterations to be done on each
339 // level of the hierarchy
340 return Npre_smooth;
341 } // End of npre_smooth
342
343 /// Pre-smoother: Perform 'max_iter' smoothing steps on the
344 /// linear system Ax=b with current RHS vector, b, starting with
345 /// current solution vector, x. Return the residual vector r=b-Ax.
346 /// Uses the default smoother (set in the MGProblem constructor)
347 /// which can be overloaded for a specific problem.
348 void pre_smooth(const unsigned& level)
349 {
350 // Run pre-smoother 'max_iter' times
351 Pre_smoothers_storage_pt[level]->smoother_solve(
353
354 // Calculate the residual r=b-Ax and assign it
355 Mg_matrices_storage_pt[level]->residual(
359 } // End of pre_smooth
360
361 /// Post-smoother: Perform max_iter smoothing steps on the
362 /// linear system Ax=b with current RHS vector, b, starting with
363 /// current solution vector, x. Uses the default smoother (set in
364 /// the MGProblem constructor) which can be overloaded for specific
365 /// problem.
366 void post_smooth(const unsigned& level)
367 {
368 // Run post-smoother 'max_iter' times
369 Post_smoothers_storage_pt[level]->smoother_solve(
371 } // End of post_smooth
372
373 /// Return norm of residual r=b-Ax and the residual vector itself
374 /// on the level-th level
375 double residual_norm(const unsigned& level)
376 {
377 // And zero the entries of residual
378 Residual_mg_vectors_storage[level].initialise(0.0);
379
380 // Get the residual
381 Mg_matrices_storage_pt[level]->residual(
385
386 // Return the norm of the residual
387 return Residual_mg_vectors_storage[level].norm();
388 } // End of residual_norm
389
390 /// Call the direct solver (SuperLU) to solve the problem exactly.
391 // The result is placed in X_mg
393 {
394 // Get solution by direct solve:
395 Mg_matrices_storage_pt[Nlevel - 1]->solve(
397 } // End of direct_solve
398
399 /// Builds a CRDoubleMatrix that is used to interpolate the
400 /// residual between levels. The transpose can be used as the full
401 /// weighting restriction.
402 void interpolation_matrix_set(const unsigned& level,
403 double* value,
404 int* col_index,
405 int* row_st,
406 unsigned& ncol,
407 unsigned& nnz)
408 {
409 // Dynamically allocate the interpolation matrix
411
412 // Build the matrix
413 Interpolation_matrices_storage_pt[level]->build_without_copy(
414 ncol, nnz, value, col_index, row_st);
415 } // End of interpolation_matrix_set
416
417 /// Builds a CRDoubleMatrix that is used to interpolate the
418 /// residual between levels. The transpose can be used as the full
419 /// weighting restriction.
420 void interpolation_matrix_set(const unsigned& level,
421 Vector<double>& value,
422 Vector<int>& col_index,
423 Vector<int>& row_st,
424 unsigned& ncol,
425 unsigned& nrow)
426 {
427 // Dynamically allocate the interpolation matrix
429
430 // Make the distribution pointer
432 Mg_hierarchy[level]->communicator_pt(), nrow, false);
433
434#ifdef PARANOID
435#ifdef OOMPH_HAS_MPI
436 // Set up the warning messages
437 std::string warning_message =
438 "Setup of interpolation matrix distribution ";
439 warning_message += "has not been tested with MPI.";
440
441 // If we're not running the code in serial
442 if (dist_pt->communicator_pt()->nproc() > 1)
443 {
444 // Throw a warning
446 warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
447 }
448#endif
449#endif
450
451 // Build the matrix itself
453 dist_pt, ncol, value, col_index, row_st);
454
455 // Delete the newly created distribution pointer
456 delete dist_pt;
457
458 // Make it a null pointer
459 dist_pt = 0;
460 } // End of interpolation_matrix_set
461
462 /// Builds a CRDoubleMatrix on each level that is used to
463 /// restrict the residual between levels. The transpose can be used
464 /// as the interpolation matrix
466 {
467 for (unsigned i = 0; i < Nlevel - 1; i++)
468 {
469 // Dynamically allocate the restriction matrix
471
472 // Set the restriction matrix to be the transpose of the
473 // interpolation matrix
474 Interpolation_matrices_storage_pt[i]->get_matrix_transpose(
476 }
477 } // End of set_restriction_matrices_as_interpolation_transposes
478
479 /// Restrict residual (computed on level-th MG level) to the next
480 /// coarser mesh and stick it into the coarse mesh RHS vector.
481 void restrict_residual(const unsigned& level);
482
483 /// Interpolate solution at current level onto next finer mesh
484 /// and correct the solution x at that level
485 void interpolate_and_correct(const unsigned& level);
486
487 /// Given the son_type of an element and a local node number
488 /// j in that element with nnode_1d nodes per coordinate direction,
489 /// return the local coordinate s in its father element. Needed in
490 /// the setup of the interpolation matrices
491 void level_up_local_coord_of_node(const int& son_type, Vector<double>& s);
492
493 /// Setup the interpolation matrix on each level
495
496 /// Setup the interpolation matrix on each level (used for
497 /// unstructured meshes)
499
500 /// Setup the transfer matrices on each level
502
503 /// Do a full setup (assumes everything will be setup around the
504 /// MGProblem pointer given in the constructor)
505 void full_setup();
506
507 /// Virtual function in the base class that needs to be implemented
508 /// later but for now just leave it empty
509 void solve(Problem* const& problem_pt, DoubleVector& result)
510 {
511 // Dynamically cast problem_pt of type Problem to a MGProblem pointer
512 MGProblem* mg_problem_pt = dynamic_cast<MGProblem*>(problem_pt);
513
514#ifdef PARANOID
515 // PARANOID check - If the dynamic_cast produces a null pointer the
516 // input was not a MGProblem
517 if (0 == mg_problem_pt)
518 {
519 throw OomphLibError("Input problem must be of type MGProblem.",
520 OOMPH_CURRENT_FUNCTION,
521 OOMPH_EXCEPTION_LOCATION);
522 }
523 // PARANOID check - If a node in the input problem has more than
524 // one value we cannot deal with it so arbitarily check the first one
525 if (problem_pt->mesh_pt()->node_pt(0)->nvalue() != 1)
526 {
527 // How many dofs are there in the first node
528 unsigned n_value = problem_pt->mesh_pt()->node_pt(0)->nvalue();
529
530 // Make the error message
531 std::ostringstream error_message_stream;
532 error_message_stream << "Cannot currently deal with more than 1 dof"
533 << " per node. This problem has " << n_value
534 << " dofs per node." << std::endl;
535
536 // Throw the error message
537 throw OomphLibError(error_message_stream.str(),
538 OOMPH_CURRENT_FUNCTION,
539 OOMPH_EXCEPTION_LOCATION);
540 }
541#endif
542
543 // Assign the new MGProblem pointer to Mg_problem_pt
544 Mg_problem_pt = mg_problem_pt;
545
546 // Set up all of the required MG structures
547 full_setup();
548
549 // Run the MG method and assign the solution to result
550 mg_solve(result);
551
552 // Only output if the V-cycle output isn't suppressed
554 {
555 // Notify user that the hierarchy of levels is complete
556 oomph_info << "\n================="
557 << "Multigrid Solve Complete"
558 << "=================\n"
559 << std::endl;
560 }
561
562 // If the user did not request all output be suppressed
564 {
565 // If the user requested all V-cycle output be suppressed
567 {
568 // Create an extra line spacing
569 oomph_info << std::endl;
570 }
571
572 // Output number of V-cycles taken to solve
573 if (Has_been_solved)
574 {
575 oomph_info << "Total number of V-cycles required for solve: "
576 << V_cycle_counter << std::endl;
577 }
578 else
579 {
580 oomph_info << "Total number of V-cycles used: " << V_cycle_counter
581 << std::endl;
582 }
583 } // if (!Suppress_all_output)
584
585 // Only enable and assign the stream pointer again if we originally
586 // suppressed everything otherwise it won't be set yet
588 {
589 // Now enable the stream pointer again
591 }
592 } // End of solve
593
594 /// Number of iterations
595 unsigned iterations() const
596 {
597 // Return the number of V-cycles which have been done
598 return V_cycle_counter;
599 } // End of iterations
600
601 /// Number of iterations
602 unsigned& max_iter()
603 {
604 // Return the number of V-cycles which have been done
605 return Nvcycle;
606 } // End of iterations
607
608 protected:
609 /// Do the actual solve -- this is called through the pure virtual
610 /// solve function in the LinearSolver base class. The function is stored
611 /// as protected to allow the MGPreconditioner derived class to use the
612 /// solver
613 void mg_solve(DoubleVector& result);
614
615 /// Normalise the rows of the restriction matrices to avoid
616 /// amplifications when projecting to the coarser level
618
619 /// Maximum number of V-cycles (this is set as a protected variable
620 /// so
621 // that it can be changed in the MGPreconditioner class)
622 unsigned Nvcycle;
623
624 /// Pointer to the MG problem (deep copy). This is protected to
625 /// provide access to the MG preconditioner
627
628 /// Vector to store the RHS vectors (Rhs_mg). This is protected to
629 /// allow the multigrid preconditioner to assign the RHS vector during
630 /// preconditioner_solve()
632
633 /// Indicates whether or not the V-cycle output should be
634 /// suppressed. Needs to be protected member data for the multigrid
635 /// preconditioner to know whether or not to output information
636 /// with each preconditioning step
638
639 /// If this is set to true then all output from the solver is
640 /// suppressed. This is protected member data so that the multigrid
641 /// preconditioner knows whether or not to restore the stream pointer
643
644 /// Pointer to the output stream -- defaults to std::cout.
645 /// This is protected member data to allow the preconditioner to
646 /// restore normal output if everything was chosen to be
647 /// suppressed by the user
648 std::ostream* Stream_pt;
649
650 private:
651 /// Function to create pre-smoothers
653
654 /// Function to create post-smoothers
656
657 /// Function to set up the hierachy of levels. Creates a vector
658 /// of pointers to each MG level
659 void setup_mg_hierarchy();
660
661 /// Function to set up the hierachy of levels. Creates a vector
662 /// of pointers to each MG level
663 void setup_mg_structures();
664
665 /// Function to set up all of the smoothers once the system matrices
666 /// have been set up
667 void setup_smoothers();
668
669 /// The number of levels in the multigrid heirachy
670 unsigned Nlevel;
671
672 /// Vector containing pointers to problems in hierarchy
674
675 /// Vector to store the system matrices
677
678 /// Vector to store the interpolation matrices
680
681 /// Vector to store the restriction matrices
683
684 /// Vector to store the solution vectors (X_mg)
686
687 /// Vector to store the residual vectors
689
690 /// Vector to store the result of interpolation on each level (only
691 /// required if the user wishes to document the output of interpolation
692 /// and restriction on each level)
694
695 /// Vector to store the result of restriction on each level (only
696 /// required if the user wishes to document the output of interpolation
697 /// and restriction on each level)
699
700 /// Vector to store the pre-smoothers
702
703 /// Vector to store the post-smoothers
705
706 /// Number of pre-smoothing steps
707 unsigned Npre_smooth;
708
709 /// Number of post-smoothing steps
710 unsigned Npost_smooth;
711
712 /// If this is set to true we document everything. In addition
713 /// to outputting the information of the setup timings and V-cycle
714 /// data we document the refinement and unrefinement patterns given
715 /// by the transfer operators which is done by keeping the coarser
716 /// MG problem pointers alive
718
719 /// Boolean variable to indicate whether or not the solver has been setup
721
722 /// Boolean variable to indicate whether or not the problem was
723 /// successfully solved
725
726 /// Pointer to counter for V-cycles
728 };
729
730 //====================================================================
731 /// An interface to allow scalar MG to be used as a Preconditioner
732 //====================================================================
733 template<unsigned DIM>
734 class MGPreconditioner : public MGSolver<DIM>, public Preconditioner
735 {
736 public:
737 /// Constructor.
738 MGPreconditioner(MGProblem* mg_problem_pt) : MGSolver<DIM>(mg_problem_pt)
739 {
740 // Set the number of V-cycles to be 1 (as expected as a preconditioner)
741 this->Nvcycle = 2;
742 } // End of MGPreconditioner (constructor)
743
744 /// Destructor (empty)
746
747 /// Broken copy constructor.
749
750 /// Broken assignment operator.
751 void operator=(const MGPreconditioner&) = delete;
752
753 /// Function to set up a preconditioner for the linear system
754 void setup()
755 {
756#ifdef OOMPH_HAS_MPI
757 // Make sure that this is running in serial. Can't guarantee it'll
758 // work when the problem is distributed over several processors
759 if (MPI_Helpers::communicator_pt()->nproc() > 1)
760 {
761 // Throw a warning
762 OomphLibWarning("Can't guarantee the MG solver will work in parallel!",
763 OOMPH_CURRENT_FUNCTION,
764 OOMPH_EXCEPTION_LOCATION);
765 }
766#endif
767
768 // Call the helper function that actually does all the work
769 this->full_setup();
770
771 // Only enable and assign the stream pointer again if we originally
772 // suppressed everything otherwise it won't be set yet
773 if (this->Suppress_all_output)
774 {
775 // Now enable the stream pointer again
777 }
778 } // End of setup
779
780 /// Function applies MG to the vector r for a full solve
781 virtual void preconditioner_solve(const DoubleVector& rhs, DoubleVector& z)
782 {
783#ifdef PARANOID
784 if (this->Mg_problem_pt->ndof() != rhs.nrow())
785 {
786 throw OomphLibError("Matrix and RHS vector sizes incompatible.",
787 OOMPH_CURRENT_FUNCTION,
788 OOMPH_EXCEPTION_LOCATION);
789 }
790#endif
791
792 // Set the right-hand side vector on the finest level to r
793 this->Rhs_mg_vectors_storage[0] = rhs;
794
795 // Run the MG method and assign the solution to z
796 this->mg_solve(z);
797
798 // Only output if the V-cycle output isn't suppressed
799 if (!(this->Suppress_v_cycle_output))
800 {
801 // Notify user that the hierarchy of levels is complete
803 << "\n==========Multigrid Preconditioner Solve Complete========="
804 << "\n"
805 << std::endl;
806 }
807
808 // Only enable and assign the stream pointer again if we originally
809 // suppressed everything otherwise it won't be set yet
810 if (this->Suppress_all_output)
811 {
812 // Now enable the stream pointer again
814 }
815 } // End of preconditioner_solve
816
817 /// Clean up memory
819 };
820
821
822 //===================================================================
823 /// Runs a full setup of the MG solver
824 //===================================================================
825 template<unsigned DIM>
827 {
828 // Initialise the timer start variable
829 double t_fs_start = 0.0;
830
831 // If we're allowed to output
832 if (!Suppress_all_output)
833 {
834 // Start the timer
835 t_fs_start = TimingHelpers::timer();
836
837 // Notify user that the hierarchy of levels is complete
839 << "\n===============Starting Multigrid Full Setup=============="
840 << std::endl;
841
842 // Notify user that the hierarchy of levels is complete
843 oomph_info << "\nStarting the full setup of the multigrid solver."
844 << std::endl;
845 }
846
847#ifdef PARANOID
848 // PARANOID check - Make sure the dimension of the solver matches the
849 // dimension of the elements used in the problems mesh
850 if (dynamic_cast<FiniteElement*>(Mg_problem_pt->mesh_pt()->element_pt(0))
851 ->dim() != DIM)
852 {
853 std::string err_strng = "The dimension of the elements used in the mesh ";
854 err_strng += "does not match the dimension of the solver.";
855 throw OomphLibError(
856 err_strng, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
857 }
858
859 // PARANOID check - The elements of the bulk mesh must all be refineable
860 // elements otherwise we cannot deal with this
861 if (Mg_problem_pt->mg_bulk_mesh_pt() != 0)
862 {
863 // Find the number of elements in the bulk mesh
864 unsigned n_elements = Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
865
866 // Loop over the elements in the mesh and ensure that they are
867 // all refineable elements
868 for (unsigned el_counter = 0; el_counter < n_elements; el_counter++)
869 {
870 // Upcast global mesh element to a refineable element
871 RefineableElement* el_pt = dynamic_cast<RefineableElement*>(
872 Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));
873
874 // Check if the upcast worked or not; if el_pt is a null pointer the
875 // element is not refineable
876 if (el_pt == 0)
877 {
878 throw OomphLibError(
879 "Element in global mesh could not be upcast to a refineable "
880 "element. We cannot deal with elements that are not refineable.",
881 OOMPH_CURRENT_FUNCTION,
882 OOMPH_EXCEPTION_LOCATION);
883 }
884 }
885 }
886 else
887 {
888 throw OomphLibError(
889 "The provided bulk mesh pointer is set to be a null pointer. "
890 "The multigrid solver operates on the bulk mesh thus a pointer "
891 "to the correct mesh must be given.",
892 OOMPH_CURRENT_FUNCTION,
893 OOMPH_EXCEPTION_LOCATION);
894 }
895#endif
896
897 // If this is not the first Newton step then we will already have things
898 // in storage. If this is the case, delete them
900
901 // Resize the Mg_hierarchy vector
902 Mg_hierarchy.resize(1, 0);
903
904 // Set the pointer to the finest level as the first entry in Mg_hierarchy
905 Mg_hierarchy[0] = Mg_problem_pt;
906
907 // Create the hierarchy of levels
908 setup_mg_hierarchy();
909
910 // Set up the interpolation and restriction matrices
911 setup_transfer_matrices();
912
913 // Set up the data structures on each level, i.e. the system matrix,
914 // LHS and RHS vectors
915 setup_mg_structures();
916
917 // Set up the smoothers on all of the levels
918 setup_smoothers();
919
920 // If we do not want to document everything we want to delete all the
921 // coarse-grid problems
922 if (!Doc_everything)
923 {
924 // Loop over all of the coarser levels
925 for (unsigned i = 1; i < Nlevel; i++)
926 {
927 // Delete the i-th coarse-grid MGProblem
928 delete Mg_hierarchy[i];
929
930 // Set it to be a null pointer
931 Mg_hierarchy[i] = 0;
932 }
933 }
934 // Otherwise, document everything!
935 else
936 {
937 // If the user wishes to document everything we run the self-test
938 self_test();
939 } // if (!Doc_everything)
940
941 // Indicate that the full setup has been completed
942 Has_been_setup = true;
943
944 // If we're allowed to output to the screen
945 if (!Suppress_all_output)
946 {
947 // Output the time taken to complete the full setup
948 double t_fs_end = TimingHelpers::timer();
949 double full_setup_time = t_fs_end - t_fs_start;
950
951 // Output the CPU time
952 oomph_info << "\nCPU time for full setup [sec]: " << full_setup_time
953 << std::endl;
954
955 // Notify user that the hierarchy of levels is complete
957 << "\n===============Multigrid Full Setup Complete=============="
958 << std::endl;
959 } // if (!Suppress_all_output)
960 } // End of full_setup
961
962 //===================================================================
963 /// Set up the MG hierarchy
964 //===================================================================
965 // Function to set up the hierachy of levels. Creates a vector of
966 // pointers to each MG level
967 template<unsigned DIM>
969 {
970 // Initialise the timer start variable
971 double t_m_start = 0.0;
972
973 // Notify the user if it is allowed
974 if (!Suppress_all_output)
975 {
976 // Notify user of progress
978 << "\n===============Creating Multigrid Hierarchy==============="
979 << std::endl;
980
981 // Start clock
982 t_m_start = TimingHelpers::timer();
983 }
984
985 // Create a bool to indicate whether or not we could create an unrefined
986 // copy. This bool will be assigned the value FALSE when the current copy
987 // is the last level of the multigrid hierarchy
988 bool managed_to_create_unrefined_copy = true;
989
990 // Now keep making copies and try to make an unrefined copy of
991 // the mesh
992 unsigned level = 0;
993
994 // Set up all of the levels by making a completely unrefined copy
995 // of the problem using the function make_new_problem
996 while (managed_to_create_unrefined_copy)
997 {
998 // Make a new object of the same type as the derived problem
999 MGProblem* new_problem_pt = Mg_problem_pt->make_new_problem();
1000
1001 // Do anything that needs to be done before we can refine the mesh
1002 new_problem_pt->actions_before_adapt();
1003
1004 // To create the next level in the hierarchy we need to create a mesh
1005 // which matches the refinement of the current problem and then unrefine
1006 // the mesh. This can alternatively be done using the function
1007 // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
1008 // reference mesh to do precisely the above with
1009 managed_to_create_unrefined_copy =
1010 new_problem_pt->mg_bulk_mesh_pt()
1012 Mg_hierarchy[level]->mg_bulk_mesh_pt());
1013
1014 // If we were able to unrefine the problem on the current level
1015 // then add the unrefined problem to a vector of the levels
1016 if (managed_to_create_unrefined_copy)
1017 {
1018 // Another level has been created so increment the level counter
1019 level++;
1020
1021 // If the documentation of everything has not been suppressed
1022 // then tell the user we managed to create another level
1023 if (!Suppress_all_output)
1024 {
1025 // Notify user that unrefinement was successful
1026 oomph_info << "\nSuccess! Level " << level << " has been created."
1027 << std::endl;
1028 }
1029
1030 // Do anything that needs to be done after refinement
1031 new_problem_pt->actions_after_adapt();
1032
1033 // Do the equation numbering for the new problem
1034 oomph_info << "\n - Number of equations: "
1035 << new_problem_pt->assign_eqn_numbers() << std::endl;
1036
1037 // Add the new problem pointer onto the vector of MG levels
1038 // and increment the value of level by 1
1039 Mg_hierarchy.push_back(new_problem_pt);
1040 }
1041 // If we weren't able to create an unrefined copy
1042 else
1043 {
1044 // Delete the new problem
1045 delete new_problem_pt;
1046
1047 // Make it a null pointer
1048 new_problem_pt = 0;
1049
1050 // Assign the number of levels to Nlevel
1051 Nlevel = Mg_hierarchy.size();
1052
1053 // If we're allowed to document then tell the user we've reached
1054 // the coarsest level of the hierarchy
1055 if (!Suppress_all_output)
1056 {
1057 // Notify the user
1058 oomph_info << "\n Reached the coarsest level! "
1059 << "Number of levels: " << Nlevel << std::endl;
1060 }
1061 } // if (managed_to_unrefine)
1062 } // while (managed_to_unrefine)
1063
1064 // Given that we know the number of levels in the hierarchy we can resize
1065 // the vectors which will store all the information required for our solver:
1066 // Resize the vector storing all of the system matrices
1067 Mg_matrices_storage_pt.resize(Nlevel, 0);
1068
1069 // Resize the vector storing all of the solution vectors (X_mg)
1070 X_mg_vectors_storage.resize(Nlevel);
1071
1072 // Resize the vector storing all of the RHS vectors (Rhs_mg)
1073 Rhs_mg_vectors_storage.resize(Nlevel);
1074
1075 // Resize the vector storing all of the residual vectors
1076 Residual_mg_vectors_storage.resize(Nlevel);
1077
1078 // Allocate space for the pre-smoother storage vector (remember, we do
1079 // not need a smoother on the coarsest level we use a direct solve there)
1080 Pre_smoothers_storage_pt.resize(Nlevel - 1, 0);
1081
1082 // Allocate space for the post-smoother storage vector (remember, we do
1083 // not need a smoother on the coarsest level we use a direct solve there)
1084 Post_smoothers_storage_pt.resize(Nlevel - 1, 0);
1085
1086 // Resize the vector storing all of the interpolation matrices
1087 Interpolation_matrices_storage_pt.resize(Nlevel - 1, 0);
1088
1089 // Resize the vector storing all of the restriction matrices
1090 Restriction_matrices_storage_pt.resize(Nlevel - 1, 0);
1091
1092 if (!Suppress_all_output)
1093 {
1094 // Stop clock
1095 double t_m_end = TimingHelpers::timer();
1096 double total_setup_time = double(t_m_end - t_m_start);
1098 << "\nCPU time for creation of hierarchy of MG problems [sec]: "
1099 << total_setup_time << std::endl;
1100
1101 // Notify user that the hierarchy of levels is complete
1103 << "\n===============Hierarchy Creation Complete================"
1104 << "\n"
1105 << std::endl;
1106 }
1107 } // End of setup_mg_hierarchy
1108
1109 //===================================================================
1110 /// Set up the transfer matrices. Both the pure injection and
1111 /// full weighting method have been implemented here but it is highly
1112 /// recommended that full weighting is used in general. In both
1113 /// methods the transpose of the transfer matrix is used to transfer
1114 /// a vector back
1115 //===================================================================
1116 template<unsigned DIM>
1118 {
1119 // Initialise the timer start variable
1120 double t_r_start = 0.0;
1121
1122 // Notify the user (if we're allowed)
1123 if (!Suppress_all_output)
1124 {
1125 // Notify user of progress
1126 oomph_info << "Creating the transfer matrices ";
1127
1128 // Start the clock!
1129 t_r_start = TimingHelpers::timer();
1130 }
1131
1132 // If we're allowed to output information
1133 if (!Suppress_all_output)
1134 {
1135 // Say what method we're using
1136 oomph_info << "using full weighting (recommended).\n" << std::endl;
1137 }
1138
1139 // Using full weighting so use setup_interpolation_matrices.
1140 // Note: There are two methods to choose from here, the ideal choice is
1141 // setup_interpolation_matrices() but that requires a refineable mesh base
1142 if (dynamic_cast<TreeBasedRefineableMeshBase*>(
1143 Mg_problem_pt->mg_bulk_mesh_pt()))
1144 {
1145 setup_interpolation_matrices();
1146 }
1147 // If the mesh is unstructured we have to use the locate_zeta function
1148 // to set up the interpolation matrices
1149 else
1150 {
1151 setup_interpolation_matrices_unstructured();
1152 }
1153
1154 // Loop over all levels that will be assigned a restriction matrix
1155 set_restriction_matrices_as_interpolation_transposes();
1156
1157 // If we're allowed
1158 if (!Suppress_all_output)
1159 {
1160 // Stop the clock
1161 double t_r_end = TimingHelpers::timer();
1162 double total_G_setup_time = double(t_r_end - t_r_start);
1163 oomph_info << "CPU time for transfer matrices setup [sec]: "
1164 << total_G_setup_time << std::endl;
1165
1166 // Notify user that the hierarchy of levels is complete
1168 << "\n============Transfer Matrices Setup Complete=============="
1169 << "\n"
1170 << std::endl;
1171 }
1172 } // End of setup_transfer_matrices function
1173
1174 //===================================================================
1175 /// Set up the MG hierarchy structures
1176 //===================================================================
1177 // Function to set up the hierachy of levels. Creates a vector of
1178 // pointers to each MG level
1179 template<unsigned DIM>
1181 {
1182 // Initialise the timer start variable
1183 double t_m_start = 0.0;
1184
1185 // Start the clock (if we're allowed to time things)
1186 if (!Suppress_all_output)
1187 {
1188 // Start the clock
1189 t_m_start = TimingHelpers::timer();
1190 }
1191
1192 // Allocate space for the system matrix on each level
1193 for (unsigned i = 0; i < Nlevel; i++)
1194 {
1195 // Dynamically allocate a new CRDoubleMatrix
1196 Mg_matrices_storage_pt[i] = new CRDoubleMatrix;
1197 }
1198
1199 // Loop over each level and extract the system matrix, solution vector
1200 // right-hand side vector and residual vector (to store the value of r=b-Ax)
1201 for (unsigned i = 0; i < Nlevel; i++)
1202 {
1203 // If we're allowed to output
1204 if (!Suppress_all_output)
1205 {
1206 // Output the level we're working on
1207 oomph_info << "Setting up MG structures on level: " << i << "\n"
1208 << std::endl;
1209 }
1210
1211 // Resize the solution and RHS vector
1212 unsigned n_dof = Mg_hierarchy[i]->ndof();
1214 Mg_hierarchy[i]->communicator_pt(), n_dof, false);
1215
1216#ifdef PARANOID
1217#ifdef OOMPH_HAS_MPI
1218 // Set up the warning messages
1219 std::string warning_message = "Setup of distribution has not been ";
1220 warning_message += "tested with MPI.";
1221
1222 // If we're not running the code in serial
1223 if (dist_pt->communicator_pt()->nproc() > 1)
1224 {
1225 // Throw a warning
1227 warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1228 }
1229#endif
1230#endif
1231
1232 // Build the approximate solution
1233 X_mg_vectors_storage[i].clear();
1234 X_mg_vectors_storage[i].build(dist_pt);
1235
1236 // Build the point source function
1237 Rhs_mg_vectors_storage[i].clear();
1238 Rhs_mg_vectors_storage[i].build(dist_pt);
1239
1240 // Build the residual vector (r=b-Ax)
1241 Residual_mg_vectors_storage[i].clear();
1242 Residual_mg_vectors_storage[i].build(dist_pt);
1243
1244 // Delete the distribution pointer
1245 delete dist_pt;
1246
1247 // Make it a null pointer
1248 dist_pt = 0;
1249
1250 // Build the matrix distribution
1251 Mg_matrices_storage_pt[i]->clear();
1252 Mg_matrices_storage_pt[i]->distribution_pt()->build(
1253 Mg_hierarchy[i]->communicator_pt(), n_dof, false);
1254
1255 // Compute system matrix on the current level. On the finest level of the
1256 // hierarchy the system matrix and RHS vector is given by the Jacobian and
1257 // vector of residuals which define the original problem which the
1258 // Galerkin approximation to the system matrix is used on the subsequent
1259 // levels so that the correct contributions are taken from each dof on the
1260 // level above (that is to say, it should match the contribution taken
1261 // from the solution vector and RHS vector on the level above)
1262 if (i == 0)
1263 {
1264 // Initialise the timer start variable
1265 double t_jac_start = 0.0;
1266
1267 // If we're allowed to output things
1268 if (!Suppress_all_output)
1269 {
1270 // Start timer for Jacobian setup
1271 t_jac_start = TimingHelpers::timer();
1272 }
1273
1274 // The system matrix on the finest level is the Jacobian and the RHS
1275 // vector is given by the residual vector which accompanies the Jacobian
1276 Mg_hierarchy[0]->get_jacobian(Rhs_mg_vectors_storage[0],
1277 *Mg_matrices_storage_pt[0]);
1278
1279 if (!Suppress_all_output)
1280 {
1281 // Document the time taken
1282 double t_jac_end = TimingHelpers::timer();
1283 double jacobian_setup_time = t_jac_end - t_jac_start;
1284 oomph_info << " - Time for setup of Jacobian [sec]: "
1285 << jacobian_setup_time << "\n"
1286 << std::endl;
1287 }
1288 }
1289 else
1290 {
1291 // Initialise the timer start variable
1292 double t_gal_start = 0.0;
1293
1294 // If we're allowed
1295 if (!Suppress_all_output)
1296 {
1297 // Start timer for Galerkin matrix calculation
1298 t_gal_start = TimingHelpers::timer();
1299 }
1300
1301 // The system matrix on the coarser levels must be formed using the
1302 // Galerkin approximation which we do by calculating the product
1303 // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
1304 // finer grid system matrix is formed by multiplying by the (fine grid)
1305 // restriction matrix from the left and the (fine grid) interpolation
1306 // matrix from the left
1307
1308 // First we need to calculate A^h * I^h_2h which we store as A^2h
1309 Mg_matrices_storage_pt[i - 1]->multiply(
1310 *Interpolation_matrices_storage_pt[i - 1],
1311 *Mg_matrices_storage_pt[i]);
1312
1313 // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1314 // was just calculated. This updates A^2h to give us the true
1315 // Galerkin approximation to the finer grid matrix
1316 Restriction_matrices_storage_pt[i - 1]->multiply(
1317 *Mg_matrices_storage_pt[i], *Mg_matrices_storage_pt[i]);
1318
1319 // If the user did not choose to suppress everything
1320 if (!Suppress_all_output)
1321 {
1322 // End timer for Galerkin matrix calculation
1323 double t_gal_end = TimingHelpers::timer();
1324
1325 // Calculate setup time
1326 double galerkin_matrix_calculation_time = t_gal_end - t_gal_start;
1327
1328 // Document the time taken
1330 << " - Time for system matrix formation using the Galerkin "
1331 << "approximation [sec]: " << galerkin_matrix_calculation_time
1332 << "\n"
1333 << std::endl;
1334 }
1335 } // if (i==0) else
1336 } // for (unsigned i=0;i<Nlevel;i++)
1337
1338 // If we're allowed to output
1339 if (!Suppress_all_output)
1340 {
1341 // Stop clock
1342 double t_m_end = TimingHelpers::timer();
1343 double total_setup_time = double(t_m_end - t_m_start);
1344 oomph_info << "Total CPU time for setup of MG structures [sec]: "
1345 << total_setup_time << std::endl;
1346
1347 // Notify user that the hierarchy of levels is complete
1348 oomph_info << "\n============"
1349 << "Multigrid Structures Setup Complete"
1350 << "===========\n"
1351 << std::endl;
1352 }
1353 } // End of setup_mg_structures
1354
1355
1356 //=========================================================================
1357 /// Set up the smoothers on all levels
1358 //=========================================================================
1359 template<unsigned DIM>
1361 {
1362 // Initialise the timer start variable
1363 double t_m_start = 0.0;
1364
1365 // Start the clock (if we're allowed to time things)
1366 if (!Suppress_all_output)
1367 {
1368 // Notify user
1369 oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
1370
1371 // Start the clock
1372 t_m_start = TimingHelpers::timer();
1373 }
1374
1375 // Loop over the levels and assign the pre- and post-smoother points
1376 for (unsigned i = 0; i < Nlevel - 1; i++)
1377 {
1378 // If the pre-smoother factory function pointer hasn't been assigned
1379 // then we simply create a new instance of the DampedJacobi smoother
1380 // which is the default pre-smoother
1381 if (0 == Pre_smoother_factory_function_pt)
1382 {
1383 Pre_smoothers_storage_pt[i] = new DampedJacobi<CRDoubleMatrix>;
1384 }
1385 // Otherwise we use the pre-smoother factory function pointer to
1386 // generate a new pre-smoother
1387 else
1388 {
1389 // Get a pointer to an object of the same type as the pre-smoother
1390 Pre_smoothers_storage_pt[i] = (*Pre_smoother_factory_function_pt)();
1391 }
1392
1393 // If the post-smoother factory function pointer hasn't been assigned
1394 // then we simply create a new instance of the DampedJacobi smoother
1395 // which is the default post-smoother
1396 if (0 == Post_smoother_factory_function_pt)
1397 {
1398 Post_smoothers_storage_pt[i] = new DampedJacobi<CRDoubleMatrix>;
1399 }
1400 // Otherwise we use the post-smoother factory function pointer to
1401 // generate a new post-smoother
1402 else
1403 {
1404 // Get a pointer to an object of the same type as the post-smoother
1405 Post_smoothers_storage_pt[i] = (*Post_smoother_factory_function_pt)();
1406 }
1407 }
1408
1409 // Set the tolerance for the pre- and post-smoothers. The norm of the
1410 // solution which is compared against the tolerance is calculated
1411 // differently to the multigrid solver. To ensure that the smoother
1412 // continues to solve for the prescribed number of iterations we
1413 // lower the tolerance
1414 for (unsigned i = 0; i < Nlevel - 1; i++)
1415 {
1416 // Set the tolerance of the i-th level pre-smoother
1417 Pre_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
1418
1419 // Set the tolerance of the i-th level post-smoother
1420 Post_smoothers_storage_pt[i]->tolerance() = 1.0e-16;
1421 }
1422
1423 // Set the number of pre- and post-smoothing iterations in each
1424 // pre- and post-smoother
1425 for (unsigned i = 0; i < Nlevel - 1; i++)
1426 {
1427 // Set the number of pre-smoothing iterations as the value of Npre_smooth
1428 Pre_smoothers_storage_pt[i]->max_iter() = Npre_smooth;
1429
1430 // Set the number of pre-smoothing iterations as the value of Npost_smooth
1431 Post_smoothers_storage_pt[i]->max_iter() = Npost_smooth;
1432 }
1433
1434 // Complete the setup of all of the smoothers
1435 for (unsigned i = 0; i < Nlevel - 1; i++)
1436 {
1437 // Pass a pointer to the system matrix on the i-th level to the i-th
1438 // level pre-smoother
1439 Pre_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
1440
1441 // Pass a pointer to the system matrix on the i-th level to the i-th
1442 // level post-smoother
1443 Post_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
1444 }
1445
1446 // Set up the distributions of each smoother
1447 for (unsigned i = 0; i < Nlevel - 1; i++)
1448 {
1449 // Get the number of dofs on the i-th level of the hierarchy.
1450 // This will be the same as the size of the solution vector
1451 // associated with the i-th level
1452 unsigned n_dof = X_mg_vectors_storage[i].nrow();
1453
1454 // Create a LinearAlgebraDistribution which will be passed to the
1455 // linear solver
1457 Mg_hierarchy[i]->communicator_pt(), n_dof, false);
1458
1459#ifdef PARANOID
1460#ifdef OOMPH_HAS_MPI
1461 // Set up the warning messages
1462 std::string warning_message =
1463 "Setup of pre- and post-smoother distribution ";
1464 warning_message += "has not been tested with MPI.";
1465
1466 // If we're not running the code in serial
1467 if (dist.communicator_pt()->nproc() > 1)
1468 {
1469 // Throw a warning
1471 warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1472 }
1473#endif
1474#endif
1475
1476 // Build the distribution of the pre-smoother
1477 Pre_smoothers_storage_pt[i]->build_distribution(dist);
1478
1479 // Build the distribution of the post-smoother
1480 Post_smoothers_storage_pt[i]->build_distribution(dist);
1481 }
1482
1483 // Disable the smoother output on this level
1484 if (!Doc_time)
1485 {
1486 disable_smoother_and_superlu_doc_time();
1487 }
1488
1489 // If we're allowed to output
1490 if (!Suppress_all_output)
1491 {
1492 // Stop clock
1493 double t_m_end = TimingHelpers::timer();
1494 double total_setup_time = double(t_m_end - t_m_start);
1495 oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
1496 << total_setup_time << std::endl;
1497
1498 // Notify user that the extraction is complete
1500 << "\n==================Smoother Setup Complete================="
1501 << std::endl;
1502 }
1503 } // End of setup_smoothers
1504
1505
1506 //===================================================================
1507 /// Setup the interpolation matrices
1508 //===================================================================
1509 template<unsigned DIM>
1511 {
1512 // Variable to hold the number of sons in an element
1513 unsigned n_sons;
1514
1515 // Number of son elements
1516 if (DIM == 2)
1517 {
1518 n_sons = 4;
1519 }
1520 else if (DIM == 3)
1521 {
1522 n_sons = 8;
1523 }
1524 else
1525 {
1526 std::ostringstream error_message_stream;
1527 error_message_stream << "DIM should be 2 or 3 not " << DIM << std::endl;
1528 throw OomphLibError(error_message_stream.str(),
1529 OOMPH_CURRENT_FUNCTION,
1530 OOMPH_EXCEPTION_LOCATION);
1531 }
1532
1533 // Vector of local coordinates in the element
1534 Vector<double> s(DIM, 0.0);
1535
1536 // Loop over each level (apart from the coarsest level; an interpolation
1537 // matrix and thus a restriction matrix is not needed here), with 0 being
1538 // the finest level and Nlevel-1 being the coarsest level
1539 for (unsigned level = 0; level < Nlevel - 1; level++)
1540 {
1541 // Assign values to a couple of variables to help out
1542 unsigned coarse_level = level + 1;
1543 unsigned fine_level = level;
1544
1545 // Make a pointer to the mesh on the finer level and dynamic_cast
1546 // it as an object of the refineable mesh class
1547 RefineableMeshBase* ref_fine_mesh_pt = dynamic_cast<RefineableMeshBase*>(
1548 Mg_hierarchy[fine_level]->mg_bulk_mesh_pt());
1549
1550 // Make a pointer to the mesh on the coarse level and dynamic_cast
1551 // it as an object of the refineable mesh class
1552 RefineableMeshBase* ref_coarse_mesh_pt =
1553 dynamic_cast<RefineableMeshBase*>(
1554 Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
1555
1556 // Access information about the number of elements in the fine mesh
1557 // from the pointer to the fine mesh (to loop over the rows of the
1558 // interpolation matrix)
1559 unsigned fine_n_element = ref_fine_mesh_pt->nelement();
1560
1561 // The numbers of rows and columns in the interpolation matrix. The
1562 // number of unknowns has been divided by 2 since there are 2 dofs at
1563 // each node in the mesh corresponding to the real and imaginary part
1564 // of the solution
1565 unsigned n_rows = Mg_hierarchy[fine_level]->ndof();
1566 unsigned n_cols = Mg_hierarchy[coarse_level]->ndof();
1567
1568 // Mapping relating the pointers to related elements in the coarse and
1569 // fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]
1570 // If the element in the fine mesh has been unrefined between these two
1571 // levels, this map returns the father element in the coarsened mesh.
1572 // If this element in the fine mesh has not been unrefined,
1573 // the map returns the pointer to the same-sized equivalent
1574 // element in the coarsened mesh.
1575 std::map<RefineableQElement<DIM>*, RefineableQElement<DIM>*>
1576 coarse_mesh_reference_element_pt;
1577
1578 // Counter of elements in coarse and fine meshes: Start with element
1579 // zero in both meshes.
1580 unsigned e_coarse = 0;
1581 unsigned e_fine = 0;
1582
1583 // While loop over fine elements (while because we're
1584 // incrementing the counter internally if the element was
1585 // unrefined...)
1586 while (e_fine < fine_n_element)
1587 {
1588 // Pointer to element in fine mesh
1589 RefineableQElement<DIM>* el_fine_pt =
1590 dynamic_cast<RefineableQElement<DIM>*>(
1591 ref_fine_mesh_pt->finite_element_pt(e_fine));
1592
1593 // Pointer to element in coarse mesh
1594 RefineableQElement<DIM>* el_coarse_pt =
1595 dynamic_cast<RefineableQElement<DIM>*>(
1596 ref_coarse_mesh_pt->finite_element_pt(e_coarse));
1597
1598 // If the levels are different then the element in the fine
1599 // mesh has been unrefined between these two levels
1600 if (el_fine_pt->tree_pt()->level() != el_coarse_pt->tree_pt()->level())
1601 {
1602 // The element in the fine mesh has been unrefined between these two
1603 // levels. Hence it and its three brothers (ASSUMED to be stored
1604 // consecutively in the Mesh's vector of pointers to its constituent
1605 // elements -- we'll check this!) share the same father element in
1606 // the coarse mesh, currently pointed to by el_coarse_pt.
1607 for (unsigned i = 0; i < n_sons; i++)
1608 {
1609 // Set mapping to father element in coarse mesh
1610 coarse_mesh_reference_element_pt
1611 [dynamic_cast<RefineableQElement<DIM>*>(
1612 ref_fine_mesh_pt->finite_element_pt(e_fine))] = el_coarse_pt;
1613
1614 // Increment counter for elements in fine mesh
1615 e_fine++;
1616 }
1617 }
1618 // The element in the fine mesh has not been unrefined between
1619 // these two levels, so the reference element is the same-sized
1620 // equivalent element in the coarse mesh
1621 else
1622 {
1623 // Set the mapping between the two elements since they are
1624 // the same element
1625 coarse_mesh_reference_element_pt[el_fine_pt] = el_coarse_pt;
1626
1627 // Increment counter for elements in fine mesh
1628 e_fine++;
1629 }
1630 // Increment counter for elements in coarse mesh
1631 e_coarse++;
1632 } // End of while loop for setting up element-coincidences
1633
1634 // To allow update of a row only once we use stl vectors for bools
1635 std::vector<bool> contribution_made(n_rows, false);
1636
1637 // Make storage vectors to form the interpolation matrix using a
1638 // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
1639 // row_start contains entries which tells us at which entry of the
1640 // vector column_start we start on the i-th row of the actual matrix.
1641 // The entries of column_start indicate the column position of each
1642 // non-zero entry in every row. This runs through the first row
1643 // from left to right then the second row (again, left to right)
1644 // and so on until the end. The entries in value are the entries in
1645 // the interpolation matrix. The vector column_start has the same length
1646 // as the vector value.
1647 Vector<double> value;
1648 Vector<int> column_index;
1649 Vector<int> row_start(n_rows + 1);
1650
1651 // The value of index will tell us which row of the interpolation matrix
1652 // we're working on in the following for loop
1653 unsigned index = 0;
1654
1655 // New loop to go over each element in the fine mesh
1656 for (unsigned k = 0; k < fine_n_element; k++)
1657 {
1658 // Set a pointer to the element in the fine mesh
1659 RefineableQElement<DIM>* el_fine_pt =
1660 dynamic_cast<RefineableQElement<DIM>*>(
1661 ref_fine_mesh_pt->finite_element_pt(k));
1662
1663 // Get the reference element (either the father element or the
1664 // same-sized element) in the coarse mesh
1665 RefineableQElement<DIM>* el_coarse_pt =
1666 coarse_mesh_reference_element_pt[el_fine_pt];
1667
1668 // Find out what type of son it is (set to OMEGA if no unrefinement
1669 // took place)
1670 int son_type = 0;
1671
1672 // Check if the elements are different on both levels (i.e. to check
1673 // if any unrefinement took place)
1674 if (el_fine_pt->tree_pt()->level() != el_coarse_pt->tree_pt()->level())
1675 {
1676 // If there was refinement we need to find the son type
1677 son_type = el_fine_pt->tree_pt()->son_type();
1678 }
1679 else
1680 {
1681 // If there was no refinement then the son_type is given by the
1682 // value of Tree::OMEGA
1683 son_type = Tree::OMEGA;
1684 }
1685
1686 // Find the number of nodes in the fine mesh element
1687 unsigned nnod_fine = el_fine_pt->nnode();
1688
1689 // Loop through all the nodes in an element in the fine mesh
1690 for (unsigned i = 0; i < nnod_fine; i++)
1691 {
1692 // Row number in interpolation matrix: Global equation number
1693 // of the d.o.f. stored at this node in the fine element
1694 int ii = el_fine_pt->node_pt(i)->eqn_number(0);
1695
1696 // Check whether or not the node is a proper d.o.f.
1697 if (ii >= 0)
1698 {
1699 // Only assign values to the given row of the interpolation
1700 // matrix if they haven't already been assigned
1701 if (contribution_made[ii] == false)
1702 {
1703 // The value of index was initialised when we allocated space
1704 // for the three vectors to store the CRDoubleMatrix. We use
1705 // index to go through the entries of the row_start vector.
1706 // The next row starts at the value.size() (draw out the entries
1707 // in value if this doesn't make sense noting that the storage
1708 // for the vector 'value' is dynamically allocated)
1709 row_start[index] = value.size();
1710
1711 // Calculate the local coordinates of the given node
1712 el_fine_pt->local_coordinate_of_node(i, s);
1713
1714 // Find the local coordinates s, of the present node, in the
1715 // reference element in the coarse mesh, given the element's
1716 // son_type (e.g. SW,SE... )
1717 level_up_local_coord_of_node(son_type, s);
1718
1719 // Allocate space for shape functions in the coarse mesh
1720 Shape psi(el_coarse_pt->nnode());
1721
1722 // Set the shape function in the reference element
1723 el_coarse_pt->shape(s, psi);
1724
1725 // Auxiliary storage
1726 std::map<unsigned, double> contribution;
1727 Vector<unsigned> keys;
1728
1729 // Find the number of nodes in an element in the coarse mesh
1730 unsigned nnod_coarse = el_coarse_pt->nnode();
1731
1732 // Loop through all the nodes in the reference element
1733 for (unsigned j = 0; j < nnod_coarse; j++)
1734 {
1735 // Column number in interpolation matrix: Global equation
1736 // number of the d.o.f. stored at this node in the coarse
1737 // element
1738 int jj = el_coarse_pt->node_pt(j)->eqn_number(0);
1739
1740 // If the value stored at this node is pinned or hanging
1741 if (jj < 0)
1742 {
1743 // Hanging node: In this case we need to accumulate the
1744 // contributions from the master nodes
1745 if (el_coarse_pt->node_pt(j)->is_hanging())
1746 {
1747 // Find the number of master nodes of the hanging
1748 // the node in the reference element
1749 HangInfo* hang_info_pt =
1750 el_coarse_pt->node_pt(j)->hanging_pt();
1751 unsigned nmaster = hang_info_pt->nmaster();
1752
1753 // Loop over the master nodes
1754 for (unsigned i_master = 0; i_master < nmaster; i_master++)
1755 {
1756 // Set up a pointer to the master node
1757 Node* master_node_pt =
1758 hang_info_pt->master_node_pt(i_master);
1759
1760 // The column number in the interpolation matrix: the
1761 // global equation number of the d.o.f. stored at this
1762 // master node for the coarse element
1763 int master_jj = master_node_pt->eqn_number(0);
1764
1765 // Is the master node a proper d.o.f.?
1766 if (master_jj >= 0)
1767 {
1768 // If the weight of the master node is non-zero
1769 if (psi(j) * hang_info_pt->master_weight(i_master) !=
1770 0.0)
1771 {
1772 contribution[master_jj] +=
1773 psi(j) * hang_info_pt->master_weight(i_master);
1774 }
1775 } // if (master_jj>=0)
1776 } // for (unsigned i_master=0;i_master<nmaster;i_master++)
1777 } // if (el_coarse_pt->node_pt(j)->is_hanging())
1778 }
1779 // In the case that the node is not pinned or hanging
1780 else
1781 {
1782 // If we can get a nonzero contribution from the shape
1783 // function
1784 if (psi(j) != 0.0)
1785 {
1786 contribution[jj] += psi(j);
1787 }
1788 } // if (jj<0) else
1789 } // for (unsigned j=0;j<nnod_coarse;j++)
1790
1791 // Put the contributions into the value vector
1792 for (std::map<unsigned, double>::iterator it =
1793 contribution.begin();
1794 it != contribution.end();
1795 ++it)
1796 {
1797 if (it->second != 0)
1798 {
1799 column_index.push_back(it->first);
1800 value.push_back(it->second);
1801 }
1802 } // for (std::map<unsigned,double>::iterator it=...)
1803
1804 // Increment the value of index by 1 to indicate that we have
1805 // finished the row, or equivalently, covered another global
1806 // node in the fine mesh
1807 index++;
1808
1809 // Change the entry in contribution_made to true now to indicate
1810 // that the row has been filled
1811 contribution_made[ii] = true;
1812 } // if(contribution_made[ii]==false)
1813 } // if (ii>=0)
1814 } // for(unsigned i=0;i<nnod_element;i++)
1815 } // for (unsigned k=0;k<fine_n_element;k++)
1816
1817 // Set the last entry in the row_start vector
1818 row_start[n_rows] = value.size();
1819
1820 // Set the interpolation matrix to be that formed as the CRDoubleMatrix
1821 // using the vectors value, row_start, column_index and the value
1822 // of fine_n_unknowns and coarse_n_unknowns
1823 interpolation_matrix_set(
1824 level, value, column_index, row_start, n_cols, n_rows);
1825 } // for (unsigned level=0;level<Nlevel-1;level++)
1826 } // End of setup_interpolation_matrices
1827
1828 //===================================================================
1829 /// Setup the interpolation matrices
1830 //===================================================================
1831 template<unsigned DIM>
1833 {
1834 // Vector of local coordinates in the element
1835 Vector<double> s(DIM, 0.0);
1836
1837 // Loop over each level (apart from the coarsest level; an interpolation
1838 // matrix and thus a restriction matrix is not needed here), with 0 being
1839 // the finest level and Nlevel-1 being the coarsest level
1840 for (unsigned level = 0; level < Nlevel - 1; level++)
1841 {
1842 // Assign values to a couple of variables to help out
1843 unsigned coarse_level = level + 1;
1844 unsigned fine_level = level;
1845
1846 // Make a pointer to the mesh on the finer level and dynamic_cast
1847 // it as an object of the refineable mesh class
1848 Mesh* ref_fine_mesh_pt = Mg_hierarchy[fine_level]->mg_bulk_mesh_pt();
1849
1850 // To use the locate zeta functionality the coarse mesh must be of the
1851 // type MeshAsGeomObject
1852 MeshAsGeomObject* coarse_mesh_from_obj_pt =
1853 new MeshAsGeomObject(Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
1854
1855 // Access information about the number of degrees of freedom
1856 // from the pointers to the problem on each level
1857 unsigned coarse_n_unknowns = Mg_hierarchy[coarse_level]->ndof();
1858 unsigned fine_n_unknowns = Mg_hierarchy[fine_level]->ndof();
1859
1860 // Make storage vectors to form the interpolation matrix using a
1861 // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
1862 // row_start contains entries which tells us at which entry of the
1863 // vector column_start we start on the i-th row of the actual matrix.
1864 // The entries of column_start indicate the column position of each
1865 // non-zero entry in every row. This runs through the first row
1866 // from left to right then the second row (again, left to right)
1867 // and so on until the end. The entries in value are the entries in
1868 // the interpolation matrix. The vector column_start has the same length
1869 // as the vector value.
1870 Vector<double> value;
1871 Vector<int> column_index;
1872 Vector<int> row_start(fine_n_unknowns + 1);
1873
1874 // Vector to contain the (Eulerian) spatial location of the fine node
1875 Vector<double> fine_node_position(DIM);
1876
1877 // Find the number of nodes in the mesh
1878 unsigned n_node_fine_mesh = ref_fine_mesh_pt->nnode();
1879
1880 // Loop over the unknowns in the mesh
1881 for (unsigned i_fine_node = 0; i_fine_node < n_node_fine_mesh;
1882 i_fine_node++)
1883 {
1884 // Set a pointer to the i_fine_unknown-th node in the fine mesh
1885 Node* fine_node_pt = ref_fine_mesh_pt->node_pt(i_fine_node);
1886
1887 // Get the global equation number
1888 int i_fine = fine_node_pt->eqn_number(0);
1889
1890 // If the node is a proper d.o.f.
1891 if (i_fine >= 0)
1892 {
1893 // Row number in interpolation matrix: Global equation number
1894 // of the d.o.f. stored at this node in the fine element
1895 row_start[i_fine] = value.size();
1896
1897 // Get the (Eulerian) spatial location of the fine node
1898 fine_node_pt->position(fine_node_position);
1899
1900 // Create a null pointer to the GeomObject class
1901 GeomObject* el_pt = 0;
1902
1903 // Get the reference element (either the father element or the
1904 // same-sized element) in the coarse mesh using locate_zeta
1905 coarse_mesh_from_obj_pt->locate_zeta(fine_node_position, el_pt, s);
1906
1907 // Upcast GeomElement as a FiniteElement
1908 FiniteElement* el_coarse_pt = dynamic_cast<FiniteElement*>(el_pt);
1909
1910 // Find the number of nodes in the element
1911 unsigned n_node = el_coarse_pt->nnode();
1912
1913 // Allocate space for shape functions in the coarse mesh
1914 Shape psi(n_node);
1915
1916 // Calculate the geometric shape functions at local coordinate s
1917 el_coarse_pt->shape(s, psi);
1918
1919 // Auxiliary storage
1920 std::map<unsigned, double> contribution;
1921 Vector<unsigned> keys;
1922
1923 // Loop through all the nodes in the (coarse mesh) element containing
1924 // the node pointed to by fine_node_pt (fine mesh)
1925 for (unsigned j_node = 0; j_node < n_node; j_node++)
1926 {
1927 // Get the j_coarse_unknown-th node in the coarse element
1928 Node* coarse_node_pt = el_coarse_pt->node_pt(j_node);
1929
1930 // Column number in interpolation matrix: Global equation number of
1931 // the d.o.f. stored at this node in the coarse element
1932 int j_coarse = coarse_node_pt->eqn_number(0);
1933
1934 // If the value stored at this node is pinned or hanging
1935 if (j_coarse < 0)
1936 {
1937 // Hanging node: In this case we need to accumulate the
1938 // contributions from the master nodes
1939 if (el_coarse_pt->node_pt(j_node)->is_hanging())
1940 {
1941 // Find the number of master nodes of the hanging
1942 // the node in the reference element
1943 HangInfo* hang_info_pt = coarse_node_pt->hanging_pt();
1944 unsigned nmaster = hang_info_pt->nmaster();
1945
1946 // Loop over the master nodes
1947 for (unsigned i_master = 0; i_master < nmaster; i_master++)
1948 {
1949 // Set up a pointer to the master node
1950 Node* master_node_pt = hang_info_pt->master_node_pt(i_master);
1951
1952 // The column number in the interpolation matrix: the
1953 // global equation number of the d.o.f. stored at this master
1954 // node for the coarse element
1955 int master_jj = master_node_pt->eqn_number(0);
1956
1957 // Is the master node a proper d.o.f.?
1958 if (master_jj >= 0)
1959 {
1960 // If the weight of the master node is non-zero
1961 if (psi(j_node) * hang_info_pt->master_weight(i_master) !=
1962 0.0)
1963 {
1964 contribution[master_jj] +=
1965 psi(j_node) * hang_info_pt->master_weight(i_master);
1966 }
1967 } // End of if statement (check it's not a boundary node)
1968 } // End of the loop over the master nodes
1969 } // End of the if statement for only hanging nodes
1970 } // End of the if statement for pinned or hanging nodes
1971 // In the case that the node is not pinned or hanging
1972 else
1973 {
1974 // If we can get a nonzero contribution from the shape function
1975 // at the j_node-th node in the element
1976 if (psi(j_node) != 0.0)
1977 {
1978 contribution[j_coarse] += psi(j_node);
1979 }
1980 } // End of the if-else statement (check if the node was
1981 // pinned/hanging)
1982 } // Finished loop over the nodes j in the reference element (coarse)
1983
1984 // Put the contributions into the value vector
1985 for (std::map<unsigned, double>::iterator it = contribution.begin();
1986 it != contribution.end();
1987 ++it)
1988 {
1989 if (it->second != 0)
1990 {
1991 value.push_back(it->second);
1992 column_index.push_back(it->first);
1993 }
1994 } // End of putting contributions into the value vector
1995 } // End check (whether or not the fine node was a d.o.f.)
1996 } // End of the for-loop over nodes in the fine mesh
1997
1998 // Set the last entry of row_start
1999 row_start[fine_n_unknowns] = value.size();
2000
2001 // Set the interpolation matrix to be that formed as the CRDoubleMatrix
2002 // using the vectors value, row_start, column_index and the value
2003 // of fine_n_unknowns and coarse_n_unknowns
2004 interpolation_matrix_set(level,
2005 value,
2006 column_index,
2007 row_start,
2008 coarse_n_unknowns,
2009 fine_n_unknowns);
2010 } // End of loop over each level
2011 } // End of setup_interpolation_matrices_unstructured
2012
2013 //===================================================================
2014 /// Restrict residual (computed on current MG level) to
2015 /// next coarser mesh and stick it into the coarse mesh RHS vector
2016 /// using the restriction matrix (if restrict_flag=1) or the transpose
2017 /// of the interpolation matrix (if restrict_flag=2)
2018 //===================================================================
2019 template<unsigned DIM>
2020 void MGSolver<DIM>::restrict_residual(const unsigned& level)
2021 {
2022#ifdef PARANOID
2023 // Check to make sure we can actually restrict the vector
2024 if (!(level < Nlevel - 1))
2025 {
2026 throw OomphLibError("Input exceeds the possible parameter choice.",
2027 OOMPH_CURRENT_FUNCTION,
2028 OOMPH_EXCEPTION_LOCATION);
2029 }
2030#endif
2031
2032 // Multiply the residual vector by the restriction matrix on the level-th
2033 // level (to restrict the vector down to the next coarser level)
2034 Restriction_matrices_storage_pt[level]->multiply(
2035 Residual_mg_vectors_storage[level], Rhs_mg_vectors_storage[level + 1]);
2036 } // End of restrict_residual
2037
2038 //===================================================================
2039 /// Interpolate solution at current level onto
2040 /// next finer mesh and correct the solution x at that level
2041 //===================================================================
2042 template<unsigned DIM>
2043 void MGSolver<DIM>::interpolate_and_correct(const unsigned& level)
2044 {
2045#ifdef PARANOID
2046 // Check to make sure we can actually restrict the vector
2047 if (!(level > 0))
2048 {
2049 throw OomphLibError("Input level exceeds the possible parameter choice.",
2050 OOMPH_CURRENT_FUNCTION,
2051 OOMPH_EXCEPTION_LOCATION);
2052 }
2053#endif
2054
2055 // Build distribution of a temporary vector
2056 DoubleVector temp_soln(X_mg_vectors_storage[level - 1].distribution_pt());
2057
2058 // Interpolate the solution vector
2059 Interpolation_matrices_storage_pt[level - 1]->multiply(
2060 X_mg_vectors_storage[level], temp_soln);
2061
2062 // Update
2063 X_mg_vectors_storage[level - 1] += temp_soln;
2064 } // End of interpolate_and_correct
2065
2066 //===================================================================
2067 /// Modify the restriction matrices
2068 //===================================================================
2069 template<unsigned DIM>
2071 {
2072 // Create a null pointer
2073 CRDoubleMatrix* restriction_matrix_pt = 0;
2074
2075 // Loop over the levels
2076 for (unsigned level = 0; level < Nlevel - 1; level++)
2077 {
2078 // Store a pointer to the (level)-th restriction matrix
2079 restriction_matrix_pt = Restriction_matrices_storage_pt[level];
2080
2081 // Get access to the row start data
2082 const int* row_start_pt = restriction_matrix_pt->row_start();
2083
2084 // Get access to the matrix entries
2085 double* value_pt = restriction_matrix_pt->value();
2086
2087 // Initialise an auxiliary variable to store the index of the start
2088 // of the i-th row
2089 unsigned start_index = 0;
2090
2091 // Initialise an auxiliary variable to store the index of the start
2092 // of the (i+1)-th row
2093 unsigned end_index = 0;
2094
2095 // Store the number of rows in the matrix
2096 unsigned n_row = restriction_matrix_pt->nrow();
2097
2098 // Loop over the rows of the matrix
2099 for (unsigned i = 0; i < n_row; i++)
2100 {
2101 // Index associated with the start of the i-th row
2102 start_index = row_start_pt[i];
2103
2104 // Index associated with the start of the (i+1)-th row
2105 end_index = row_start_pt[i + 1];
2106
2107 // Variable to store the sum of the absolute values of the off-diagonal
2108 // entries in the i-th row
2109 double row_sum = 0.0;
2110
2111 // Loop over the entries in the row
2112 for (unsigned j = start_index; j < end_index; j++)
2113 {
2114 // Add the absolute value of this entry to the off-diagonal sum
2115 row_sum += value_pt[j];
2116 }
2117
2118 // Loop over the entries in the row
2119 for (unsigned j = start_index; j < end_index; j++)
2120 {
2121 // Normalise the row entry
2122 value_pt[j] /= row_sum;
2123 }
2124 } // for (unsigned i=0;i<n_row;i++)
2125 } // for (unsigned level=0;level<Nlevel-1;level++)
2126 } // End of modify_restriction_matrices
2127
2128 //===================================================================
2129 /// Makes a vector, restricts it down the levels of the hierarchy
2130 /// and documents it at each level. After this is done the vector is
2131 /// interpolated up the levels of the hierarchy with the output
2132 /// being documented at each level
2133 //===================================================================
2134 template<unsigned DIM>
2136 {
2137 // Start the timer
2138 double t_st_start = TimingHelpers::timer();
2139
2140 // Notify user that the hierarchy of levels is complete
2141 oomph_info << "\nStarting the multigrid solver self-test." << std::endl;
2142
2143 // Resize the vector storing all of the restricted vectors
2144 Restriction_self_test_vectors_storage.resize(Nlevel);
2145
2146 // Resize the vector storing all of the interpolated vectors
2147 Interpolation_self_test_vectors_storage.resize(Nlevel);
2148
2149 // Find the number of dofs on the finest level
2150 unsigned n_dof = X_mg_vectors_storage[0].nrow();
2151
2152 // Need to set up the distribution of the finest level vector
2154 Mg_problem_pt->communicator_pt(), n_dof, false);
2155
2156#ifdef PARANOID
2157#ifdef OOMPH_HAS_MPI
2158 // Set up the warning messages
2159 std::string warning_message = "Setup of distribution has not been ";
2160 warning_message += "tested with MPI.";
2161
2162 // If we're not running the code in serial
2163 if (dist_pt->communicator_pt()->nproc() > 1)
2164 {
2165 // Throw a warning
2167 warning_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2168 }
2169#endif
2170#endif
2171
2172 // Build the vector using the distribution of the restriction matrix
2173 Restriction_self_test_vectors_storage[0].build(dist_pt);
2174
2175 // Delete the distribution pointer
2176 delete dist_pt;
2177
2178 // Make it a null pointer
2179 dist_pt = 0;
2180
2181 // Now assign the values to the first vector. This will be restricted down
2182 // the levels of the hierarchy then back up
2183 set_self_test_vector();
2184
2185 // Call the restriction self-test
2186 restriction_self_test();
2187
2188 // Set the coarest level vector in Restriction_self_test_vectors_storage
2189 // to be the last entry in Interpolation_self_test_vectors_storage. This
2190 // will be interpolated up to the finest level in interpolation_self_test()
2191 Interpolation_self_test_vectors_storage[Nlevel - 1] =
2192 Restriction_self_test_vectors_storage[Nlevel - 1];
2193
2194 // Call the interpolation self-test
2195 interpolation_self_test();
2196
2197 // Destroy all of the created restriction self-test vectors
2198 Restriction_self_test_vectors_storage.resize(0);
2199
2200 // Destroy all of the created interpolation self-test vectors
2201 Interpolation_self_test_vectors_storage.resize(0);
2202
2203 // Stop the clock
2204 double t_st_end = TimingHelpers::timer();
2205 double total_st_time = double(t_st_end - t_st_start);
2206 oomph_info << "\nCPU time for self-test [sec]: " << total_st_time
2207 << std::endl;
2208
2209 // Notify user that the hierarchy of levels is complete
2210 oomph_info << "\n====================Self-Test Complete===================="
2211 << std::endl;
2212 } // End of self_test
2213
2214 //===================================================================
2215 /// Sets the initial vector to be used in the restriction and
2216 /// interpolation self-tests
2217 //===================================================================
2218 template<unsigned DIM>
2220 {
2221 // Set up a pointer to the refineable mesh
2222 TreeBasedRefineableMeshBase* bulk_mesh_pt =
2223 Mg_problem_pt->mg_bulk_mesh_pt();
2224
2225 // Find the number of elements in the bulk mesh
2226 unsigned n_el = bulk_mesh_pt->nelement();
2227
2228 // Get the dimension of the problem (assumed to be the dimension of any
2229 // node in the mesh)
2230 unsigned n_dim = bulk_mesh_pt->node_pt(0)->ndim();
2231
2232 // Find the number of nodes in an element
2233 unsigned nnod = bulk_mesh_pt->finite_element_pt(0)->nnode();
2234
2235 // Loop over all elements
2236 for (unsigned e = 0; e < n_el; e++)
2237 {
2238 // Get pointer to element
2239 FiniteElement* el_pt = bulk_mesh_pt->finite_element_pt(e);
2240
2241 // Loop over nodes
2242 for (unsigned j = 0; j < nnod; j++)
2243 {
2244 // Pointer to node
2245 Node* nod_pt = el_pt->node_pt(j);
2246
2247 // Sanity check
2248 if (nod_pt->nvalue() != 1)
2249 {
2250 // Set up an output stream
2251 std::ostringstream error_stream;
2252
2253 // Output the error message
2254 error_stream << "Sorry, not sure what to do here! I can't deal with "
2255 << nod_pt->nvalue() << " values!" << std::endl;
2256
2257 // Throw an error to indicate that it was not possible to plot the
2258 // data
2259 throw OomphLibError(error_stream.str(),
2260 OOMPH_CURRENT_FUNCTION,
2261 OOMPH_EXCEPTION_LOCATION);
2262 }
2263
2264 // Free or pinned?
2265 int eqn_num = nod_pt->eqn_number(0);
2266
2267 // If it's a free node
2268 if (eqn_num >= 0)
2269 {
2270 // Initialise the variable coordinate_sum
2271 double coordinate_sum = 0.0;
2272
2273 // Loop over the coordinates of the node
2274 for (unsigned i = 0; i < n_dim; i++)
2275 {
2276 // Increment coordinate_sum by the value of x(i)
2277 coordinate_sum += nod_pt->x(i);
2278 }
2279
2280 // Store the value of pi in cache
2281 double pi = MathematicalConstants::Pi;
2282
2283 // Make the vector represent a constant function
2284 Restriction_self_test_vectors_storage[0][eqn_num] =
2285 sin(pi * (nod_pt->x(0))) * sin(pi * (nod_pt->x(1)));
2286 }
2287 } // for (unsigned j=0;j<nnod;j++)
2288 } // for (unsigned e=0;e<n_el;e++)
2289
2290 // // Get the size of the vector
2291 // unsigned n_row=Restriction_self_test_vectors_storage[0].nrow();
2292
2293 // // Loop over the entries of the vector
2294 // for (unsigned i=0;i<n_row;i++)
2295 // {
2296 // // Make the vector represent a constant function
2297 // Restriction_self_test_vectors_storage[0][i]=1.0;
2298 // }
2299 } // End of set_self_test_vector
2300
2301 //===================================================================
2302 /// Function which implements a self-test to restrict the
2303 /// residual vector down all of the coarser grids and output
2304 /// the restricted vectors to file
2305 //===================================================================
2306 template<unsigned DIM>
2308 {
2309 // Set the start of the output file name. The full names will
2310 // set after we calculate the restricted vectors
2311 std::string outputfile = "RESLT/restriction_self_test";
2312
2313 // Loop over the levels of the hierarchy
2314 for (unsigned level = 0; level < Nlevel - 1; level++)
2315 {
2316 // Restrict the vector down to the next level
2317 Restriction_matrices_storage_pt[level]->multiply(
2318 Restriction_self_test_vectors_storage[level],
2319 Restriction_self_test_vectors_storage[level + 1]);
2320 } // End of the for loop over the hierarchy levels
2321
2322 // Loop over the levels of hierarchy to plot the restricted vectors
2323 for (unsigned level = 0; level < Nlevel; level++)
2324 {
2325 // Set output file name
2326 std::string filename = outputfile;
2327 std::stringstream string;
2328 string << level;
2329 filename += string.str();
2330 filename += ".dat";
2331
2332 // Plot the RHS vector on the current level
2333 plot(level, Restriction_self_test_vectors_storage[level], filename);
2334
2335 } // End of for-loop to output the RHS vector on each level
2336 } // End of restriction_self_test
2337
2338 //=======================================================================
2339 /// Function which implements a self-test to interpolate a
2340 /// vector up all of levels of the MG hierarchy and outputs the
2341 /// restricted vectors to file
2342 //=======================================================================
2343 template<unsigned DIM>
2345 {
2346 // Set the start of the output file name. The full names will
2347 // set after we calculate the interpolated vectors
2348 std::string outputfile = "RESLT/interpolation_self_test";
2349
2350 // Loop over the levels of the hierarchy in reverse order
2351 for (unsigned level = Nlevel - 1; level > 0; level--)
2352 {
2353 // Interpolate the vector up a level
2354 Interpolation_matrices_storage_pt[level - 1]->multiply(
2355 Interpolation_self_test_vectors_storage[level],
2356 Interpolation_self_test_vectors_storage[level - 1]);
2357 } // End of the for loop over the hierarchy levels
2358
2359 for (unsigned level = 0; level < Nlevel; level++)
2360 {
2361 // Set output file name
2362 std::string filename = outputfile;
2363 std::stringstream string;
2364 string << level;
2365 filename += string.str();
2366 filename += ".dat";
2367
2368 // Plot the RHS vector on the current level
2369 plot(level, Interpolation_self_test_vectors_storage[level], filename);
2370
2371 } // End of for-loop to output the RHS vector on each level
2372 } // End of interpolation_self_test
2373
2374 //===================================================================
2375 /// Plots the input vector (assuming we're dealing with scalar
2376 /// nodal data, otherwise I don't know how to implement this...)
2377 //===================================================================
2378 template<unsigned DIM>
2379 void MGSolver<DIM>::plot(const unsigned& hierarchy_level,
2380 const DoubleVector& input_vector,
2381 const std::string& filename)
2382 {
2383 // Setup an output file
2384 std::ofstream some_file;
2385 some_file.open(filename.c_str());
2386
2387 // Set up a pointer to the refineable mesh
2388 TreeBasedRefineableMeshBase* bulk_mesh_pt =
2389 Mg_hierarchy[hierarchy_level]->mg_bulk_mesh_pt();
2390
2391 // Find the number of elements in the bulk mesh
2392 unsigned n_el = bulk_mesh_pt->nelement();
2393
2394 // Get the dimension of the problem (assumed to be the dimension of any
2395 // node in the mesh)
2396 unsigned n_dim = bulk_mesh_pt->node_pt(0)->ndim();
2397
2398 // Find the number of nodes in an element
2399 unsigned nnod = bulk_mesh_pt->finite_element_pt(0)->nnode();
2400
2401 // Find the number of nodes in an element in one direction
2402 unsigned nnod_1d = bulk_mesh_pt->finite_element_pt(0)->nnode_1d();
2403
2404 // Loop over all elements
2405 for (unsigned e = 0; e < n_el; e++)
2406 {
2407 // Get pointer to element
2408 FiniteElement* el_pt = bulk_mesh_pt->finite_element_pt(e);
2409
2410 // Input string for tecplot zone header (when plotting nnode_1d
2411 // points in each coordinate direction)
2412 some_file << el_pt->tecplot_zone_string(nnod_1d) << std::endl;
2413
2414 // Loop over nodes
2415 for (unsigned j = 0; j < nnod; j++)
2416 {
2417 // Pointer to node
2418 Node* nod_pt = el_pt->node_pt(j);
2419
2420 // Sanity check
2421 if (nod_pt->nvalue() != 1)
2422 {
2423 std::ostringstream error_stream;
2424
2425 error_stream << "Sorry, not sure what to do here!" << nod_pt->nvalue()
2426 << std::endl;
2427
2428 // Output the value of dimension to screen
2429 error_stream << "The dimension is: " << n_dim << std::endl;
2430
2431 // Output the values of x_i for all i
2432 for (unsigned i = 0; i < n_dim; i++)
2433 {
2434 error_stream << nod_pt->x(i) << " ";
2435 }
2436
2437 // End line
2438 error_stream << std::endl;
2439
2440 // Throw an error to indicate that it was
2441 // not possible to plot the data
2442 throw OomphLibError(error_stream.str(),
2443 OOMPH_CURRENT_FUNCTION,
2444 OOMPH_EXCEPTION_LOCATION);
2445 }
2446
2447 // Output the coordinates of the node
2448 for (unsigned i = 0; i < n_dim; i++)
2449 {
2450 some_file << nod_pt->x(i) << " ";
2451 }
2452
2453 // Free or pinned?
2454 int e = nod_pt->eqn_number(0);
2455 if (e >= 0)
2456 {
2457 // Output the e-th entry if the node is free
2458 some_file << input_vector[e] << std::endl;
2459 }
2460 else
2461 {
2462 // Boundary node
2463 if (nod_pt->is_on_boundary())
2464 {
2465 // On the finest level the boundary data is pinned so set to 0
2466 if (hierarchy_level == 0)
2467 {
2468 some_file << 0.0 << std::endl;
2469 }
2470 // On any other level we output this data since it corresponds
2471 // to the correction on the boundary nodes
2472 else
2473 {
2474 some_file << input_vector[e] << std::endl;
2475 }
2476 }
2477 // Hanging node
2478 else if (nod_pt->is_hanging())
2479 {
2480 // Allocate space for a double to store the weighted sum
2481 // of the surrounding master nodes of the hanging node
2482 double hang_sum = 0.0;
2483
2484 // Find the number of master nodes of the hanging
2485 // the node in the reference element
2486 HangInfo* hang_info_pt = nod_pt->hanging_pt();
2487 unsigned nmaster = hang_info_pt->nmaster();
2488
2489 // Loop over the master nodes
2490 for (unsigned i_master = 0; i_master < nmaster; i_master++)
2491 {
2492 // Set up a pointer to the master node
2493 Node* master_node_pt = hang_info_pt->master_node_pt(i_master);
2494
2495 // Get the global equation number of this node
2496 int master_jj = master_node_pt->eqn_number(0);
2497
2498 // Is the master node a proper d.o.f.?
2499 if (master_jj >= 0)
2500 {
2501 // If the weight of the master node is non-zero
2502 if (hang_info_pt->master_weight(i_master) != 0.0)
2503 {
2504 hang_sum += hang_info_pt->master_weight(i_master) *
2505 input_vector[master_jj];
2506 }
2507 } // if (master_jj>=0)
2508 } // for (unsigned i_master=0;i_master<nmaster;i_master++)
2509
2510 // Output the weighted sum of the surrounding master nodes
2511 some_file << hang_sum << std::endl;
2512 }
2513 } // if (e>=0)
2514 } // for (unsigned j=0;j<nnod;j++)
2515 } // for (unsigned e=0;e<n_el;e++)
2516
2517 // Close output file
2518 some_file.close();
2519 } // End of plot
2520
2521 //===================================================================
2522 /// Linear solver
2523 //===================================================================
2524 template<unsigned DIM>
2526 {
2527 // If we're allowed to time things
2528 double t_mg_start = 0.0;
2529 if (!Suppress_v_cycle_output)
2530 {
2531 // Start the clock!
2532 t_mg_start = TimingHelpers::timer();
2533 }
2534
2535 // Current level
2536 unsigned finest_level = 0;
2537
2538 // Initialise the V-cycle counter
2539 V_cycle_counter = 0;
2540
2541 // Calculate the norm of the residual then output
2542 double normalised_residual_norm = residual_norm(finest_level);
2543 if (!Suppress_v_cycle_output)
2544 {
2545 oomph_info << "\nResidual on finest level for V-cycle: "
2546 << normalised_residual_norm << std::endl;
2547 }
2548
2549 // Outer loop over V-cycles
2550 //-------------------------
2551 // While the tolerance is not met and the maximum number of
2552 // V-cycles has not been completed
2553 while ((normalised_residual_norm > (this->Tolerance)) &&
2554 (V_cycle_counter != Nvcycle))
2555 {
2556 if (!Suppress_v_cycle_output)
2557 {
2558 // Output the V-cycle counter
2559 oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
2560 }
2561
2562 // Loop downwards over all levels that have coarser levels beneath them
2563 //---------------------------------------------------------------------
2564 for (unsigned i = 0; i < Nlevel - 1; i++)
2565 {
2566 // Initialise x_mg and residual_mg to 0.0 except for the finest level
2567 // since x_mg contains the current approximation to the solution and
2568 // residual_mg contains the RHS vector on the finest level
2569 if (i != 0)
2570 {
2571 // Loop over the entries in x_mg and residual_mg
2572 X_mg_vectors_storage[i].initialise(0.0);
2573 Residual_mg_vectors_storage[i].initialise(0.0);
2574 }
2575
2576 // Perform a few pre-smoothing steps and return vector that contains
2577 // the residuals of the linear system at this level.
2578 pre_smooth(i);
2579
2580 // Restrict the residual to the next coarser mesh and
2581 // assign it to the RHS vector at that level
2582 restrict_residual(i);
2583 } // Moving down the V-cycle
2584
2585 // Reached the lowest level: Do a direct solve, using the RHS vector
2586 //------------------------------------------------------------------
2587 // obtained by restriction from above.
2588 //------------------------------------
2589 direct_solve();
2590
2591 // Loop upwards over all levels that have finer levels above them
2592 //---------------------------------------------------------------
2593 for (unsigned i = Nlevel - 1; i > 0; i--)
2594 {
2595 // Interpolate solution at current level onto
2596 // next finer mesh and correct the solution x at that level
2597 interpolate_and_correct(i);
2598
2599 // Perform a few post-smoothing steps (ignore
2600 // vector that contains the residuals of the linear system
2601 // at this level)
2602 post_smooth(i - 1);
2603 }
2604
2605 // Calculate the new residual norm then output (if allowed)
2606 normalised_residual_norm = residual_norm(finest_level);
2607
2608 // If required, we will document the convergence history to screen or file
2609 // (if stream open)
2610 if (Doc_convergence_history)
2611 {
2612 if (!Output_file_stream.is_open())
2613 {
2614 oomph_info << V_cycle_counter << " " << normalised_residual_norm
2615 << std::endl;
2616 }
2617 else
2618 {
2619 Output_file_stream << V_cycle_counter << " "
2620 << normalised_residual_norm << std::endl;
2621 }
2622 } // if (Doc_convergence_history)
2623
2624 // Set counter for number of cycles (for doc)
2625 V_cycle_counter++;
2626
2627 // Print the residual on the finest level
2628 if (!Suppress_v_cycle_output)
2629 {
2630 oomph_info << "Residual on finest level of V-cycle: "
2631 << normalised_residual_norm << std::endl;
2632 }
2633 } // End of the V-cycles
2634
2635 // Copy the solution into the result vector
2636 result = X_mg_vectors_storage[finest_level];
2637
2638 // Need an extra line space if V-cycle output is suppressed
2639 if (!Suppress_v_cycle_output)
2640 {
2641 oomph_info << std::endl;
2642 }
2643
2644 // If all output is to be suppressed
2645 if (!Suppress_all_output)
2646 {
2647 // Output number of V-cycles taken to solve
2648 if (normalised_residual_norm < (this->Tolerance))
2649 {
2650 // Indicate that the problem has been solved
2651 Has_been_solved = true;
2652 }
2653 } // if (!Suppress_all_output)
2654
2655 // If the V-cycle output isn't suppressed
2656 if (!Suppress_v_cycle_output)
2657 {
2658 // Stop the clock
2659 double t_mg_end = TimingHelpers::timer();
2660 double total_G_setup_time = double(t_mg_end - t_mg_start);
2661 oomph_info << "CPU time for MG solver [sec]: " << total_G_setup_time
2662 << std::endl;
2663 }
2664 } // End of mg_solve
2665} // End of namespace oomph
2666
2667#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
unsigned nrow() const
access function to the number of global rows.
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 std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
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
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
/////////////////////////////////////////////////////////////////////
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
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
double Tolerance
Convergence tolerance.
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
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:77
An interface to allow scalar MG to be used as a Preconditioner.
MGPreconditioner(const MGPreconditioner &)=delete
Broken copy constructor.
void setup()
Function to set up a preconditioner for the linear system.
MGPreconditioner(MGProblem *mg_problem_pt)
Constructor.
void operator=(const MGPreconditioner &)=delete
Broken assignment operator.
void clean_up_memory()
Clean up memory.
~MGPreconditioner()
Destructor (empty)
virtual void preconditioner_solve(const DoubleVector &rhs, DoubleVector &z)
Function applies MG to the vector r for a full solve.
MGProblem class; subclass of Problem.
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 ...
MGProblem()
Constructor. Initialise pointers to coarser and finer levels.
virtual MGProblem * make_new_problem()=0
This function needs to be implemented in the derived problem: Returns a pointer to a new object of th...
virtual ~MGProblem()
Destructor (empty)
/////////////////////////////////////////////////////// /////////////////////////////////////////////...
void set_self_test_vector()
Makes a vector which will be used in the self-test. Is currently set to make the entries of the vecto...
Smoother *(* PostSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class Smoother to be used as the po...
void set_pre_smoother_factory_function(PreSmootherFactoryFctPt pre_smoother_fn)
Access function to set the pre-smoother creation function.
bool Suppress_all_output
If this is set to true then all output from the solver is suppressed. This is protected member data s...
void full_setup()
Do a full setup (assumes everything will be setup around the MGProblem pointer given in the construct...
Vector< DoubleVector > Rhs_mg_vectors_storage
Vector to store the RHS vectors (Rhs_mg). This is protected to allow the multigrid preconditioner to ...
void setup_mg_structures()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level.
void interpolation_self_test()
Make a self-test to make sure that the interpolation matrices are doing the same thing to interpolate...
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...
unsigned V_cycle_counter
Pointer to counter for V-cycles.
Vector< Smoother * > Post_smoothers_storage_pt
Vector to store the post-smoothers.
PostSmootherFactoryFctPt Post_smoother_factory_function_pt
Function to create post-smoothers.
Vector< DoubleVector > Residual_mg_vectors_storage
Vector to store the residual vectors.
bool Doc_everything
If this is set to true we document everything. In addition to outputting the information of the setup...
void setup_transfer_matrices()
Setup the transfer matrices on each level.
void enable_v_cycle_output()
Enable the output of the V-cycle timings and other output.
void setup_interpolation_matrices()
Setup the interpolation matrix on each level.
void post_smooth(const unsigned &level)
Post-smoother: Perform max_iter smoothing steps on the linear system Ax=b with current RHS vector,...
void set_post_smoother_factory_function(PostSmootherFactoryFctPt post_smoother_fn)
Access function to set the post-smoother creation function.
unsigned Nlevel
The number of levels in the multigrid heirachy.
void plot(const unsigned &hierarchy_level, const DoubleVector &input_vector, const std::string &filename)
Given a level in the hierarchy, an input vector and a filename this function will document the given ...
unsigned & npost_smooth()
Return the number of post-smoothing iterations (lvalue)
Vector< Smoother * > Pre_smoothers_storage_pt
Vector to store the pre-smoothers.
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout. This is protected member data to allow the prec...
void set_restriction_matrices_as_interpolation_transposes()
Builds a CRDoubleMatrix on each level that is used to restrict the residual between levels....
Smoother *(* PreSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class Smoother to be used as the pr...
MGSolver(MGProblem *mg_problem_pt)
Constructor: Set up default values for number of V-cycles and pre- and post-smoothing steps.
void solve(Problem *const &problem_pt, DoubleVector &result)
Virtual function in the base class that needs to be implemented later but for now just leave it empty...
void mg_solve(DoubleVector &result)
Do the actual solve – this is called through the pure virtual solve function in the LinearSolver base...
void disable_output()
Suppress anything that can be suppressed, i.e. any timings. Things like mesh adaptation can not howev...
void restrict_residual(const unsigned &level)
Restrict residual (computed on level-th MG level) to the next coarser mesh and stick it into the coar...
Vector< DoubleVector > Restriction_self_test_vectors_storage
Vector to store the result of restriction on each level (only required if the user wishes to document...
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...
void disable_smoother_and_superlu_doc_time()
Suppress the output of both smoothers and SuperLU.
void clean_up_memory()
Clean up anything that needs to be cleaned up.
void self_test()
Makes a vector, restricts it down the levels of the hierarchy and documents it at each level....
Vector< CRDoubleMatrix * > Restriction_matrices_storage_pt
Vector to store the restriction matrices.
void setup_mg_hierarchy()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level.
Vector< CRDoubleMatrix * > Mg_matrices_storage_pt
Vector to store the system matrices.
void setup_interpolation_matrices_unstructured()
Setup the interpolation matrix on each level (used for unstructured meshes)
void direct_solve()
Call the direct solver (SuperLU) to solve the problem exactly.
void setup_smoothers()
Function to set up all of the smoothers once the system matrices have been set up.
unsigned Npost_smooth
Number of post-smoothing steps.
void disable_v_cycle_output()
Disable all output from mg_solve apart from the number of V-cycles used to solve the problem.
~MGSolver()
Delete any dynamically allocated data.
unsigned Nvcycle
Maximum number of V-cycles (this is set as a protected variable so.
Vector< MGProblem * > Mg_hierarchy
Vector containing pointers to problems in hierarchy.
Vector< DoubleVector > Interpolation_self_test_vectors_storage
Vector to store the result of interpolation on each level (only required if the user wishes to docume...
Vector< CRDoubleMatrix * > Interpolation_matrices_storage_pt
Vector to store the interpolation matrices.
void enable_output()
Enable the output from anything that could have been suppressed.
void enable_doc_everything()
Enable the output from anything that could have been suppressed.
void pre_smooth(const unsigned &level)
Pre-smoother: Perform 'max_iter' smoothing steps on the linear system Ax=b with current RHS vector,...
bool Suppress_v_cycle_output
Indicates whether or not the V-cycle output should be suppressed. Needs to be protected member data f...
unsigned Npre_smooth
Number of pre-smoothing steps.
unsigned & npre_smooth()
Return the number of pre-smoothing iterations (lvalue)
unsigned iterations() const
Number of iterations.
void restriction_self_test()
Make a self-test to make sure that the interpolation matrices are doing the same thing to restrict th...
Vector< DoubleVector > X_mg_vectors_storage
Vector to store the solution vectors (X_mg)
void modify_restriction_matrices()
Normalise the rows of the restriction matrices to avoid amplifications when projecting to the coarser...
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...
void interpolate_and_correct(const unsigned &level)
Interpolate solution at current level onto next finer mesh and correct the solution x at that level.
double residual_norm(const unsigned &level)
Return norm of residual r=b-Ax and the residual vector itself on the level-th level.
bool Has_been_setup
Boolean variable to indicate whether or not the solver has been setup.
MGProblem * Mg_problem_pt
Pointer to the MG problem (deep copy). This is protected to provide access to the MG preconditioner.
bool Has_been_solved
Boolean variable to indicate whether or not the problem was successfully solved.
PreSmootherFactoryFctPt Pre_smoother_factory_function_pt
Function to create pre-smoothers.
unsigned & max_iter()
Number of iterations.
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
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
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2499
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
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
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....
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1674
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
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1280
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Base class for refineable meshes. Provides standardised interfaces for the following standard mesh ad...
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
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const double Pi
50 digits from maple
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...