Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-63.h
Go to the documentation of this file.
1
520 *   ParameterHandler prm;
521 *  
522 *   prm.declare_entry("Epsilon",
523 *   "0.005",
524 *   Patterns::Double(0),
525 *   "Diffusion parameter");
526 *  
527 *   prm.declare_entry("Fe degree",
528 *   "1",
529 *   Patterns::Integer(1),
530 *   "Finite Element degree");
531 *   prm.declare_entry("Smoother type",
532 *   "block SOR",
533 *   Patterns::Selection("SOR|Jacobi|block SOR|block Jacobi"),
534 *   "Select smoother: SOR|Jacobi|block SOR|block Jacobi");
535 *   prm.declare_entry("Smoothing steps",
536 *   "2",
537 *   Patterns::Integer(1),
538 *   "Number of smoothing steps");
539 *   prm.declare_entry(
540 *   "DoF renumbering",
541 *   "downstream",
542 *   Patterns::Selection("none|downstream|upstream|random"),
543 *   "Select DoF renumbering: none|downstream|upstream|random");
544 *   prm.declare_entry("With streamline diffusion",
545 *   "true",
546 *   Patterns::Bool(),
547 *   "Enable streamline diffusion stabilization: true|false");
548 *   prm.declare_entry("Output",
549 *   "true",
550 *   Patterns::Bool(),
551 *   "Generate graphical output: true|false");
552 *  
553 *   /* ...and then try to read their values from the input file: */
554 *   if (prm_filename.empty())
555 *   {
556 *   prm.print_parameters(std::cout, ParameterHandler::Text);
557 *   AssertThrow(
558 *   false, ExcMessage("Please pass a .prm file as the first argument!"));
559 *   }
560 *  
561 *   prm.parse_input(prm_filename);
562 *  
563 *   epsilon = prm.get_double("Epsilon");
564 *   fe_degree = prm.get_integer("Fe degree");
565 *   smoother_type = prm.get("Smoother type");
566 *   smoothing_steps = prm.get_integer("Smoothing steps");
567 *  
568 *   const std::string renumbering = prm.get("DoF renumbering");
569 *   if (renumbering == "none")
570 *   dof_renumbering = DoFRenumberingStrategy::none;
571 *   else if (renumbering == "downstream")
572 *   dof_renumbering = DoFRenumberingStrategy::downstream;
573 *   else if (renumbering == "upstream")
574 *   dof_renumbering = DoFRenumberingStrategy::upstream;
575 *   else if (renumbering == "random")
576 *   dof_renumbering = DoFRenumberingStrategy::random;
577 *   else
578 *   AssertThrow(false,
579 *   ExcMessage("The <DoF renumbering> parameter has "
580 *   "an invalid value."));
581 *  
582 *   with_streamline_diffusion = prm.get_bool("With streamline diffusion");
583 *   output = prm.get_bool("Output");
584 *   }
585 *  
586 *  
587 * @endcode
588 *
589 *
590 * <a name="Cellpermutations"></a>
591 * <h3>Cell permutations</h3>
592 *
593
594 *
595 * The ordering in which cells and degrees of freedom are traversed
596 * will play a role in the speed of convergence for multiplicative
597 * methods. Here we define functions which return a specific ordering
598 * of cells to be used by the block smoothers.
599 *
600
601 *
602 * For each type of cell ordering, we define a function for the
603 * active mesh and one for a level mesh (i.e., for the cells at one
604 * level of a multigrid hierarchy). While the only reordering
605 * necessary for solving the system will be on the level meshes, we
606 * include the active reordering for visualization purposes in
607 * output_results().
608 *
609
610 *
611 * For the two downstream ordering functions, we first create an
612 * array with all of the relevant cells that we then sort in
613 * downstream direction using a "comparator" object. The output of
614 * the functions is then simply an array of the indices of the cells
615 * in the just computed order.
616 *
617 * @code
618 *   template <int dim>
619 *   std::vector<unsigned int>
620 *   create_downstream_cell_ordering(const DoFHandler<dim> &dof_handler,
621 *   const Tensor<1, dim> direction,
622 *   const unsigned int level)
623 *   {
624 *   std::vector<typename DoFHandler<dim>::level_cell_iterator> ordered_cells;
625 *   ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level));
626 *   for (const auto &cell : dof_handler.cell_iterators_on_level(level))
627 *   ordered_cells.push_back(cell);
628 *  
629 *   const DoFRenumbering::
630 *   CompareDownstream<typename DoFHandler<dim>::level_cell_iterator, dim>
631 *   comparator(direction);
632 *   std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
633 *  
634 *   std::vector<unsigned> ordered_indices;
635 *   ordered_indices.reserve(dof_handler.get_triangulation().n_cells(level));
636 *  
637 *   for (const auto &cell : ordered_cells)
638 *   ordered_indices.push_back(cell->index());
639 *  
640 *   return ordered_indices;
641 *   }
642 *  
643 *  
644 *  
645 *   template <int dim>
646 *   std::vector<unsigned int>
647 *   create_downstream_cell_ordering(const DoFHandler<dim> &dof_handler,
648 *   const Tensor<1, dim> direction)
649 *   {
650 *   std::vector<typename DoFHandler<dim>::active_cell_iterator> ordered_cells;
651 *   ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
652 *   for (const auto &cell : dof_handler.active_cell_iterators())
653 *   ordered_cells.push_back(cell);
654 *  
655 *   const DoFRenumbering::
656 *   CompareDownstream<typename DoFHandler<dim>::active_cell_iterator, dim>
657 *   comparator(direction);
658 *   std::sort(ordered_cells.begin(), ordered_cells.end(), comparator);
659 *  
660 *   std::vector<unsigned int> ordered_indices;
661 *   ordered_indices.reserve(dof_handler.get_triangulation().n_active_cells());
662 *  
663 *   for (const auto &cell : ordered_cells)
664 *   ordered_indices.push_back(cell->index());
665 *  
666 *   return ordered_indices;
667 *   }
668 *  
669 *  
670 * @endcode
671 *
672 * The functions that produce a random ordering are similar in
673 * spirit in that they first put information about all cells into an
674 * array. But then, instead of sorting them, they shuffle the
675 * elements randomly using the facilities C++ offers to generate
676 * random numbers. The way this is done is by iterating over all
677 * elements of the array, drawing a random number for another
678 * element before that, and then exchanging these elements. The
679 * result is a random shuffle of the elements of the array.
680 *
681 * @code
682 *   template <int dim>
683 *   std::vector<unsigned int>
684 *   create_random_cell_ordering(const DoFHandler<dim> &dof_handler,
685 *   const unsigned int level)
686 *   {
687 *   std::vector<unsigned int> ordered_cells;
688 *   ordered_cells.reserve(dof_handler.get_triangulation().n_cells(level));
689 *   for (const auto &cell : dof_handler.cell_iterators_on_level(level))
690 *   ordered_cells.push_back(cell->index());
691 *  
692 *   std::mt19937 random_number_generator;
693 *   std::shuffle(ordered_cells.begin(),
694 *   ordered_cells.end(),
695 *   random_number_generator);
696 *  
697 *   return ordered_cells;
698 *   }
699 *  
700 *  
701 *  
702 *   template <int dim>
703 *   std::vector<unsigned int>
704 *   create_random_cell_ordering(const DoFHandler<dim> &dof_handler)
705 *   {
706 *   std::vector<unsigned int> ordered_cells;
707 *   ordered_cells.reserve(dof_handler.get_triangulation().n_active_cells());
708 *   for (const auto &cell : dof_handler.active_cell_iterators())
709 *   ordered_cells.push_back(cell->index());
710 *  
711 *   std::mt19937 random_number_generator;
712 *   std::shuffle(ordered_cells.begin(),
713 *   ordered_cells.end(),
714 *   random_number_generator);
715 *  
716 *   return ordered_cells;
717 *   }
718 *  
719 *  
720 * @endcode
721 *
722 *
723 * <a name="Righthandsideandboundaryvalues"></a>
724 * <h3>Right-hand side and boundary values</h3>
725 *
726
727 *
728 * The problem solved in this tutorial is an adaptation of Ex. 3.1.3 found
729 * on pg. 118 of <a
730 * href="https://global.oup.com/academic/product/finite-elements-and-fast-iterative-solvers-9780199678808">
731 * Finite Elements and Fast Iterative Solvers: with Applications in
732 * Incompressible Fluid Dynamics by Elman, Silvester, and Wathen</a>. The
733 * main difference being that we add a hole in the center of our domain with
734 * zero Dirichlet boundary conditions.
735 *
736
737 *
738 * For a complete description, we need classes that implement the
739 * zero right-hand side first (we could of course have just used
741 *
742 * @code
743 *   template <int dim>
744 *   class RightHandSide : public Function<dim>
745 *   {
746 *   public:
747 *   virtual double value(const Point<dim> & p,
748 *   const unsigned int component = 0) const override;
749 *  
750 *   virtual void value_list(const std::vector<Point<dim>> &points,
751 *   std::vector<double> & values,
752 *   const unsigned int component = 0) const override;
753 *   };
754 *  
755 *  
756 *  
757 *   template <int dim>
758 *   double RightHandSide<dim>::value(const Point<dim> &,
759 *   const unsigned int component) const
760 *   {
761 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
762 *   (void)component;
763 *  
764 *   return 0.0;
765 *   }
766 *  
767 *  
768 *  
769 *   template <int dim>
770 *   void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
771 *   std::vector<double> & values,
772 *   const unsigned int component) const
773 *   {
774 *   AssertDimension(values.size(), points.size());
775 *  
776 *   for (unsigned int i = 0; i < points.size(); ++i)
777 *   values[i] = RightHandSide<dim>::value(points[i], component);
778 *   }
779 *  
780 *  
781 * @endcode
782 *
783 * We also have Dirichlet boundary conditions. On a connected portion of the
784 * outer, square boundary we set the value to 1, and we set the value to 0
785 * everywhere else (including the inner, circular boundary):
786 *
787 * @code
788 *   template <int dim>
789 *   class BoundaryValues : public Function<dim>
790 *   {
791 *   public:
792 *   virtual double value(const Point<dim> & p,
793 *   const unsigned int component = 0) const override;
794 *  
795 *   virtual void value_list(const std::vector<Point<dim>> &points,
796 *   std::vector<double> & values,
797 *   const unsigned int component = 0) const override;
798 *   };
799 *  
800 *  
801 *  
802 *   template <int dim>
803 *   double BoundaryValues<dim>::value(const Point<dim> & p,
804 *   const unsigned int component) const
805 *   {
806 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
807 *   (void)component;
808 *  
809 * @endcode
810 *
811 * Set boundary to 1 if @f$x=1@f$, or if @f$x>0.5@f$ and @f$y=-1@f$.
812 *
813 * @code
814 *   if (std::fabs(p[0] - 1) < 1e-8 ||
815 *   (std::fabs(p[1] + 1) < 1e-8 && p[0] >= 0.5))
816 *   {
817 *   return 1.0;
818 *   }
819 *   else
820 *   {
821 *   return 0.0;
822 *   }
823 *   }
824 *  
825 *  
826 *  
827 *   template <int dim>
828 *   void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
829 *   std::vector<double> & values,
830 *   const unsigned int component) const
831 *   {
832 *   AssertDimension(values.size(), points.size());
833 *  
834 *   for (unsigned int i = 0; i < points.size(); ++i)
835 *   values[i] = BoundaryValues<dim>::value(points[i], component);
836 *   }
837 *  
838 *  
839 *  
840 * @endcode
841 *
842 *
843 * <a name="Streamlinediffusionimplementation"></a>
844 * <h3>Streamline diffusion implementation</h3>
845 *
846
847 *
848 * The streamline diffusion method has a stabilization constant that
849 * we need to be able to compute. The choice of how this parameter
850 * is computed is taken from <a
851 * href="https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">On
852 * Discontinuity-Capturing Methods for Convection-Diffusion
853 * Equations by Volker John and Petr Knobloch</a>.
854 *
855 * @code
856 *   template <int dim>
857 *   double compute_stabilization_delta(const double hk,
858 *   const double eps,
859 *   const Tensor<1, dim> dir,
860 *   const double pk)
861 *   {
862 *   const double Peclet = dir.norm() * hk / (2.0 * eps * pk);
863 *   const double coth =
864 *   (1.0 + std::exp(-2.0 * Peclet)) / (1.0 - std::exp(-2.0 * Peclet));
865 *  
866 *   return hk / (2.0 * dir.norm() * pk) * (coth - 1.0 / Peclet);
867 *   }
868 *  
869 *  
870 * @endcode
871 *
872 *
873 * <a name="codeAdvectionProlemcodeclass"></a>
874 * <h3><code>AdvectionProlem</code> class</h3>
875 *
876
877 *
878 * This is the main class of the program, and should look very similar to
879 * @ref step_16 "step-16". The major difference is that, since we are defining our multigrid
880 * smoother at runtime, we choose to define a function `create_smoother()` and
881 * a class object `mg_smoother` which is a `std::unique_ptr` to a smoother
882 * that is derived from MGSmoother. Note that for smoothers derived from
883 * RelaxationBlock, we must include a `smoother_data` object for each level.
884 * This will contain information about the cell ordering and the method of
885 * inverting cell matrices.
886 *
887
888 *
889 *
890 * @code
891 *   template <int dim>
892 *   class AdvectionProblem
893 *   {
894 *   public:
895 *   AdvectionProblem(const Settings &settings);
896 *   void run();
897 *  
898 *   private:
899 *   void setup_system();
900 *  
901 *   template <class IteratorType>
902 *   void assemble_cell(const IteratorType &cell,
903 *   ScratchData<dim> & scratch_data,
904 *   CopyData & copy_data);
905 *   void assemble_system_and_multigrid();
906 *  
907 *   void setup_smoother();
908 *  
909 *   void solve();
910 *   void refine_grid();
911 *   void output_results(const unsigned int cycle) const;
912 *  
914 *   DoFHandler<dim> dof_handler;
915 *  
916 *   const FE_Q<dim> fe;
917 *   const MappingQ<dim> mapping;
918 *  
919 *   AffineConstraints<double> constraints;
920 *  
921 *   SparsityPattern sparsity_pattern;
922 *   SparseMatrix<double> system_matrix;
923 *  
924 *   Vector<double> solution;
925 *   Vector<double> system_rhs;
926 *  
927 *   MGLevelObject<SparsityPattern> mg_sparsity_patterns;
928 *   MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
929 *  
930 *   MGLevelObject<SparseMatrix<double>> mg_matrices;
931 *   MGLevelObject<SparseMatrix<double>> mg_interface_in;
932 *   MGLevelObject<SparseMatrix<double>> mg_interface_out;
933 *  
934 *   mg::Matrix<Vector<double>> mg_matrix;
935 *   mg::Matrix<Vector<double>> mg_interface_matrix_in;
936 *   mg::Matrix<Vector<double>> mg_interface_matrix_out;
937 *  
938 *   std::unique_ptr<MGSmoother<Vector<double>>> mg_smoother;
939 *  
940 *   using SmootherType =
942 *   using SmootherAdditionalDataType = SmootherType::AdditionalData;
944 *  
945 *   MGConstrainedDoFs mg_constrained_dofs;
946 *  
947 *   Tensor<1, dim> advection_direction;
948 *  
949 *   const Settings settings;
950 *   };
951 *  
952 *  
953 *  
954 *   template <int dim>
955 *   AdvectionProblem<dim>::AdvectionProblem(const Settings &settings)
957 *   , dof_handler(triangulation)
958 *   , fe(settings.fe_degree)
959 *   , mapping(settings.fe_degree)
960 *   , settings(settings)
961 *   {
962 *   advection_direction[0] = -std::sin(numbers::PI / 6.0);
963 *   if (dim >= 2)
964 *   advection_direction[1] = std::cos(numbers::PI / 6.0);
965 *   if (dim >= 3)
966 *   AssertThrow(false, ExcNotImplemented());
967 *   }
968 *  
969 *  
970 * @endcode
971 *
972 *
973 * <a name="codeAdvectionProblemsetup_systemcode"></a>
974 * <h4><code>AdvectionProblem::setup_system()</code></h4>
975 *
976
977 *
978 * Here we first set up the DoFHandler, AffineConstraints, and
979 * SparsityPattern objects for both active and multigrid level meshes.
980 *
981
982 *
983 * We could renumber the active DoFs with the DoFRenumbering class,
984 * but the smoothers only act on multigrid levels and as such, this
985 * would not matter for the computations. Instead, we will renumber the
986 * DoFs on each multigrid level below.
987 *
988 * @code
989 *   template <int dim>
990 *   void AdvectionProblem<dim>::setup_system()
991 *   {
992 *   const unsigned int n_levels = triangulation.n_levels();
993 *  
994 *   dof_handler.distribute_dofs(fe);
995 *  
996 *   solution.reinit(dof_handler.n_dofs());
997 *   system_rhs.reinit(dof_handler.n_dofs());
998 *  
999 *   constraints.clear();
1000 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1001 *  
1003 *   mapping, dof_handler, 0, BoundaryValues<dim>(), constraints);
1005 *   mapping, dof_handler, 1, BoundaryValues<dim>(), constraints);
1006 *   constraints.close();
1007 *  
1008 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1009 *   DoFTools::make_sparsity_pattern(dof_handler,
1010 *   dsp,
1011 *   constraints,
1012 *   /*keep_constrained_dofs = */ false);
1013 *  
1014 *   sparsity_pattern.copy_from(dsp);
1015 *   system_matrix.reinit(sparsity_pattern);
1016 *  
1017 *   dof_handler.distribute_mg_dofs();
1018 *  
1019 * @endcode
1020 *
1021 * Having enumerated the global degrees of freedom as well as (in
1022 * the last line above) the level degrees of freedom, let us
1023 * renumber the level degrees of freedom to get a better smoother
1024 * as explained in the introduction. The first block below
1025 * renumbers DoFs on each level in downstream or upstream
1026 * direction if needed. This is only necessary for point smoothers
1027 * (SOR and Jacobi) as the block smoothers operate on cells (see
1028 * `create_smoother()`). The blocks below then also implement
1029 * random numbering.
1030 *
1031 * @code
1032 *   if (settings.smoother_type == "SOR" || settings.smoother_type == "Jacobi")
1033 *   {
1034 *   if (settings.dof_renumbering ==
1035 *   Settings::DoFRenumberingStrategy::downstream ||
1036 *   settings.dof_renumbering ==
1037 *   Settings::DoFRenumberingStrategy::upstream)
1038 *   {
1039 *   const Tensor<1, dim> direction =
1040 *   (settings.dof_renumbering ==
1041 *   Settings::DoFRenumberingStrategy::upstream ?
1042 *   -1.0 :
1043 *   1.0) *
1044 *   advection_direction;
1045 *  
1046 *   for (unsigned int level = 0; level < n_levels; ++level)
1047 *   DoFRenumbering::downstream(dof_handler,
1048 *   level,
1049 *   direction,
1050 *   /*dof_wise_renumbering = */ true);
1051 *   }
1052 *   else if (settings.dof_renumbering ==
1053 *   Settings::DoFRenumberingStrategy::random)
1054 *   {
1055 *   for (unsigned int level = 0; level < n_levels; ++level)
1056 *   DoFRenumbering::random(dof_handler, level);
1057 *   }
1058 *   else
1059 *   Assert(false, ExcNotImplemented());
1060 *   }
1061 *  
1062 * @endcode
1063 *
1064 * The rest of the function just sets up data structures. The last
1065 * lines of the code below is unlike the other GMG tutorials, as
1066 * it sets up both the interface in and out matrices. We need this
1067 * since our problem is non-symmetric.
1068 *
1069 * @code
1070 *   mg_constrained_dofs.clear();
1071 *   mg_constrained_dofs.initialize(dof_handler);
1072 *  
1073 *   mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, {0, 1});
1074 *  
1075 *   mg_matrices.resize(0, n_levels - 1);
1076 *   mg_matrices.clear_elements();
1077 *   mg_interface_in.resize(0, n_levels - 1);
1078 *   mg_interface_in.clear_elements();
1079 *   mg_interface_out.resize(0, n_levels - 1);
1080 *   mg_interface_out.clear_elements();
1081 *   mg_sparsity_patterns.resize(0, n_levels - 1);
1082 *   mg_interface_sparsity_patterns.resize(0, n_levels - 1);
1083 *  
1084 *   for (unsigned int level = 0; level < n_levels; ++level)
1085 *   {
1086 *   {
1087 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
1088 *   dof_handler.n_dofs(level));
1089 *   MGTools::make_sparsity_pattern(dof_handler, dsp, level);
1090 *   mg_sparsity_patterns[level].copy_from(dsp);
1091 *   mg_matrices[level].reinit(mg_sparsity_patterns[level]);
1092 *   }
1093 *   {
1094 *   DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
1095 *   dof_handler.n_dofs(level));
1097 *   mg_constrained_dofs,
1098 *   dsp,
1099 *   level);
1100 *   mg_interface_sparsity_patterns[level].copy_from(dsp);
1101 *  
1102 *   mg_interface_in[level].reinit(mg_interface_sparsity_patterns[level]);
1103 *   mg_interface_out[level].reinit(mg_interface_sparsity_patterns[level]);
1104 *   }
1105 *   }
1106 *   }
1107 *  
1108 *  
1109 * @endcode
1110 *
1111 *
1112 * <a name="codeAdvectionProblemassemble_cellcode"></a>
1113 * <h4><code>AdvectionProblem::assemble_cell()</code></h4>
1114 *
1115
1116 *
1117 * Here we define the assembly of the linear system on each cell to
1118 * be used by the mesh_loop() function below. This one function
1119 * assembles the cell matrix for either an active or a level cell
1120 * (whatever it is passed as its first argument), and only assembles
1121 * a right-hand side if called with an active cell.
1122 *
1123
1124 *
1125 *
1126 * @code
1127 *   template <int dim>
1128 *   template <class IteratorType>
1129 *   void AdvectionProblem<dim>::assemble_cell(const IteratorType &cell,
1130 *   ScratchData<dim> & scratch_data,
1131 *   CopyData & copy_data)
1132 *   {
1133 *   copy_data.level = cell->level();
1134 *  
1135 *   const unsigned int dofs_per_cell =
1136 *   scratch_data.fe_values.get_fe().n_dofs_per_cell();
1137 *   copy_data.dofs_per_cell = dofs_per_cell;
1138 *   copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
1139 *  
1140 *   const unsigned int n_q_points =
1141 *   scratch_data.fe_values.get_quadrature().size();
1142 *  
1143 *   if (cell->is_level_cell() == false)
1144 *   copy_data.cell_rhs.reinit(dofs_per_cell);
1145 *  
1146 *   copy_data.local_dof_indices.resize(dofs_per_cell);
1147 *   cell->get_active_or_mg_dof_indices(copy_data.local_dof_indices);
1148 *  
1149 *   scratch_data.fe_values.reinit(cell);
1150 *  
1151 *   RightHandSide<dim> right_hand_side;
1152 *   std::vector<double> rhs_values(n_q_points);
1153 *  
1154 *   right_hand_side.value_list(scratch_data.fe_values.get_quadrature_points(),
1155 *   rhs_values);
1156 *  
1157 * @endcode
1158 *
1159 * If we are using streamline diffusion we must add its contribution
1160 * to both the cell matrix and the cell right-hand side. If we are not
1161 * using streamline diffusion, setting @f$\delta=0@f$ negates this contribution
1162 * below and we are left with the standard, Galerkin finite element
1163 * assembly.
1164 *
1165 * @code
1166 *   const double delta = (settings.with_streamline_diffusion ?
1167 *   compute_stabilization_delta(cell->diameter(),
1168 *   settings.epsilon,
1169 *   advection_direction,
1170 *   settings.fe_degree) :
1171 *   0.0);
1172 *  
1173 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1174 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1175 *   {
1176 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1177 *   {
1178 * @endcode
1179 *
1180 * The assembly of the local matrix has two parts. First
1181 * the Galerkin contribution:
1182 *
1183 * @code
1184 *   copy_data.cell_matrix(i, j) +=
1185 *   (settings.epsilon *
1186 *   scratch_data.fe_values.shape_grad(i, q_point) *
1187 *   scratch_data.fe_values.shape_grad(j, q_point) *
1188 *   scratch_data.fe_values.JxW(q_point)) +
1189 *   (scratch_data.fe_values.shape_value(i, q_point) *
1190 *   (advection_direction *
1191 *   scratch_data.fe_values.shape_grad(j, q_point)) *
1192 *   scratch_data.fe_values.JxW(q_point))
1193 * @endcode
1194 *
1195 * and then the streamline diffusion contribution:
1196 *
1197 * @code
1198 *   + delta *
1199 *   (advection_direction *
1200 *   scratch_data.fe_values.shape_grad(j, q_point)) *
1201 *   (advection_direction *
1202 *   scratch_data.fe_values.shape_grad(i, q_point)) *
1203 *   scratch_data.fe_values.JxW(q_point) -
1204 *   delta * settings.epsilon *
1205 *   trace(scratch_data.fe_values.shape_hessian(j, q_point)) *
1206 *   (advection_direction *
1207 *   scratch_data.fe_values.shape_grad(i, q_point)) *
1208 *   scratch_data.fe_values.JxW(q_point);
1209 *   }
1210 *   if (cell->is_level_cell() == false)
1211 *   {
1212 * @endcode
1213 *
1214 * The same applies to the right hand side. First the
1215 * Galerkin contribution:
1216 *
1217 * @code
1218 *   copy_data.cell_rhs(i) +=
1219 *   scratch_data.fe_values.shape_value(i, q_point) *
1220 *   rhs_values[q_point] * scratch_data.fe_values.JxW(q_point)
1221 * @endcode
1222 *
1223 * and then the streamline diffusion contribution:
1224 *
1225 * @code
1226 *   + delta * rhs_values[q_point] * advection_direction *
1227 *   scratch_data.fe_values.shape_grad(i, q_point) *
1228 *   scratch_data.fe_values.JxW(q_point);
1229 *   }
1230 *   }
1231 *   }
1232 *  
1233 *  
1234 * @endcode
1235 *
1236 *
1237 * <a name="codeAdvectionProblemassemble_system_and_multigridcode"></a>
1238 * <h4><code>AdvectionProblem::assemble_system_and_multigrid()</code></h4>
1239 *
1240
1241 *
1242 * Here we employ MeshWorker::mesh_loop() to go over cells and assemble the
1243 * system_matrix, system_rhs, and all mg_matrices for us.
1244 *
1245
1246 *
1247 *
1248 * @code
1249 *   template <int dim>
1250 *   void AdvectionProblem<dim>::assemble_system_and_multigrid()
1251 *   {
1252 *   const auto cell_worker_active =
1253 *   [&](const decltype(dof_handler.begin_active()) &cell,
1254 *   ScratchData<dim> & scratch_data,
1255 *   CopyData & copy_data) {
1256 *   this->assemble_cell(cell, scratch_data, copy_data);
1257 *   };
1258 *  
1259 *   const auto copier_active = [&](const CopyData &copy_data) {
1260 *   constraints.distribute_local_to_global(copy_data.cell_matrix,
1261 *   copy_data.cell_rhs,
1262 *   copy_data.local_dof_indices,
1263 *   system_matrix,
1264 *   system_rhs);
1265 *   };
1266 *  
1267 *  
1268 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1269 *   dof_handler.end(),
1270 *   cell_worker_active,
1271 *   copier_active,
1272 *   ScratchData<dim>(fe, fe.degree + 1),
1273 *   CopyData(),
1275 *  
1276 * @endcode
1277 *
1278 * Unlike the constraints for the active level, we choose to create
1279 * constraint objects for each multigrid level local to this function
1280 * since they are never needed elsewhere in the program.
1281 *
1282 * @code
1283 *   std::vector<AffineConstraints<double>> boundary_constraints(
1284 *   triangulation.n_global_levels());
1285 *   for (unsigned int level = 0; level < triangulation.n_global_levels();
1286 *   ++level)
1287 *   {
1288 *   const IndexSet locally_owned_level_dof_indices =
1290 *   boundary_constraints[level].reinit(locally_owned_level_dof_indices);
1291 *   boundary_constraints[level].add_lines(
1292 *   mg_constrained_dofs.get_refinement_edge_indices(level));
1293 *   boundary_constraints[level].add_lines(
1294 *   mg_constrained_dofs.get_boundary_indices(level));
1295 *   boundary_constraints[level].close();
1296 *   }
1297 *  
1298 *   const auto cell_worker_mg =
1299 *   [&](const decltype(dof_handler.begin_mg()) &cell,
1300 *   ScratchData<dim> & scratch_data,
1301 *   CopyData & copy_data) {
1302 *   this->assemble_cell(cell, scratch_data, copy_data);
1303 *   };
1304 *  
1305 *   const auto copier_mg = [&](const CopyData &copy_data) {
1306 *   boundary_constraints[copy_data.level].distribute_local_to_global(
1307 *   copy_data.cell_matrix,
1308 *   copy_data.local_dof_indices,
1309 *   mg_matrices[copy_data.level]);
1310 *  
1311 * @endcode
1312 *
1313 * If @f$(i,j)@f$ is an `interface_out` dof pair, then @f$(j,i)@f$ is an
1314 * `interface_in` dof pair. Note: For `interface_in`, we load
1315 * the transpose of the interface entries, i.e., the entry for
1316 * dof pair @f$(j,i)@f$ is stored in `interface_in(i,j)`. This is an
1317 * optimization for the symmetric case which allows only one
1318 * matrix to be used when setting the edge_matrices in
1319 * solve(). Here, however, since our problem is non-symmetric,
1320 * we must store both `interface_in` and `interface_out`
1321 * matrices.
1322 *
1323 * @code
1324 *   for (unsigned int i = 0; i < copy_data.dofs_per_cell; ++i)
1325 *   for (unsigned int j = 0; j < copy_data.dofs_per_cell; ++j)
1326 *   if (mg_constrained_dofs.is_interface_matrix_entry(
1327 *   copy_data.level,
1328 *   copy_data.local_dof_indices[i],
1329 *   copy_data.local_dof_indices[j]))
1330 *   {
1331 *   mg_interface_out[copy_data.level].add(
1332 *   copy_data.local_dof_indices[i],
1333 *   copy_data.local_dof_indices[j],
1334 *   copy_data.cell_matrix(i, j));
1335 *   mg_interface_in[copy_data.level].add(
1336 *   copy_data.local_dof_indices[i],
1337 *   copy_data.local_dof_indices[j],
1338 *   copy_data.cell_matrix(j, i));
1339 *   }
1340 *   };
1341 *  
1342 *   MeshWorker::mesh_loop(dof_handler.begin_mg(),
1343 *   dof_handler.end_mg(),
1344 *   cell_worker_mg,
1345 *   copier_mg,
1346 *   ScratchData<dim>(fe, fe.degree + 1),
1347 *   CopyData(),
1349 *   }
1350 *  
1351 *  
1352 * @endcode
1353 *
1354 *
1355 * <a name="codeAdvectionProblemsetup_smoothercode"></a>
1356 * <h4><code>AdvectionProblem::setup_smoother()</code></h4>
1357 *
1358
1359 *
1360 * Next, we set up the smoother based on the settings in the `.prm` file. The
1361 * two options that are of significance is the number of pre- and
1362 * post-smoothing steps on each level of the multigrid v-cycle and the
1363 * relaxation parameter.
1364 *
1365
1366 *
1367 * Since multiplicative methods tend to be more powerful than additive method,
1368 * fewer smoothing steps are required to see convergence independent of mesh
1369 * size. The same holds for block smoothers over point smoothers. This is
1370 * reflected in the choice for the number of smoothing steps for each type of
1371 * smoother below.
1372 *
1373
1374 *
1375 * The relaxation parameter for point smoothers is chosen based on trial and
1376 * error, and reflects values necessary to keep the iteration counts in
1377 * the GMRES solve constant (or as close as possible) as we refine the mesh.
1378 * The two values given for both "Jacobi" and "SOR" in the `.prm` files are
1379 * for degree 1 and degree 3 finite elements. If the user wants to change to
1380 * another degree, they may need to adjust these numbers. For block smoothers,
1381 * this parameter has a more straightforward interpretation, namely that for
1382 * additive methods in 2d, a DoF can have a repeated contribution from up to 4
1383 * cells, therefore we must relax these methods by 0.25 to compensate. This is
1384 * not an issue for multiplicative methods as each cell's inverse application
1385 * carries new information to all its DoFs.
1386 *
1387
1388 *
1389 * Finally, as mentioned above, the point smoothers only operate on DoFs, and
1390 * the block smoothers on cells, so only the block smoothers need to be given
1391 * information regarding cell orderings. DoF ordering for point smoothers has
1392 * already been taken care of in `setup_system()`.
1393 *
1394
1395 *
1396 *
1397 * @code
1398 *   template <int dim>
1399 *   void AdvectionProblem<dim>::setup_smoother()
1400 *   {
1401 *   if (settings.smoother_type == "SOR")
1402 *   {
1403 *   using Smoother = PreconditionSOR<SparseMatrix<double>>;
1404 *  
1405 *   auto smoother =
1406 *   std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1407 *   Smoother,
1408 *   Vector<double>>>();
1409 *   smoother->initialize(mg_matrices,
1410 *   Smoother::AdditionalData(fe.degree == 1 ? 1.0 :
1411 *   0.62));
1412 *   smoother->set_steps(settings.smoothing_steps);
1413 *   mg_smoother = std::move(smoother);
1414 *   }
1415 *   else if (settings.smoother_type == "Jacobi")
1416 *   {
1417 *   using Smoother = PreconditionJacobi<SparseMatrix<double>>;
1418 *   auto smoother =
1419 *   std::make_unique<MGSmootherPrecondition<SparseMatrix<double>,
1420 *   Smoother,
1421 *   Vector<double>>>();
1422 *   smoother->initialize(mg_matrices,
1423 *   Smoother::AdditionalData(fe.degree == 1 ? 0.6667 :
1424 *   0.47));
1425 *   smoother->set_steps(settings.smoothing_steps);
1426 *   mg_smoother = std::move(smoother);
1427 *   }
1428 *   else if (settings.smoother_type == "block SOR" ||
1429 *   settings.smoother_type == "block Jacobi")
1430 *   {
1431 *   smoother_data.resize(0, triangulation.n_levels() - 1);
1432 *  
1433 *   for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
1434 *   {
1435 *   DoFTools::make_cell_patches(smoother_data[level].block_list,
1436 *   dof_handler,
1437 *   level);
1438 *  
1439 *   smoother_data[level].relaxation =
1440 *   (settings.smoother_type == "block SOR" ? 1.0 : 0.25);
1441 *   smoother_data[level].inversion = PreconditionBlockBase<double>::svd;
1442 *  
1443 *   std::vector<unsigned int> ordered_indices;
1444 *   switch (settings.dof_renumbering)
1445 *   {
1446 *   case Settings::DoFRenumberingStrategy::downstream:
1447 *   ordered_indices =
1448 *   create_downstream_cell_ordering(dof_handler,
1449 *   advection_direction,
1450 *   level);
1451 *   break;
1452 *  
1453 *   case Settings::DoFRenumberingStrategy::upstream:
1454 *   ordered_indices =
1455 *   create_downstream_cell_ordering(dof_handler,
1456 *   -1.0 * advection_direction,
1457 *   level);
1458 *   break;
1459 *  
1460 *   case Settings::DoFRenumberingStrategy::random:
1461 *   ordered_indices =
1462 *   create_random_cell_ordering(dof_handler, level);
1463 *   break;
1464 *  
1465 *   case Settings::DoFRenumberingStrategy::none:
1466 *   break;
1467 *  
1468 *   default:
1469 *   AssertThrow(false, ExcNotImplemented());
1470 *   break;
1471 *   }
1472 *  
1473 *   smoother_data[level].order =
1474 *   std::vector<std::vector<unsigned int>>(1, ordered_indices);
1475 *   }
1476 *  
1477 *   if (settings.smoother_type == "block SOR")
1478 *   {
1479 *   auto smoother = std::make_unique<MGSmootherPrecondition<
1480 *   SparseMatrix<double>,
1481 *   RelaxationBlockSOR<SparseMatrix<double>, double, Vector<double>>,
1482 *   Vector<double>>>();
1483 *   smoother->initialize(mg_matrices, smoother_data);
1484 *   smoother->set_steps(settings.smoothing_steps);
1485 *   mg_smoother = std::move(smoother);
1486 *   }
1487 *   else if (settings.smoother_type == "block Jacobi")
1488 *   {
1489 *   auto smoother = std::make_unique<
1490 *   MGSmootherPrecondition<SparseMatrix<double>,
1491 *   RelaxationBlockJacobi<SparseMatrix<double>,
1492 *   double,
1493 *   Vector<double>>,
1494 *   Vector<double>>>();
1495 *   smoother->initialize(mg_matrices, smoother_data);
1496 *   smoother->set_steps(settings.smoothing_steps);
1497 *   mg_smoother = std::move(smoother);
1498 *   }
1499 *   }
1500 *   else
1501 *   AssertThrow(false, ExcNotImplemented());
1502 *   }
1503 *  
1504 *  
1505 * @endcode
1506 *
1507 *
1508 * <a name="codeAdvectionProblemsolvecode"></a>
1509 * <h4><code>AdvectionProblem::solve()</code></h4>
1510 *
1511
1512 *
1513 * Before we can solve the system, we must first set up the multigrid
1514 * preconditioner. This requires the setup of the transfer between levels,
1515 * the coarse matrix solver, and the smoother. This setup follows almost
1516 * identically to @ref step_16 "step-16", the main difference being the various smoothers
1517 * defined above and the fact that we need different interface edge matrices
1518 * for in and out since our problem is non-symmetric. (In reality, for this
1519 * tutorial these interface matrices are empty since we are only using global
1520 * refinement, and thus have no refinement edges. However, we have still
1521 * included both here since if one made the simple switch to an adaptively
1522 * refined method, the program would still run correctly.)
1523 *
1524
1525 *
1526 * The last thing to note is that since our problem is non-symmetric, we must
1527 * use an appropriate Krylov subspace method. We choose here to
1528 * use GMRES since it offers the guarantee of residual reduction in each
1529 * iteration. The major disadvantage of GMRES is that, for each iteration,
1530 * the number of stored temporary vectors increases by one, and one also needs
1531 * to compute a scalar product with all previously stored vectors. This is
1532 * rather expensive. This requirement is relaxed by using the restarted GMRES
1533 * method which puts a cap on the number of vectors we are required to store
1534 * at any one time (here we restart after 50 temporary vectors, or 48
1535 * iterations). This then has the disadvantage that we lose information we
1536 * have gathered throughout the iteration and therefore we could see slower
1537 * convergence. As a consequence, where to restart is a question of balancing
1538 * memory consumption, CPU effort, and convergence speed.
1539 * However, the goal of this tutorial is to have very low
1540 * iteration counts by using a powerful GMG preconditioner, so we have picked
1541 * the restart length such that all of the results shown below converge prior
1542 * to restart happening, and thus we have a standard GMRES method. If the user
1543 * is interested, another suitable method offered in deal.II would be
1544 * BiCGStab.
1545 *
1546
1547 *
1548 *
1549 * @code
1550 *   template <int dim>
1551 *   void AdvectionProblem<dim>::solve()
1552 *   {
1553 *   const unsigned int max_iters = 200;
1554 *   const double solve_tolerance = 1e-8 * system_rhs.l2_norm();
1555 *   SolverControl solver_control(max_iters, solve_tolerance, true, true);
1556 *   solver_control.enable_history_data();
1557 *  
1558 *   using Transfer = MGTransferPrebuilt<Vector<double>>;
1559 *   Transfer mg_transfer(mg_constrained_dofs);
1560 *   mg_transfer.build(dof_handler);
1561 *  
1562 *   FullMatrix<double> coarse_matrix;
1563 *   coarse_matrix.copy_from(mg_matrices[0]);
1564 *   MGCoarseGridHouseholder<double, Vector<double>> coarse_grid_solver;
1565 *   coarse_grid_solver.initialize(coarse_matrix);
1566 *  
1567 *   setup_smoother();
1568 *  
1569 *   mg_matrix.initialize(mg_matrices);
1570 *   mg_interface_matrix_in.initialize(mg_interface_in);
1571 *   mg_interface_matrix_out.initialize(mg_interface_out);
1572 *  
1573 *   Multigrid<Vector<double>> mg(
1574 *   mg_matrix, coarse_grid_solver, mg_transfer, *mg_smoother, *mg_smoother);
1575 *   mg.set_edge_matrices(mg_interface_matrix_out, mg_interface_matrix_in);
1576 *  
1577 *   PreconditionMG<dim, Vector<double>, Transfer> preconditioner(dof_handler,
1578 *   mg,
1579 *   mg_transfer);
1580 *  
1581 *   std::cout << " Solving with GMRES to tol " << solve_tolerance << "..."
1582 *   << std::endl;
1583 *   SolverGMRES<Vector<double>> solver(
1584 *   solver_control, SolverGMRES<Vector<double>>::AdditionalData(50, true));
1585 *  
1586 *   Timer time;
1587 *   time.start();
1588 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
1589 *   time.stop();
1590 *  
1591 *   std::cout << " converged in " << solver_control.last_step()
1592 *   << " iterations"
1593 *   << " in " << time.last_wall_time() << " seconds " << std::endl;
1594 *  
1595 *   constraints.distribute(solution);
1596 *  
1597 *   mg_smoother.release();
1598 *   }
1599 *  
1600 *  
1601 * @endcode
1602 *
1603 *
1604 * <a name="codeAdvectionProblemoutput_resultscode"></a>
1605 * <h4><code>AdvectionProblem::output_results()</code></h4>
1606 *
1607
1608 *
1609 * The final function of interest generates graphical output.
1610 * Here we output the solution and cell ordering in a .vtu format.
1611 *
1612
1613 *
1614 * At the top of the function, we generate an index for each cell to
1615 * visualize the ordering used by the smoothers. Note that we do
1616 * this only for the active cells instead of the levels, where the
1617 * smoothers are actually used. For the point smoothers we renumber
1618 * DoFs instead of cells, so this is only an approximation of what
1619 * happens in reality. Finally, the random ordering is not the
1620 * random ordering we actually use (see `create_smoother()` for that).
1621 *
1622
1623 *
1624 * The (integer) ordering of cells is then copied into a (floating
1625 * point) vector for graphical output.
1626 *
1627 * @code
1628 *   template <int dim>
1629 *   void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
1630 *   {
1631 *   const unsigned int n_active_cells = triangulation.n_active_cells();
1632 *   Vector<double> cell_indices(n_active_cells);
1633 *   {
1634 *   std::vector<unsigned int> ordered_indices;
1635 *   switch (settings.dof_renumbering)
1636 *   {
1637 *   case Settings::DoFRenumberingStrategy::downstream:
1638 *   ordered_indices =
1639 *   create_downstream_cell_ordering(dof_handler, advection_direction);
1640 *   break;
1641 *  
1642 *   case Settings::DoFRenumberingStrategy::upstream:
1643 *   ordered_indices =
1644 *   create_downstream_cell_ordering(dof_handler,
1645 *   -1.0 * advection_direction);
1646 *   break;
1647 *  
1648 *   case Settings::DoFRenumberingStrategy::random:
1649 *   ordered_indices = create_random_cell_ordering(dof_handler);
1650 *   break;
1651 *  
1652 *   case Settings::DoFRenumberingStrategy::none:
1653 *   ordered_indices.resize(n_active_cells);
1654 *   for (unsigned int i = 0; i < n_active_cells; ++i)
1655 *   ordered_indices[i] = i;
1656 *   break;
1657 *  
1658 *   default:
1659 *   AssertThrow(false, ExcNotImplemented());
1660 *   break;
1661 *   }
1662 *  
1663 *   for (unsigned int i = 0; i < n_active_cells; ++i)
1664 *   cell_indices(ordered_indices[i]) = static_cast<double>(i);
1665 *   }
1666 *  
1667 * @endcode
1668 *
1669 * The remainder of the function is then straightforward, given
1670 * previous tutorial programs:
1671 *
1672 * @code
1673 *   DataOut<dim> data_out;
1674 *   data_out.attach_dof_handler(dof_handler);
1675 *   data_out.add_data_vector(solution, "solution");
1676 *   data_out.add_data_vector(cell_indices, "cell_index");
1677 *   data_out.build_patches();
1678 *  
1679 *   const std::string filename =
1680 *   "solution-" + Utilities::int_to_string(cycle) + ".vtu";
1681 *   std::ofstream output(filename.c_str());
1682 *   data_out.write_vtu(output);
1683 *   }
1684 *  
1685 *  
1686 * @endcode
1687 *
1688 *
1689 * <a name="codeAdvectionProblemruncode"></a>
1690 * <h4><code>AdvectionProblem::run()</code></h4>
1691 *
1692
1693 *
1694 * As in most tutorials, this function creates/refines the mesh and calls
1695 * the various functions defined above to set up, assemble, solve, and output
1696 * the results.
1697 *
1698
1699 *
1700 * In cycle zero, we generate the mesh for the on the square
1701 * <code>[-1,1]^dim</code> with a hole of radius 3/10 units centered
1702 * at the origin. For objects with `manifold_id` equal to one
1703 * (namely, the faces adjacent to the hole), we assign a spherical
1704 * manifold.
1705 *
1706
1707 *
1708 *
1709 * @code
1710 *   template <int dim>
1711 *   void AdvectionProblem<dim>::run()
1712 *   {
1713 *   for (unsigned int cycle = 0; cycle < (settings.fe_degree == 1 ? 7 : 5);
1714 *   ++cycle)
1715 *   {
1716 *   std::cout << " Cycle " << cycle << ':' << std::endl;
1717 *  
1718 *   if (cycle == 0)
1719 *   {
1720 *   GridGenerator::hyper_cube_with_cylindrical_hole(triangulation,
1721 *   0.3,
1722 *   1.0);
1723 *  
1724 *   const SphericalManifold<dim> manifold_description(Point<dim>(0, 0));
1725 *   triangulation.set_manifold(1, manifold_description);
1726 *   }
1727 *  
1728 *   triangulation.refine_global();
1729 *  
1730 *   setup_system();
1731 *  
1732 *   std::cout << " Number of active cells: "
1733 *   << triangulation.n_active_cells() << " ("
1734 *   << triangulation.n_levels() << " levels)" << std::endl;
1735 *   std::cout << " Number of degrees of freedom: "
1736 *   << dof_handler.n_dofs() << std::endl;
1737 *  
1738 *   assemble_system_and_multigrid();
1739 *  
1740 *   solve();
1741 *  
1742 *   if (settings.output)
1743 *   output_results(cycle);
1744 *  
1745 *   std::cout << std::endl;
1746 *   }
1747 *   }
1748 *   } // namespace Step63
1749 *  
1750 *  
1751 * @endcode
1752 *
1753 *
1754 * <a name="Thecodemaincodefunction"></a>
1755 * <h3>The <code>main</code> function</h3>
1756 *
1757
1758 *
1759 * Finally, the main function is like most tutorials. The only
1760 * interesting bit is that we require the user to pass a `.prm` file
1761 * as a sole command line argument. If no parameter file is given, the
1762 * program will output the contents of a sample parameter file with
1763 * all default values to the screen that the user can then copy and
1764 * paste into their own `.prm` file.
1765 *
1766
1767 *
1768 *
1769 * @code
1770 *   int main(int argc, char *argv[])
1771 *   {
1772 *   try
1773 *   {
1774 *   Step63::Settings settings;
1775 *   settings.get_parameters((argc > 1) ? (argv[1]) : "");
1776 *  
1777 *   Step63::AdvectionProblem<2> advection_problem_2d(settings);
1778 *   advection_problem_2d.run();
1779 *   }
1780 *   catch (std::exception &exc)
1781 *   {
1782 *   std::cerr << std::endl
1783 *   << std::endl
1784 *   << "----------------------------------------------------"
1785 *   << std::endl;
1786 *   std::cerr << "Exception on processing: " << std::endl
1787 *   << exc.what() << std::endl
1788 *   << "Aborting!" << std::endl
1789 *   << "----------------------------------------------------"
1790 *   << std::endl;
1791 *   return 1;
1792 *   }
1793 *   catch (...)
1794 *   {
1795 *   std::cerr << std::endl
1796 *   << std::endl
1797 *   << "----------------------------------------------------"
1798 *   << std::endl;
1799 *   std::cerr << "Unknown exception!" << std::endl
1800 *   << "Aborting!" << std::endl
1801 *   << "----------------------------------------------------"
1802 *   << std::endl;
1803 *   return 1;
1804 *   }
1805 *  
1806 *   return 0;
1807 *   }
1808 * @endcode
1809<a name="Results"></a><h1>Results</h1>
1810
1811
1812<a name="GMRESIterationNumbers"></a><h3> GMRES Iteration Numbers </h3>
1813
1814
1815The major advantage for GMG is that it is an @f$\mathcal{O}(n)@f$ method,
1816that is, the complexity of the problem increases linearly with the
1817problem size. To show then that the linear solver presented in this
1818tutorial is in fact @f$\mathcal{O}(n)@f$, all one needs to do is show that
1819the iteration counts for the GMRES solve stay roughly constant as we
1820refine the mesh.
1821
1822Each of the following tables gives the GMRES iteration counts to reduce the
1823initial residual by a factor of @f$10^8@f$. We selected a sufficient number of smoothing steps
1824(based on the method) to get iteration numbers independent of mesh size. As
1825can be seen from the tables below, the method is indeed @f$\mathcal{O}(n)@f$.
1826
1827<a name="DoFCellRenumbering"></a><h4> DoF/Cell Renumbering </h4>
1828
1829
1830The point-wise smoothers ("Jacobi" and "SOR") get applied in the order the
1831DoFs are numbered on each level. We can influence this using the
1832DoFRenumbering namespace. The block smoothers are applied based on the
1833ordering we set in `setup_smoother()`. We can visualize this numbering. The
1834following pictures show the cell numbering of the active cells in downstream,
1835random, and upstream numbering (left to right):
1836
1837<img src="https://www.dealii.org/images/steps/developer/step-63-cell-order.png" alt="">
1838
1839Let us start with the additive smoothers. The following table shows
1840the number of iterations necessary to obtain convergence from GMRES:
1841
1842<table align="center" class="doxtable">
1843<tr>
1844 <th></th>
1845 <th></th>
1846 <th colspan="1">@f$Q_1@f$</th>
1847 <th colspan="7">Smoother (smoothing steps)</th>
1848</tr>
1849<tr>
1850 <th></th>
1851 <th></th>
1852 <th></th>
1853 <th colspan="3">Jacobi (6)</th>
1854 <th></th>
1855 <th colspan="3">Block Jacobi (3)</th>
1856</tr>
1857<tr>
1858 <th></th>
1859 <th></th>
1860 <th></th>
1861 <th colspan="3">Renumbering Strategy</th>
1862 <th></th>
1863 <th colspan="3">Renumbering Strategy</th>
1864</tr>
1865<tr>
1866 <th>Cells</th>
1867 <th></th>
1868 <th>DoFs</th>
1869 <th>Downstream</th>
1870 <th>Random</th>
1871 <th>Upstream</th>
1872 <th></th>
1873 <th>Downstream</th>
1874 <th>Random</th>
1875 <th>Upstream</th>
1876</tr>
1877<tr>
1878 <th>32</th>
1879 <th></th>
1880 <th>48</th>
1881 <td>3</th>
1882 <td>3</th>
1883 <td>3</th>
1884 <th></th>
1885 <td>3</th>
1886 <td>3</th>
1887 <td>3</th>
1888</tr>
1889<tr>
1890 <th>128</th>
1891 <th></th>
1892 <th>160</th>
1893 <td>6</th>
1894 <td>6</th>
1895 <td>6</th>
1896 <th></th>
1897 <td>6</th>
1898 <td>6</th>
1899 <td>6</th>
1900</tr>
1901<tr>
1902 <th>512</th>
1903 <th></th>
1904 <th>576</th>
1905 <td>11</th>
1906 <td>11</th>
1907 <td>11</th>
1908 <th></th>
1909 <td>9</th>
1910 <td>9</th>
1911 <td>9</th>
1912</tr>
1913<tr>
1914 <th>2048</th>
1915 <th></th>
1916 <th>2176</th>
1917 <td>15</th>
1918 <td>15</th>
1919 <td>15</th>
1920 <th></th>
1921 <td>13</th>
1922 <td>13</th>
1923 <td>13</th>
1924</tr>
1925<tr>
1926 <th>8192</th>
1927 <th></th>
1928 <th>8448</th>
1929 <td>18</th>
1930 <td>18</th>
1931 <td>18</th>
1932 <th></th>
1933 <td>15</th>
1934 <td>15</th>
1935 <td>15</th>
1936</tr>
1937<tr>
1938 <th>32768</th>
1939 <th></th>
1940 <th>33280</th>
1941 <td>20</th>
1942 <td>20</th>
1943 <td>20</th>
1944 <th></th>
1945 <td>16</th>
1946 <td>16</th>
1947 <td>16</th>
1948</tr>
1949<tr>
1950 <th>131072</th>
1951 <th></th>
1952 <th>132096</th>
1953 <td>20</th>
1954 <td>20</th>
1955 <td>20</th>
1956 <th></th>
1957 <td>16</th>
1958 <td>16</th>
1959 <td>16</th>
1960</tr>
1961</table>
1962
1963We see that renumbering the
1964DoFs/cells has no effect on convergence speed. This is because these
1965smoothers compute operations on each DoF (point-smoother) or cell
1966(block-smoother) independently and add up the results. Since we can
1967define these smoothers as an application of a sum of matrices, and
1968matrix addition is commutative, the order at which we sum the
1969different components will not affect the end result.
1970
1971On the other hand, the situation is different for multiplicative smoothers:
1972
1973<table align="center" class="doxtable">
1974<tr>
1975 <th></th>
1976 <th></th>
1977 <th colspan="1">@f$Q_1@f$</th>
1978 <th colspan="7">Smoother (smoothing steps)</th>
1979</tr>
1980<tr>
1981 <th></th>
1982 <th></th>
1983 <th></th>
1984 <th colspan="3">SOR (3)</th>
1985 <th></th>
1986 <th colspan="3">Block SOR (1)</th>
1987</tr>
1988<tr>
1989 <th></th>
1990 <th></th>
1991 <th></th>
1992 <th colspan="3">Renumbering Strategy</th>
1993 <th></th>
1994 <th colspan="3">Renumbering Strategy</th>
1995</tr>
1996<tr>
1997 <th>Cells</th>
1998 <th></th>
1999 <th>DoFs</th>
2000 <th>Downstream</th>
2001 <th>Random</th>
2002 <th>Upstream</th>
2003 <th></th>
2004 <th>Downstream</th>
2005 <th>Random</th>
2006 <th>Upstream</th>
2007</tr>
2008<tr>
2009 <th>32</th>
2010 <th></th>
2011 <th>48</th>
2012 <td>2</th>
2013 <td>2</th>
2014 <td>3</th>
2015 <th></th>
2016 <td>2</th>
2017 <td>2</th>
2018 <td>3</th>
2019</tr>
2020<tr>
2021 <th>128</th>
2022 <th></th>
2023 <th>160</th>
2024 <td>5</th>
2025 <td>5</th>
2026 <td>7</th>
2027 <th></th>
2028 <td>5</th>
2029 <td>5</th>
2030 <td>7</th>
2031</tr>
2032<tr>
2033 <th>512</th>
2034 <th></th>
2035 <th>576</th>
2036 <td>7</th>
2037 <td>9</th>
2038 <td>11</th>
2039 <th></th>
2040 <td>7</th>
2041 <td>7</th>
2042 <td>12</th>
2043</tr>
2044<tr>
2045 <th>2048</th>
2046 <th></th>
2047 <th>2176</th>
2048 <td>10</th>
2049 <td>12</th>
2050 <td>15</th>
2051 <th></th>
2052 <td>8</th>
2053 <td>10</th>
2054 <td>17</th>
2055</tr>
2056<tr>
2057 <th>8192</th>
2058 <th></th>
2059 <th>8448</th>
2060 <td>11</th>
2061 <td>15</th>
2062 <td>19</th>
2063 <th></th>
2064 <td>10</th>
2065 <td>11</th>
2066 <td>20</th>
2067</tr>
2068<tr>
2069 <th>32768</th>
2070 <th></th>
2071 <th>33280</th>
2072 <td>12</th>
2073 <td>16</th>
2074 <td>20</th>
2075 <th></th>
2076 <td>10</th>
2077 <td>12</th>
2078 <td>21</th>
2079</tr>
2080<tr>
2081 <th>131072</th>
2082 <th></th>
2083 <th>132096</th>
2084 <td>12</th>
2085 <td>16</th>
2086 <td>19</th>
2087 <th></th>
2088 <td>11</th>
2089 <td>12</th>
2090 <td>21</th>
2091</tr>
2092</table>
2093
2094Here, we can speed up
2095convergence by renumbering the DoFs/cells in the advection direction,
2096and similarly, we can slow down convergence if we do the renumbering
2097in the opposite direction. This is because advection-dominated
2098problems have a directional flow of information (in the advection
2099direction) which, given the right renumbering of DoFs/cells,
2100multiplicative methods are able to capture.
2101
2102This feature of multiplicative methods is, however, dependent on the
2103value of @f$\varepsilon@f$. As we increase @f$\varepsilon@f$ and the problem
2104becomes more diffusion-dominated, we have a more uniform propagation
2105of information over the mesh and there is a diminished advantage for
2106renumbering in the advection direction. On the opposite end, in the
2107extreme case of @f$\varepsilon=0@f$ (advection-only), we have a 1st-order
2108PDE and multiplicative methods with the right renumbering become
2109effective solvers: A correct downstream numbering may lead to methods
2110that require only a single iteration because information can be
2111propagated from the inflow boundary downstream, with no information
2112transport in the opposite direction. (Note, however, that in the case
2113of @f$\varepsilon=0@f$, special care must be taken for the boundary
2114conditions in this case).
2115
2116
2117<a name="Pointvsblocksmoothers"></a><h4> %Point vs. block smoothers </h4>
2118
2119
2120We will limit the results to runs using the downstream
2121renumbering. Here is a cross comparison of all four smoothers for both
2122@f$Q_1@f$ and @f$Q_3@f$ elements:
2123
2124<table align="center" class="doxtable">
2125<tr>
2126 <th></th>
2127 <td></th>
2128 <th colspan="1">@f$Q_1@f$</th>
2129 <th colspan="4">Smoother (smoothing steps)</th>
2130 <th></th>
2131 <th colspan="1">@f$Q_3@f$</th>
2132 <th colspan="4">Smoother (smoothing steps)</th>
2133</tr>
2134<tr>
2135 <th colspan="1">Cells</th>
2136 <td></th>
2137 <th colspan="1">DoFs</th>
2138 <th colspan="1">Jacobi (6)</th>
2139 <th colspan="1">Block Jacobi (3)</th>
2140 <th colspan="1">SOR (3)</th>
2141 <th colspan="1">Block SOR (1)</th>
2142 <th></th>
2143 <th colspan="1">DoFs</th>
2144 <th colspan="1">Jacobi (6)</th>
2145 <th colspan="1">Block Jacobi (3)</th>
2146 <th colspan="1">SOR (3)</th>
2147 <th colspan="1">Block SOR (1)</th>
2148</tr>
2149<tr>
2150 <th>32</th>
2151 <td></th>
2152 <th>48</th>
2153 <td>3</th>
2154 <td>3</th>
2155 <td>2</th>
2156 <td>2</th>
2157 <td></th>
2158 <th>336</th>
2159 <td>15</th>
2160 <td>14</th>
2161 <td>15</th>
2162 <td>6</th>
2163</tr>
2164<tr>
2165 <th>128</th>
2166 <td></th>
2167 <th>160</th>
2168 <td>6</th>
2169 <td>6</th>
2170 <td>5</th>
2171 <td>5</th>
2172 <td></th>
2173 <th>1248</th>
2174 <td>23</th>
2175 <td>18</th>
2176 <td>21</th>
2177 <td>9</th>
2178</tr>
2179<tr>
2180 <th>512</th>
2181 <td></th>
2182 <th>576</th>
2183 <td>11</th>
2184 <td>9</th>
2185 <td>7</th>
2186 <td>7</th>
2187 <td></th>
2188 <th>4800</th>
2189 <td>29</th>
2190 <td>21</th>
2191 <td>28</th>
2192 <td>9</th>
2193</tr>
2194<tr>
2195 <th>2048</th>
2196 <td></th>
2197 <th>2176</th>
2198 <td>15</th>
2199 <td>13</th>
2200 <td>10</th>
2201 <td>8</th>
2202 <td></th>
2203 <th>18816</th>
2204 <td>33</th>
2205 <td>22</th>
2206 <td>32</th>
2207 <td>9</th>
2208</tr>
2209<tr>
2210 <th>8192</th>
2211 <td></th>
2212 <th>8448</th>
2213 <td>18</th>
2214 <td>15</th>
2215 <td>11</th>
2216 <td>10</th>
2217 <td></th>
2218 <th>74496</th>
2219 <td>35</th>
2220 <td>22</th>
2221 <td>34</th>
2222 <td>10</th>
2223</tr>
2224<tr>
2225 <th>32768</th>
2226 <td></th>
2227 <th>33280</th>
2228 <td>20</th>
2229 <td>16</th>
2230 <td>12</th>
2231 <td>10</th>
2232 <td></th>
2233</tr>
2234<tr>
2235 <th>131072</th>
2236 <td></th>
2237 <th>132096</th>
2238 <td>20</th>
2239 <td>16</th>
2240 <td>12</th>
2241 <td>11</th>
2242 <td></th>
2243</tr>
2244</table>
2245
2246We see that for @f$Q_1@f$, both multiplicative smoothers require a smaller
2247combination of smoothing steps and iteration counts than either
2248additive smoother. However, when we increase the degree to a @f$Q_3@f$
2249element, there is a clear advantage for the block smoothers in terms
2250of the number of smoothing steps and iterations required to
2251solve. Specifically, the block SOR smoother gives constant iteration
2252counts over the degree, and the block Jacobi smoother only sees about
2253a 38% increase in iterations compared to 75% and 183% for Jacobi and
2254SOR respectively.
2255
2256<a name="Cost"></a><h3> Cost </h3>
2257
2258
2259Iteration counts do not tell the full story in the optimality of a one
2260smoother over another. Obviously we must examine the cost of an
2261iteration. Block smoothers here are at a disadvantage as they are
2262having to construct and invert a cell matrix for each cell. Here is a
2263comparison of solve times for a @f$Q_3@f$ element with 74,496 DoFs:
2264
2265<table align="center" class="doxtable">
2266<tr>
2267 <th colspan="1">@f$Q_3@f$</th>
2268 <th colspan="4">Smoother (smoothing steps)</th>
2269</tr>
2270<tr>
2271 <th colspan="1">DoFs</th>
2272 <th colspan="1">Jacobi (6)</th>
2273 <th colspan="1">Block Jacobi (3)</th>
2274 <th colspan="1">SOR (3)</th>
2275 <th colspan="1">Block SOR (1)</th>
2276</tr>
2277<tr>
2278 <th>74496</th>
2279 <td>0.68s</th>
2280 <td>5.82s</th>
2281 <td>1.18s</th>
2282 <td>1.02s</th>
2283</tr>
2284</table>
2285
2286The smoother that requires the most iterations (Jacobi) actually takes
2287the shortest time (roughly 2/3 the time of the next fastest
2288method). This is because all that is required to apply a Jacobi
2289smoothing step is multiplication by a diagonal matrix which is very
2290cheap. On the other hand, while SOR requires over 3x more iterations
2291(each with 3x more smoothing steps) than block SOR, the times are
2292roughly equivalent, implying that a smoothing step of block SOR is
2293roughly 9x slower than a smoothing step of SOR. Lastly, block Jacobi
2294is almost 6x more expensive than block SOR, which intuitively makes
2295sense from the fact that 1 step of each method has the same cost
2296(inverting the cell matrices and either adding or multiply them
2297together), and block Jacobi has 3 times the number of smoothing steps per
2298iteration with 2 times the iterations.
2299
2300
2301<a name="Additionalpoints"></a><h3> Additional points </h3>
2302
2303
2304There are a few more important points to mention:
2305
2306<ol>
2307<li> For a mesh distributed in parallel, multiplicative methods cannot
2308be executed over the entire domain. This is because they operate one
2309cell at a time, and downstream cells can only be handled once upstream
2310cells have already been done. This is fine on a single processor: The
2311processor just goes through the list of cells one after the
2312other. However, in parallel, it would imply that some processors are
2313idle because upstream processors have not finished doing the work on
2314cells upstream from the ones owned by the current processor. Once the
2315upstream processors are done, the downstream ones can start, but by
2316that time the upstream processors have no work left. In other words,
2317most of the time during these smoother steps, most processors are in
2318fact idle. This is not how one obtains good parallel scalability!
2319
2320One can use a hybrid method where
2321a multiplicative smoother is applied on each subdomain, but as you
2322increase the number of subdomains, the method approaches the behavior
2323of an additive method. This is a major disadvantage to these methods.
2324</li>
2325
2326<li> Current research into block smoothers suggest that soon we will be
2327able to compute the inverse of the cell matrices much cheaper than
2328what is currently being done inside deal.II. This research is based on
2329the fast diagonalization method (dating back to the 1960s) and has
2330been used in the spectral community for around 20 years (see, e.g., <a
2331href="https://doi.org/10.1007/s10915-004-4787-3"> Hybrid
2332Multigrid/Schwarz Algorithms for the Spectral Element Method by Lottes
2333and Fischer</a>). There are currently efforts to generalize these
2334methods to DG and make them more robust. Also, it seems that one
2335should be able to take advantage of matrix-free implementations and
2336the fact that, in the interior of the domain, cell matrices tend to
2337look very similar, allowing fewer matrix inverse computations.
2338</li>
2339</ol>
2340
2341Combining 1. and 2. gives a good reason for expecting that a method
2342like block Jacobi could become very powerful in the future, even
2343though currently for these examples it is quite slow.
2344
2345
2346<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
2347
2348
2349<a name="ConstantiterationsforQsub5sub"></a><h4> Constant iterations for Q<sub>5</sub> </h4>
2350
2351
2352Change the number of smoothing steps and the smoother relaxation
2353parameter (set in <code>Smoother::AdditionalData()</code> inside
2354<code>create_smoother()</code>, only necessary for point smoothers) so
2355that we maintain a constant number of iterations for a @f$Q_5@f$ element.
2356
2357<a name="Effectivenessofrenumberingforchangingepsilon"></a><h4> Effectiveness of renumbering for changing epsilon </h4>
2358
2359
2360Increase/decrease the parameter "Epsilon" in the `.prm` files of the
2361multiplicative methods and observe for which values renumbering no
2362longer influences convergence speed.
2363
2364<a name="Meshadaptivity"></a><h4> Mesh adaptivity </h4>
2365
2366
2367The code is set up to work correctly with an adaptively refined mesh (the
2368interface matrices are created and set). Devise a suitable refinement
2369criterium or try the KellyErrorEstimator class.
2370 *
2371 *
2372<a name="PlainProg"></a>
2373<h1> The plain program</h1>
2374@include "step-63.cc"
2375*/
Definition fe_q.h:551
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< RangeNumberType > &values, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
Definition point.h:112
Point< 3 > center
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:282
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Expression coth(const Expression &x)
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void random(DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true)
Definition mg_tools.cc:576
void make_interface_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, const MGConstrainedDoFs &mg_constrained_dofs, SparsityPatternBase &sparsity, const unsigned int level)
Definition mg_tools.cc:1014
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
Definition numbers.h:259
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)
const TriangulationDescription::Settings settings