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-50.h
Go to the documentation of this file.
1,
481 *   const unsigned int /*component*/ = 0) const override
482 *   {
483 *   return 1.0;
484 *   }
485 *  
486 *  
487 *   template <typename number>
489 *   value(const Point<dim, VectorizedArray<number>> & /*p*/,
490 *   const unsigned int /*component*/ = 0) const
491 *   {
492 *   return VectorizedArray<number>(1.0);
493 *   }
494 *   };
495 *  
496 *  
497 * @endcode
498 *
499 * This next class represents the diffusion coefficient. We use a variable
500 * coefficient which is 100.0 at any point where at least one coordinate is
501 * less than -0.5, and 1.0 at all other points. As above, a separate value()
502 * returning a VectorizedArray is used for the matrix-free code. An @p
503 * average() function computes the arithmetic average for a set of points.
504 *
505 * @code
506 *   template <int dim>
507 *   class Coefficient : public Function<dim>
508 *   {
509 *   public:
510 *   virtual double value(const Point<dim> &p,
511 *   const unsigned int /*component*/ = 0) const override;
512 *  
513 *   template <typename number>
515 *   const unsigned int /*component*/ = 0) const;
516 *  
517 *   template <typename number>
518 *   number average_value(const std::vector<Point<dim, number>> &points) const;
519 *  
520 * @endcode
521 *
522 * When using a coefficient in the MatrixFree framework, we also
523 * need a function that creates a Table of coefficient values for a
524 * set of cells provided by the MatrixFree operator argument here.
525 *
526 * @code
527 *   template <typename number>
528 *   std::shared_ptr<Table<2, VectorizedArray<number>>> make_coefficient_table(
529 *   const MatrixFree<dim, number, VectorizedArray<number>> &mf_storage) const;
530 *   };
531 *  
532 *  
533 *  
534 *   template <int dim>
535 *   double Coefficient<dim>::value(const Point<dim> &p, const unsigned int) const
536 *   {
537 *   for (int d = 0; d < dim; ++d)
538 *   {
539 *   if (p[d] < -0.5)
540 *   return 100.0;
541 *   }
542 *   return 1.0;
543 *   }
544 *  
545 *  
546 *  
547 *   template <int dim>
548 *   template <typename number>
550 *   Coefficient<dim>::value(const Point<dim, VectorizedArray<number>> &p,
551 *   const unsigned int) const
552 *   {
553 *   VectorizedArray<number> return_value = VectorizedArray<number>(1.0);
554 *   for (unsigned int i = 0; i < VectorizedArray<number>::size(); ++i)
555 *   {
556 *   for (int d = 0; d < dim; ++d)
557 *   if (p[d][i] < -0.5)
558 *   {
559 *   return_value[i] = 100.0;
560 *   break;
561 *   }
562 *   }
563 *  
564 *   return return_value;
565 *   }
566 *  
567 *  
568 *  
569 *   template <int dim>
570 *   template <typename number>
571 *   number Coefficient<dim>::average_value(
572 *   const std::vector<Point<dim, number>> &points) const
573 *   {
574 *   number average(0);
575 *   for (unsigned int i = 0; i < points.size(); ++i)
576 *   average += value(points[i]);
577 *   average /= points.size();
578 *  
579 *   return average;
580 *   }
581 *  
582 *  
583 *  
584 *   template <int dim>
585 *   template <typename number>
586 *   std::shared_ptr<Table<2, VectorizedArray<number>>>
587 *   Coefficient<dim>::make_coefficient_table(
588 *   const MatrixFree<dim, number, VectorizedArray<number>> &mf_storage) const
589 *   {
590 *   auto coefficient_table =
591 *   std::make_shared<Table<2, VectorizedArray<number>>>();
592 *  
593 *   FEEvaluation<dim, -1, 0, 1, number> fe_eval(mf_storage);
594 *  
595 *   const unsigned int n_cells = mf_storage.n_cell_batches();
596 *   const unsigned int n_q_points = fe_eval.n_q_points;
597 *  
598 *   coefficient_table->reinit(n_cells, 1);
599 *  
600 *   for (unsigned int cell = 0; cell < n_cells; ++cell)
601 *   {
602 *   fe_eval.reinit(cell);
603 *  
604 *   VectorizedArray<number> average_value = 0.;
605 *   for (unsigned int q = 0; q < n_q_points; ++q)
606 *   average_value += value(fe_eval.quadrature_point(q));
607 *   average_value /= n_q_points;
608 *  
609 *   (*coefficient_table)(cell, 0) = average_value;
610 *   }
611 *  
612 *   return coefficient_table;
613 *   }
614 *  
615 *  
616 *  
617 * @endcode
618 *
619 *
620 * <a name="Runtimeparameters"></a>
621 * <h3>Run time parameters</h3>
622 *
623
624 *
625 * We will use ParameterHandler to pass in parameters at runtime. The
626 * structure @p Settings parses and stores these parameters to be queried
627 * throughout the program.
628 *
629 * @code
630 *   struct Settings
631 *   {
632 *   bool try_parse(const std::string &prm_filename);
633 *  
634 *   enum SolverType
635 *   {
636 *   gmg_mb,
637 *   gmg_mf,
638 *   amg
639 *   };
640 *  
641 *   SolverType solver;
642 *  
643 *   int dimension;
644 *   double smoother_dampen;
645 *   unsigned int smoother_steps;
646 *   unsigned int n_steps;
647 *   bool output;
648 *   };
649 *  
650 *  
651 *  
652 *   bool Settings::try_parse(const std::string &prm_filename)
653 *   {
654 *   ParameterHandler prm;
655 *   prm.declare_entry("dim", "2", Patterns::Integer(), "The problem dimension.");
656 *   prm.declare_entry("n_steps",
657 *   "10",
658 *   Patterns::Integer(0),
659 *   "Number of adaptive refinement steps.");
660 *   prm.declare_entry("smoother dampen",
661 *   "1.0",
662 *   Patterns::Double(0.0),
663 *   "Dampen factor for the smoother.");
664 *   prm.declare_entry("smoother steps",
665 *   "1",
666 *   Patterns::Integer(1),
667 *   "Number of smoother steps.");
668 *   prm.declare_entry("solver",
669 *   "MF",
670 *   Patterns::Selection("MF|MB|AMG"),
671 *   "Switch between matrix-free GMG, "
672 *   "matrix-based GMG, and AMG.");
673 *   prm.declare_entry("output",
674 *   "false",
675 *   Patterns::Bool(),
676 *   "Output graphical results.");
677 *  
678 *   if (prm_filename.size() == 0)
679 *   {
680 *   std::cout << "**** Error: No input file provided!\n"
681 *   << "**** Error: Call this program as './step-50 input.prm\n"
682 *   << '\n'
683 *   << "**** You may want to use one of the input files in this\n"
684 *   << "**** directory, or use the following default values\n"
685 *   << "**** to create an input file:\n";
686 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
687 *   prm.print_parameters(std::cout, ParameterHandler::Text);
688 *   return false;
689 *   }
690 *  
691 *   try
692 *   {
693 *   prm.parse_input(prm_filename);
694 *   }
695 *   catch (std::exception &e)
696 *   {
697 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
698 *   std::cerr << e.what() << std::endl;
699 *   return false;
700 *   }
701 *  
702 *   if (prm.get("solver") == "MF")
703 *   this->solver = gmg_mf;
704 *   else if (prm.get("solver") == "MB")
705 *   this->solver = gmg_mb;
706 *   else if (prm.get("solver") == "AMG")
707 *   this->solver = amg;
708 *   else
709 *   AssertThrow(false, ExcNotImplemented());
710 *  
711 *   this->dimension = prm.get_integer("dim");
712 *   this->n_steps = prm.get_integer("n_steps");
713 *   this->smoother_dampen = prm.get_double("smoother dampen");
714 *   this->smoother_steps = prm.get_integer("smoother steps");
715 *   this->output = prm.get_bool("output");
716 *  
717 *   return true;
718 *   }
719 *  
720 *  
721 *  
722 * @endcode
723 *
724 *
725 * <a name="LaplaceProblemclass"></a>
726 * <h3>LaplaceProblem class</h3>
727 *
728
729 *
730 * This is the main class of the program. It looks very similar to
731 * @ref step_16 "step-16", @ref step_37 "step-37", and @ref step_40 "step-40". For the MatrixFree setup, we use the
732 * MatrixFreeOperators::LaplaceOperator class which defines `local_apply()`,
733 * `compute_diagonal()`, and `set_coefficient()` functions internally. Note that
734 * the polynomial degree is a template parameter of this class. This is
735 * necessary for the matrix-free code.
736 *
737 * @code
738 *   template <int dim, int degree>
739 *   class LaplaceProblem
740 *   {
741 *   public:
742 *   LaplaceProblem(const Settings &settings);
743 *   void run();
744 *  
745 *   private:
746 * @endcode
747 *
748 * We will use the following types throughout the program. First the
749 * matrix-based types, after that the matrix-free classes. For the
750 * matrix-free implementation, we use @p float for the level operators.
751 *
752 * @code
753 *   using MatrixType = LA::MPI::SparseMatrix;
754 *   using VectorType = LA::MPI::Vector;
755 *   using PreconditionAMG = LA::MPI::PreconditionAMG;
756 *  
757 *   using MatrixFreeLevelMatrix = MatrixFreeOperators::LaplaceOperator<
758 *   dim,
759 *   degree,
760 *   degree + 1,
761 *   1,
763 *   using MatrixFreeActiveMatrix = MatrixFreeOperators::LaplaceOperator<
764 *   dim,
765 *   degree,
766 *   degree + 1,
767 *   1,
769 *  
770 *   using MatrixFreeLevelVector = LinearAlgebra::distributed::Vector<float>;
771 *   using MatrixFreeActiveVector = LinearAlgebra::distributed::Vector<double>;
772 *  
773 *   void setup_system();
774 *   void setup_multigrid();
775 *   void assemble_system();
776 *   void assemble_multigrid();
777 *   void assemble_rhs();
778 *   void solve();
779 *   void estimate();
780 *   void refine_grid();
781 *   void output_results(const unsigned int cycle);
782 *  
783 *   Settings settings;
784 *  
785 *   MPI_Comm mpi_communicator;
786 *   ConditionalOStream pcout;
787 *  
789 *   const MappingQ1<dim> mapping;
790 *   FE_Q<dim> fe;
791 *  
792 *   DoFHandler<dim> dof_handler;
793 *  
794 *   IndexSet locally_owned_dofs;
795 *   IndexSet locally_relevant_dofs;
796 *   AffineConstraints<double> constraints;
797 *  
798 *   MatrixType system_matrix;
799 *   MatrixFreeActiveMatrix mf_system_matrix;
800 *   VectorType solution;
801 *   VectorType right_hand_side;
802 *   Vector<double> estimated_error_square_per_cell;
803 *  
804 *   MGLevelObject<MatrixType> mg_matrix;
805 *   MGLevelObject<MatrixType> mg_interface_in;
806 *   MGConstrainedDoFs mg_constrained_dofs;
807 *  
808 *   MGLevelObject<MatrixFreeLevelMatrix> mf_mg_matrix;
809 *  
810 *   TimerOutput computing_timer;
811 *   };
812 *  
813 *  
814 * @endcode
815 *
816 * The only interesting part about the constructor is that we construct the
817 * multigrid hierarchy unless we use AMG. For that, we need to parse the
818 * run time parameters before this constructor completes.
819 *
820 * @code
821 *   template <int dim, int degree>
822 *   LaplaceProblem<dim, degree>::LaplaceProblem(const Settings &settings)
823 *   : settings(settings)
824 *   , mpi_communicator(MPI_COMM_WORLD)
825 *   , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
826 *   , triangulation(mpi_communicator,
828 *   (settings.solver == Settings::amg) ?
830 *   parallel::distributed::Triangulation<
832 *   , mapping()
833 *   , fe(degree)
834 *   , dof_handler(triangulation)
835 *   , computing_timer(pcout, TimerOutput::never, TimerOutput::wall_times)
836 *   {
837 *   GridGenerator::hyper_L(triangulation, -1., 1., /*colorize*/ false);
838 *   triangulation.refine_global(1);
839 *   }
840 *  
841 *  
842 *  
843 * @endcode
844 *
845 *
846 * <a name="LaplaceProblemsetup_system"></a>
847 * <h4>LaplaceProblem::setup_system()</h4>
848 *
849
850 *
851 * Unlike @ref step_16 "step-16" and @ref step_37 "step-37", we split the set up into two parts,
852 * setup_system() and setup_multigrid(). Here is the typical setup_system()
853 * function for the active mesh found in most tutorials. For matrix-free, the
854 * active mesh set up is similar to @ref step_37 "step-37"; for matrix-based (GMG and AMG
855 * solvers), the setup is similar to @ref step_40 "step-40".
856 *
857 * @code
858 *   template <int dim, int degree>
859 *   void LaplaceProblem<dim, degree>::setup_system()
860 *   {
861 *   TimerOutput::Scope timing(computing_timer, "Setup");
862 *  
863 *   dof_handler.distribute_dofs(fe);
864 *  
865 *   locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handler);
866 *   locally_owned_dofs = dof_handler.locally_owned_dofs();
867 *  
868 *   solution.reinit(locally_owned_dofs, mpi_communicator);
869 *   right_hand_side.reinit(locally_owned_dofs, mpi_communicator);
870 *   constraints.reinit(locally_relevant_dofs);
871 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
872 *  
874 *   mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
875 *   constraints.close();
876 *  
877 *   switch (settings.solver)
878 *   {
879 *   case Settings::gmg_mf:
880 *   {
881 *   typename MatrixFree<dim, double>::AdditionalData additional_data;
882 *   additional_data.tasks_parallel_scheme =
884 *   additional_data.mapping_update_flags =
886 *   std::shared_ptr<MatrixFree<dim, double>> mf_storage =
887 *   std::make_shared<MatrixFree<dim, double>>();
888 *   mf_storage->reinit(mapping,
889 *   dof_handler,
890 *   constraints,
891 *   QGauss<1>(degree + 1),
892 *   additional_data);
893 *  
894 *   mf_system_matrix.initialize(mf_storage);
895 *  
896 *   const Coefficient<dim> coefficient;
897 *   mf_system_matrix.set_coefficient(
898 *   coefficient.make_coefficient_table(*mf_storage));
899 *  
900 *   break;
901 *   }
902 *  
903 *   case Settings::gmg_mb:
904 *   case Settings::amg:
905 *   {
906 *   #ifdef USE_PETSC_LA
907 *   DynamicSparsityPattern dsp(locally_relevant_dofs);
908 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
909 *  
911 *   locally_owned_dofs,
912 *   mpi_communicator,
913 *   locally_relevant_dofs);
914 *  
915 *   system_matrix.reinit(locally_owned_dofs,
916 *   locally_owned_dofs,
917 *   dsp,
918 *   mpi_communicator);
919 *   #else
920 *   TrilinosWrappers::SparsityPattern dsp(locally_owned_dofs,
921 *   locally_owned_dofs,
922 *   locally_relevant_dofs,
923 *   mpi_communicator);
924 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
925 *   dsp.compress();
926 *   system_matrix.reinit(dsp);
927 *   #endif
928 *  
929 *   break;
930 *   }
931 *  
932 *   default:
933 *   Assert(false, ExcNotImplemented());
934 *   }
935 *   }
936 *  
937 * @endcode
938 *
939 *
940 * <a name="LaplaceProblemsetup_multigrid"></a>
941 * <h4>LaplaceProblem::setup_multigrid()</h4>
942 *
943
944 *
945 * This function does the multilevel setup for both matrix-free and
946 * matrix-based GMG. The matrix-free setup is similar to that of @ref step_37 "step-37", and
947 * the matrix-based is similar to @ref step_16 "step-16", except we must use appropriate
948 * distributed sparsity patterns.
949 *
950
951 *
952 * The function is not called for the AMG approach, but to err on the
953 * safe side, the main `switch` statement of this function
954 * nevertheless makes sure that the function only operates on known
955 * multigrid settings by throwing an assertion if the function were
956 * called for anything other than the two geometric multigrid methods.
957 *
958 * @code
959 *   template <int dim, int degree>
960 *   void LaplaceProblem<dim, degree>::setup_multigrid()
961 *   {
962 *   TimerOutput::Scope timing(computing_timer, "Setup multigrid");
963 *  
964 *   dof_handler.distribute_mg_dofs();
965 *  
966 *   mg_constrained_dofs.clear();
967 *   mg_constrained_dofs.initialize(dof_handler);
968 *  
969 *   const std::set<types::boundary_id> boundary_ids = {types::boundary_id(0)};
970 *   mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, boundary_ids);
971 *  
972 *   const unsigned int n_levels = triangulation.n_global_levels();
973 *  
974 *   switch (settings.solver)
975 *   {
976 *   case Settings::gmg_mf:
977 *   {
978 *   mf_mg_matrix.resize(0, n_levels - 1);
979 *  
980 *   for (unsigned int level = 0; level < n_levels; ++level)
981 *   {
982 *   const IndexSet relevant_dofs =
984 *   level);
985 *   AffineConstraints<double> level_constraints;
986 *   level_constraints.reinit(relevant_dofs);
987 *   level_constraints.add_lines(
988 *   mg_constrained_dofs.get_boundary_indices(level));
989 *   level_constraints.close();
990 *  
991 *   typename MatrixFree<dim, float>::AdditionalData additional_data;
992 *   additional_data.tasks_parallel_scheme =
994 *   additional_data.mapping_update_flags =
997 *   additional_data.mg_level = level;
998 *   std::shared_ptr<MatrixFree<dim, float>> mf_storage_level(
999 *   new MatrixFree<dim, float>());
1000 *   mf_storage_level->reinit(mapping,
1001 *   dof_handler,
1002 *   level_constraints,
1003 *   QGauss<1>(degree + 1),
1004 *   additional_data);
1005 *  
1006 *   mf_mg_matrix[level].initialize(mf_storage_level,
1007 *   mg_constrained_dofs,
1008 *   level);
1009 *  
1010 *   const Coefficient<dim> coefficient;
1011 *   mf_mg_matrix[level].set_coefficient(
1012 *   coefficient.make_coefficient_table(*mf_storage_level));
1013 *  
1014 *   mf_mg_matrix[level].compute_diagonal();
1015 *   }
1016 *  
1017 *   break;
1018 *   }
1019 *  
1020 *   case Settings::gmg_mb:
1021 *   {
1022 *   mg_matrix.resize(0, n_levels - 1);
1023 *   mg_matrix.clear_elements();
1024 *   mg_interface_in.resize(0, n_levels - 1);
1025 *   mg_interface_in.clear_elements();
1026 *  
1027 *   for (unsigned int level = 0; level < n_levels; ++level)
1028 *   {
1029 *   const IndexSet dof_set =
1031 *   level);
1032 *  
1033 *   {
1034 *   #ifdef USE_PETSC_LA
1035 *   DynamicSparsityPattern dsp(dof_set);
1036 *   MGTools::make_sparsity_pattern(dof_handler, dsp, level);
1037 *   dsp.compress();
1039 *   dsp,
1040 *   dof_handler.locally_owned_mg_dofs(level),
1041 *   mpi_communicator,
1042 *   dof_set);
1043 *  
1044 *   mg_matrix[level].reinit(
1045 *   dof_handler.locally_owned_mg_dofs(level),
1046 *   dof_handler.locally_owned_mg_dofs(level),
1047 *   dsp,
1048 *   mpi_communicator);
1049 *   #else
1051 *   dof_handler.locally_owned_mg_dofs(level),
1052 *   dof_handler.locally_owned_mg_dofs(level),
1053 *   dof_set,
1054 *   mpi_communicator);
1055 *   MGTools::make_sparsity_pattern(dof_handler, dsp, level);
1056 *  
1057 *   dsp.compress();
1058 *   mg_matrix[level].reinit(dsp);
1059 *   #endif
1060 *   }
1061 *  
1062 *   {
1063 *   #ifdef USE_PETSC_LA
1064 *   DynamicSparsityPattern dsp(dof_set);
1066 *   mg_constrained_dofs,
1067 *   dsp,
1068 *   level);
1069 *   dsp.compress();
1071 *   dsp,
1072 *   dof_handler.locally_owned_mg_dofs(level),
1073 *   mpi_communicator,
1074 *   dof_set);
1075 *  
1076 *   mg_interface_in[level].reinit(
1077 *   dof_handler.locally_owned_mg_dofs(level),
1078 *   dof_handler.locally_owned_mg_dofs(level),
1079 *   dsp,
1080 *   mpi_communicator);
1081 *   #else
1083 *   dof_handler.locally_owned_mg_dofs(level),
1084 *   dof_handler.locally_owned_mg_dofs(level),
1085 *   dof_set,
1086 *   mpi_communicator);
1087 *  
1089 *   mg_constrained_dofs,
1090 *   dsp,
1091 *   level);
1092 *   dsp.compress();
1093 *   mg_interface_in[level].reinit(dsp);
1094 *   #endif
1095 *   }
1096 *   }
1097 *   break;
1098 *   }
1099 *  
1100 *   default:
1101 *   Assert(false, ExcNotImplemented());
1102 *   }
1103 *   }
1104 *  
1105 *  
1106 * @endcode
1107 *
1108 *
1109 * <a name="LaplaceProblemassemble_system"></a>
1110 * <h4>LaplaceProblem::assemble_system()</h4>
1111 *
1112
1113 *
1114 * The assembly is split into three parts: `assemble_system()`,
1115 * `assemble_multigrid()`, and `assemble_rhs()`. The
1116 * `assemble_system()` function here assembles and stores the (global)
1117 * system matrix and the right-hand side for the matrix-based
1118 * methods. It is similar to the assembly in @ref step_40 "step-40".
1119 *
1120
1121 *
1122 * Note that the matrix-free method does not execute this function as it does
1123 * not need to assemble a matrix, and it will instead assemble the right-hand
1124 * side in assemble_rhs().
1125 *
1126 * @code
1127 *   template <int dim, int degree>
1128 *   void LaplaceProblem<dim, degree>::assemble_system()
1129 *   {
1130 *   TimerOutput::Scope timing(computing_timer, "Assemble");
1131 *  
1132 *   const QGauss<dim> quadrature_formula(degree + 1);
1133 *  
1134 *   FEValues<dim> fe_values(fe,
1135 *   quadrature_formula,
1138 *  
1139 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1140 *   const unsigned int n_q_points = quadrature_formula.size();
1141 *  
1142 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1143 *   Vector<double> cell_rhs(dofs_per_cell);
1144 *  
1145 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1146 *  
1147 *   const Coefficient<dim> coefficient;
1148 *   RightHandSide<dim> rhs;
1149 *   std::vector<double> rhs_values(n_q_points);
1150 *  
1151 *   for (const auto &cell : dof_handler.active_cell_iterators())
1152 *   if (cell->is_locally_owned())
1153 *   {
1154 *   cell_matrix = 0;
1155 *   cell_rhs = 0;
1156 *  
1157 *   fe_values.reinit(cell);
1158 *  
1159 *   const double coefficient_value =
1160 *   coefficient.average_value(fe_values.get_quadrature_points());
1161 *   rhs.value_list(fe_values.get_quadrature_points(), rhs_values);
1162 *  
1163 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1164 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1165 *   {
1166 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1167 *   cell_matrix(i, j) +=
1168 *   coefficient_value * // epsilon(x)
1169 *   fe_values.shape_grad(i, q_point) * // * grad phi_i(x)
1170 *   fe_values.shape_grad(j, q_point) * // * grad phi_j(x)
1171 *   fe_values.JxW(q_point); // * dx
1172 *  
1173 *   cell_rhs(i) +=
1174 *   fe_values.shape_value(i, q_point) * // grad phi_i(x)
1175 *   rhs_values[q_point] * // * f(x)
1176 *   fe_values.JxW(q_point); // * dx
1177 *   }
1178 *  
1179 *   cell->get_dof_indices(local_dof_indices);
1180 *   constraints.distribute_local_to_global(cell_matrix,
1181 *   cell_rhs,
1182 *   local_dof_indices,
1183 *   system_matrix,
1184 *   right_hand_side);
1185 *   }
1186 *  
1187 *   system_matrix.compress(VectorOperation::add);
1188 *   right_hand_side.compress(VectorOperation::add);
1189 *   }
1190 *  
1191 *  
1192 * @endcode
1193 *
1194 *
1195 * <a name="LaplaceProblemassemble_multigrid"></a>
1196 * <h4>LaplaceProblem::assemble_multigrid()</h4>
1197 *
1198
1199 *
1200 * The following function assembles and stores the multilevel matrices for the
1201 * matrix-based GMG method. This function is similar to the one found in
1202 * @ref step_16 "step-16", only here it works for distributed meshes. This difference amounts
1203 * to adding a condition that we only assemble on locally owned level cells and
1204 * a call to compress() for each matrix that is built.
1205 *
1206 * @code
1207 *   template <int dim, int degree>
1208 *   void LaplaceProblem<dim, degree>::assemble_multigrid()
1209 *   {
1210 *   TimerOutput::Scope timing(computing_timer, "Assemble multigrid");
1211 *  
1212 *   QGauss<dim> quadrature_formula(degree + 1);
1213 *  
1214 *   FEValues<dim> fe_values(fe,
1215 *   quadrature_formula,
1218 *  
1219 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1220 *   const unsigned int n_q_points = quadrature_formula.size();
1221 *  
1222 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1223 *  
1224 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1225 *  
1226 *   const Coefficient<dim> coefficient;
1227 *  
1228 *   std::vector<AffineConstraints<double>> boundary_constraints(
1229 *   triangulation.n_global_levels());
1230 *   for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
1231 *   {
1232 *   const IndexSet dof_set =
1234 *   boundary_constraints[level].reinit(dof_set);
1235 *   boundary_constraints[level].add_lines(
1236 *   mg_constrained_dofs.get_refinement_edge_indices(level));
1237 *   boundary_constraints[level].add_lines(
1238 *   mg_constrained_dofs.get_boundary_indices(level));
1239 *  
1240 *   boundary_constraints[level].close();
1241 *   }
1242 *  
1243 *   for (const auto &cell : dof_handler.cell_iterators())
1244 *   if (cell->level_subdomain_id() == triangulation.locally_owned_subdomain())
1245 *   {
1246 *   cell_matrix = 0;
1247 *   fe_values.reinit(cell);
1248 *  
1249 *   const double coefficient_value =
1250 *   coefficient.average_value(fe_values.get_quadrature_points());
1251 *  
1252 *   for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1253 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1254 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1255 *   cell_matrix(i, j) +=
1256 *   coefficient_value * fe_values.shape_grad(i, q_point) *
1257 *   fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point);
1258 *  
1259 *   cell->get_mg_dof_indices(local_dof_indices);
1260 *  
1261 *   boundary_constraints[cell->level()].distribute_local_to_global(
1262 *   cell_matrix, local_dof_indices, mg_matrix[cell->level()]);
1263 *  
1264 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1265 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1266 *   if (mg_constrained_dofs.is_interface_matrix_entry(
1267 *   cell->level(), local_dof_indices[i], local_dof_indices[j]))
1268 *   mg_interface_in[cell->level()].add(local_dof_indices[i],
1269 *   local_dof_indices[j],
1270 *   cell_matrix(i, j));
1271 *   }
1272 *  
1273 *   for (unsigned int i = 0; i < triangulation.n_global_levels(); ++i)
1274 *   {
1275 *   mg_matrix[i].compress(VectorOperation::add);
1276 *   mg_interface_in[i].compress(VectorOperation::add);
1277 *   }
1278 *   }
1279 *  
1280 *  
1281 *  
1282 * @endcode
1283 *
1284 *
1285 * <a name="LaplaceProblemassemble_rhs"></a>
1286 * <h4>LaplaceProblem::assemble_rhs()</h4>
1287 *
1288
1289 *
1290 * The final function in this triptych assembles the right-hand side
1291 * vector for the matrix-free method -- because in the matrix-free
1292 * framework, we don't have to assemble the matrix and can get away
1293 * with only assembling the right hand side. We could do this by extracting the
1294 * code from the `assemble_system()` function above that deals with the right
1295 * hand side, but we decide instead to go all in on the matrix-free approach and
1296 * do the assembly using that way as well.
1297 *
1298
1299 *
1300 * The result is a function that is similar
1301 * to the one found in the "Use FEEvaluation::read_dof_values_plain()
1302 * to avoid resolving constraints" subsection in the "Possibilities
1303 * for extensions" section of @ref step_37 "step-37".
1304 *
1305
1306 *
1307 * The reason for this function is that the MatrixFree operators do not take
1308 * into account non-homogeneous Dirichlet constraints, instead treating all
1309 * Dirichlet constraints as homogeneous. To account for this, the right-hand
1310 * side here is assembled as the residual @f$r_0 = f-Au_0@f$, where @f$u_0@f$ is a
1311 * zero vector except in the Dirichlet values. Then when solving, we have that
1312 * the solution is @f$u = u_0 + A^{-1}r_0@f$. This can be seen as a Newton
1313 * iteration on a linear system with initial guess @f$u_0@f$. The CG solve in the
1314 * `solve()` function below computes @f$A^{-1}r_0@f$ and the call to
1315 * `constraints.distribute()` (which directly follows) adds the @f$u_0@f$.
1316 *
1317
1318 *
1319 * Obviously, since we are considering a problem with zero Dirichlet boundary,
1320 * we could have taken a similar approach to @ref step_37 "step-37" `assemble_rhs()`, but this
1321 * additional work allows us to change the problem declaration if we so
1322 * choose.
1323 *
1324
1325 *
1326 * This function has two parts in the integration loop: applying the negative
1327 * of matrix @f$A@f$ to @f$u_0@f$ by submitting the negative of the gradient, and adding
1328 * the right-hand side contribution by submitting the value @f$f@f$. We must be sure
1329 * to use `read_dof_values_plain()` for evaluating @f$u_0@f$ as `read_dof_values()`
1330 * would set all Dirichlet values to zero.
1331 *
1332
1333 *
1334 * Finally, the system_rhs vector is of type LA::MPI::Vector, but the
1335 * MatrixFree class only work for
1336 * LinearAlgebra::distributed::Vector. Therefore we must
1337 * compute the right-hand side using MatrixFree functionality and then
1338 * use the functions in the `ChangeVectorType` namespace to copy it to
1339 * the correct type.
1340 *
1341 * @code
1342 *   template <int dim, int degree>
1343 *   void LaplaceProblem<dim, degree>::assemble_rhs()
1344 *   {
1345 *   TimerOutput::Scope timing(computing_timer, "Assemble right-hand side");
1346 *  
1347 *   MatrixFreeActiveVector solution_copy;
1348 *   MatrixFreeActiveVector right_hand_side_copy;
1349 *   mf_system_matrix.initialize_dof_vector(solution_copy);
1350 *   mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
1351 *  
1352 *   solution_copy = 0.;
1353 *   constraints.distribute(solution_copy);
1354 *   solution_copy.update_ghost_values();
1355 *   right_hand_side_copy = 0;
1356 *   const Table<2, VectorizedArray<double>> &coefficient =
1357 *   *(mf_system_matrix.get_coefficient());
1358 *  
1359 *   RightHandSide<dim> right_hand_side_function;
1360 *  
1361 *   FEEvaluation<dim, degree, degree + 1, 1, double> phi(
1362 *   *mf_system_matrix.get_matrix_free());
1363 *  
1364 *   for (unsigned int cell = 0;
1365 *   cell < mf_system_matrix.get_matrix_free()->n_cell_batches();
1366 *   ++cell)
1367 *   {
1368 *   phi.reinit(cell);
1369 *   phi.read_dof_values_plain(solution_copy);
1370 *   phi.evaluate(EvaluationFlags::gradients);
1371 *  
1372 *   for (unsigned int q = 0; q < phi.n_q_points; ++q)
1373 *   {
1374 *   phi.submit_gradient(-1.0 *
1375 *   (coefficient(cell, 0) * phi.get_gradient(q)),
1376 *   q);
1377 *   phi.submit_value(
1378 *   right_hand_side_function.value(phi.quadrature_point(q)), q);
1379 *   }
1380 *  
1381 *   phi.integrate_scatter(EvaluationFlags::values |
1382 *   EvaluationFlags::gradients,
1383 *   right_hand_side_copy);
1384 *   }
1385 *  
1386 *   right_hand_side_copy.compress(VectorOperation::add);
1387 *  
1388 *   ChangeVectorTypes::copy(right_hand_side, right_hand_side_copy);
1389 *   }
1390 *  
1391 *  
1392 *  
1393 * @endcode
1394 *
1395 *
1396 * <a name="LaplaceProblemsolve"></a>
1397 * <h4>LaplaceProblem::solve()</h4>
1398 *
1399
1400 *
1401 * Here we set up the multigrid preconditioner, test the timing of a single
1402 * V-cycle, and solve the linear system. Unsurprisingly, this is one of the
1403 * places where the three methods differ the most.
1404 *
1405 * @code
1406 *   template <int dim, int degree>
1407 *   void LaplaceProblem<dim, degree>::solve()
1408 *   {
1409 *   TimerOutput::Scope timing(computing_timer, "Solve");
1410 *  
1411 *   SolverControl solver_control(1000, 1.e-10 * right_hand_side.l2_norm());
1412 *   solver_control.enable_history_data();
1413 *  
1414 *   solution = 0.;
1415 *  
1416 * @endcode
1417 *
1418 * The solver for the matrix-free GMG method is similar to @ref step_37 "step-37", apart
1419 * from adding some interface matrices in complete analogy to @ref step_16 "step-16".
1420 *
1421 * @code
1422 *   switch (settings.solver)
1423 *   {
1424 *   case Settings::gmg_mf:
1425 *   {
1426 *   computing_timer.enter_subsection("Solve: Preconditioner setup");
1427 *  
1428 *   MGTransferMatrixFree<dim, float> mg_transfer(mg_constrained_dofs);
1429 *   mg_transfer.build(dof_handler);
1430 *  
1431 *   SolverControl coarse_solver_control(1000, 1e-12, false, false);
1432 *   SolverCG<MatrixFreeLevelVector> coarse_solver(coarse_solver_control);
1433 *   PreconditionIdentity identity;
1434 *   MGCoarseGridIterativeSolver<MatrixFreeLevelVector,
1435 *   SolverCG<MatrixFreeLevelVector>,
1436 *   MatrixFreeLevelMatrix,
1437 *   PreconditionIdentity>
1438 *   coarse_grid_solver(coarse_solver, mf_mg_matrix[0], identity);
1439 *  
1440 *   using Smoother = PreconditionJacobi<MatrixFreeLevelMatrix>;
1441 *   MGSmootherPrecondition<MatrixFreeLevelMatrix,
1442 *   Smoother,
1443 *   MatrixFreeLevelVector>
1444 *   smoother;
1445 *   smoother.initialize(mf_mg_matrix,
1446 *   typename Smoother::AdditionalData(
1447 *   settings.smoother_dampen));
1448 *   smoother.set_steps(settings.smoother_steps);
1449 *  
1450 *   mg::Matrix<MatrixFreeLevelVector> mg_m(mf_mg_matrix);
1451 *  
1452 *   MGLevelObject<
1453 *   MatrixFreeOperators::MGInterfaceOperator<MatrixFreeLevelMatrix>>
1454 *   mg_interface_matrices;
1455 *   mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1);
1456 *   for (unsigned int level = 0; level < triangulation.n_global_levels();
1457 *   ++level)
1458 *   mg_interface_matrices[level].initialize(mf_mg_matrix[level]);
1459 *   mg::Matrix<MatrixFreeLevelVector> mg_interface(mg_interface_matrices);
1460 *  
1461 *   Multigrid<MatrixFreeLevelVector> mg(
1462 *   mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
1463 *   mg.set_edge_matrices(mg_interface, mg_interface);
1464 *  
1465 *   PreconditionMG<dim,
1466 *   MatrixFreeLevelVector,
1467 *   MGTransferMatrixFree<dim, float>>
1468 *   preconditioner(dof_handler, mg, mg_transfer);
1469 *  
1470 * @endcode
1471 *
1472 * Copy the solution vector and right-hand side from LA::MPI::Vector
1473 * to LinearAlgebra::distributed::Vector so that we can solve.
1474 *
1475 * @code
1476 *   MatrixFreeActiveVector solution_copy;
1477 *   MatrixFreeActiveVector right_hand_side_copy;
1478 *   mf_system_matrix.initialize_dof_vector(solution_copy);
1479 *   mf_system_matrix.initialize_dof_vector(right_hand_side_copy);
1480 *  
1481 *   ChangeVectorTypes::copy(solution_copy, solution);
1482 *   ChangeVectorTypes::copy(right_hand_side_copy, right_hand_side);
1483 *   computing_timer.leave_subsection("Solve: Preconditioner setup");
1484 *  
1485 * @endcode
1486 *
1487 * Timing for 1 V-cycle.
1488 *
1489 * @code
1490 *   {
1491 *   TimerOutput::Scope timing(computing_timer,
1492 *   "Solve: 1 multigrid V-cycle");
1493 *   preconditioner.vmult(solution_copy, right_hand_side_copy);
1494 *   }
1495 *   solution_copy = 0.;
1496 *  
1497 * @endcode
1498 *
1499 * Solve the linear system, update the ghost values of the solution,
1500 * copy back to LA::MPI::Vector and distribute constraints.
1501 *
1502 * @code
1503 *   {
1504 *   SolverCG<MatrixFreeActiveVector> solver(solver_control);
1505 *  
1506 *   TimerOutput::Scope timing(computing_timer, "Solve: CG");
1507 *   solver.solve(mf_system_matrix,
1508 *   solution_copy,
1509 *   right_hand_side_copy,
1510 *   preconditioner);
1511 *   }
1512 *  
1513 *   solution_copy.update_ghost_values();
1514 *   ChangeVectorTypes::copy(solution, solution_copy);
1515 *   constraints.distribute(solution);
1516 *  
1517 *   break;
1518 *   }
1519 *  
1520 * @endcode
1521 *
1522 * Solver for the matrix-based GMG method, similar to @ref step_16 "step-16", only
1523 * using a Jacobi smoother instead of a SOR smoother (which is not
1524 * implemented in parallel).
1525 *
1526 * @code
1527 *   case Settings::gmg_mb:
1528 *   {
1529 *   computing_timer.enter_subsection("Solve: Preconditioner setup");
1530 *  
1531 *   MGTransferPrebuilt<VectorType> mg_transfer(mg_constrained_dofs);
1532 *   mg_transfer.build(dof_handler);
1533 *  
1534 *   SolverControl coarse_solver_control(1000, 1e-12, false, false);
1535 *   SolverCG<VectorType> coarse_solver(coarse_solver_control);
1536 *   PreconditionIdentity identity;
1537 *   MGCoarseGridIterativeSolver<VectorType,
1538 *   SolverCG<VectorType>,
1539 *   MatrixType,
1540 *   PreconditionIdentity>
1541 *   coarse_grid_solver(coarse_solver, mg_matrix[0], identity);
1542 *  
1543 *   using Smoother = LA::MPI::PreconditionJacobi;
1544 *   MGSmootherPrecondition<MatrixType, Smoother, VectorType> smoother;
1545 *  
1546 *   #ifdef USE_PETSC_LA
1547 *   smoother.initialize(mg_matrix);
1548 *   Assert(
1549 *   settings.smoother_dampen == 1.0,
1550 *   ExcNotImplemented(
1551 *   "PETSc's PreconditionJacobi has no support for a damping parameter."));
1552 *   #else
1553 *   smoother.initialize(mg_matrix, settings.smoother_dampen);
1554 *   #endif
1555 *  
1556 *   smoother.set_steps(settings.smoother_steps);
1557 *  
1558 *   mg::Matrix<VectorType> mg_m(mg_matrix);
1559 *   mg::Matrix<VectorType> mg_in(mg_interface_in);
1560 *   mg::Matrix<VectorType> mg_out(mg_interface_in);
1561 *  
1562 *   Multigrid<VectorType> mg(
1563 *   mg_m, coarse_grid_solver, mg_transfer, smoother, smoother);
1564 *   mg.set_edge_matrices(mg_out, mg_in);
1565 *  
1566 *  
1567 *   PreconditionMG<dim, VectorType, MGTransferPrebuilt<VectorType>>
1568 *   preconditioner(dof_handler, mg, mg_transfer);
1569 *  
1570 *   computing_timer.leave_subsection("Solve: Preconditioner setup");
1571 *  
1572 * @endcode
1573 *
1574 * Timing for 1 V-cycle.
1575 *
1576 * @code
1577 *   {
1578 *   TimerOutput::Scope timing(computing_timer,
1579 *   "Solve: 1 multigrid V-cycle");
1580 *   preconditioner.vmult(solution, right_hand_side);
1581 *   }
1582 *   solution = 0.;
1583 *  
1584 * @endcode
1585 *
1586 * Solve the linear system and distribute constraints.
1587 *
1588 * @code
1589 *   {
1590 *   SolverCG<VectorType> solver(solver_control);
1591 *  
1592 *   TimerOutput::Scope timing(computing_timer, "Solve: CG");
1593 *   solver.solve(system_matrix,
1594 *   solution,
1595 *   right_hand_side,
1596 *   preconditioner);
1597 *   }
1598 *  
1599 *   constraints.distribute(solution);
1600 *  
1601 *   break;
1602 *   }
1603 *  
1604 * @endcode
1605 *
1606 * Solver for the AMG method, similar to @ref step_40 "step-40".
1607 *
1608 * @code
1609 *   case Settings::amg:
1610 *   {
1611 *   computing_timer.enter_subsection("Solve: Preconditioner setup");
1612 *  
1613 *   PreconditionAMG preconditioner;
1614 *   PreconditionAMG::AdditionalData Amg_data;
1615 *  
1616 *   #ifdef USE_PETSC_LA
1617 *   Amg_data.symmetric_operator = true;
1618 *   #else
1619 *   Amg_data.elliptic = true;
1620 *   Amg_data.smoother_type = "Jacobi";
1621 *   Amg_data.higher_order_elements = true;
1622 *   Amg_data.smoother_sweeps = settings.smoother_steps;
1623 *   Amg_data.aggregation_threshold = 0.02;
1624 *   #endif
1625 *  
1626 *   Amg_data.output_details = false;
1627 *  
1628 *   preconditioner.initialize(system_matrix, Amg_data);
1629 *   computing_timer.leave_subsection("Solve: Preconditioner setup");
1630 *  
1631 * @endcode
1632 *
1633 * Timing for 1 V-cycle.
1634 *
1635 * @code
1636 *   {
1637 *   TimerOutput::Scope timing(computing_timer,
1638 *   "Solve: 1 multigrid V-cycle");
1639 *   preconditioner.vmult(solution, right_hand_side);
1640 *   }
1641 *   solution = 0.;
1642 *  
1643 * @endcode
1644 *
1645 * Solve the linear system and distribute constraints.
1646 *
1647 * @code
1648 *   {
1649 *   SolverCG<VectorType> solver(solver_control);
1650 *  
1651 *   TimerOutput::Scope timing(computing_timer, "Solve: CG");
1652 *   solver.solve(system_matrix,
1653 *   solution,
1654 *   right_hand_side,
1655 *   preconditioner);
1656 *   }
1657 *   constraints.distribute(solution);
1658 *  
1659 *   break;
1660 *   }
1661 *  
1662 *   default:
1663 *   Assert(false, ExcInternalError());
1664 *   }
1665 *  
1666 *   pcout << " Number of CG iterations: " << solver_control.last_step()
1667 *   << std::endl;
1668 *   }
1669 *  
1670 *  
1671 * @endcode
1672 *
1673 *
1674 * <a name="Theerrorestimator"></a>
1675 * <h3>The error estimator</h3>
1676 *
1677
1678 *
1679 * We use the FEInterfaceValues class to assemble an error estimator to decide
1680 * which cells to refine. See the exact definition of the cell and face
1681 * integrals in the introduction. To use the method, we define Scratch and
1682 * Copy objects for the MeshWorker::mesh_loop() with much of the following code
1683 * being in essence as was set up in @ref step_12 "step-12" already (or at least similar in
1684 * spirit).
1685 *
1686 * @code
1687 *   template <int dim>
1688 *   struct ScratchData
1689 *   {
1690 *   ScratchData(const Mapping<dim> & mapping,
1691 *   const FiniteElement<dim> &fe,
1692 *   const unsigned int quadrature_degree,
1693 *   const UpdateFlags update_flags,
1694 *   const UpdateFlags interface_update_flags)
1695 *   : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
1696 *   , fe_interface_values(mapping,
1697 *   fe,
1698 *   QGauss<dim - 1>(quadrature_degree),
1699 *   interface_update_flags)
1700 *   {}
1701 *  
1702 *  
1703 *   ScratchData(const ScratchData<dim> &scratch_data)
1704 *   : fe_values(scratch_data.fe_values.get_mapping(),
1705 *   scratch_data.fe_values.get_fe(),
1706 *   scratch_data.fe_values.get_quadrature(),
1707 *   scratch_data.fe_values.get_update_flags())
1708 *   , fe_interface_values(scratch_data.fe_values.get_mapping(),
1709 *   scratch_data.fe_values.get_fe(),
1710 *   scratch_data.fe_interface_values.get_quadrature(),
1711 *   scratch_data.fe_interface_values.get_update_flags())
1712 *   {}
1713 *  
1714 *   FEValues<dim> fe_values;
1715 *   FEInterfaceValues<dim> fe_interface_values;
1716 *   };
1717 *  
1718 *  
1719 *  
1720 *   struct CopyData
1721 *   {
1722 *   CopyData()
1723 *   : cell_index(numbers::invalid_unsigned_int)
1724 *   , value(0.)
1725 *   {}
1726 *  
1727 *   struct FaceData
1728 *   {
1729 *   unsigned int cell_indices[2];
1730 *   double values[2];
1731 *   };
1732 *  
1733 *   unsigned int cell_index;
1734 *   double value;
1735 *   std::vector<FaceData> face_data;
1736 *   };
1737 *  
1738 *  
1739 *   template <int dim, int degree>
1740 *   void LaplaceProblem<dim, degree>::estimate()
1741 *   {
1742 *   TimerOutput::Scope timing(computing_timer, "Estimate");
1743 *  
1744 *   VectorType temp_solution;
1745 *   temp_solution.reinit(locally_owned_dofs,
1746 *   locally_relevant_dofs,
1747 *   mpi_communicator);
1748 *   temp_solution = solution;
1749 *  
1750 *   const Coefficient<dim> coefficient;
1751 *  
1752 *   estimated_error_square_per_cell.reinit(triangulation.n_active_cells());
1753 *  
1754 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
1755 *  
1756 * @endcode
1757 *
1758 * Assembler for cell residual @f$h^2 \| f + \epsilon \triangle u \|_K^2@f$
1759 *
1760 * @code
1761 *   auto cell_worker = [&](const Iterator & cell,
1762 *   ScratchData<dim> &scratch_data,
1763 *   CopyData & copy_data) {
1764 *   FEValues<dim> &fe_values = scratch_data.fe_values;
1765 *   fe_values.reinit(cell);
1766 *  
1767 *   RightHandSide<dim> rhs;
1768 *   const double rhs_value = rhs.value(cell->center());
1769 *  
1770 *   const double nu = coefficient.value(cell->center());
1771 *  
1772 *   std::vector<Tensor<2, dim>> hessians(fe_values.n_quadrature_points);
1773 *   fe_values.get_function_hessians(temp_solution, hessians);
1774 *  
1775 *   copy_data.cell_index = cell->active_cell_index();
1776 *  
1777 *   double residual_norm_square = 0.;
1778 *   for (unsigned k = 0; k < fe_values.n_quadrature_points; ++k)
1779 *   {
1780 *   const double residual = (rhs_value + nu * trace(hessians[k]));
1781 *   residual_norm_square += residual * residual * fe_values.JxW(k);
1782 *   }
1783 *  
1784 *   copy_data.value =
1785 *   cell->diameter() * cell->diameter() * residual_norm_square;
1786 *   };
1787 *  
1788 * @endcode
1789 *
1790 * Assembler for face term @f$\sum_F h_F \| \jump{\epsilon \nabla u \cdot n}
1791 * \|_F^2@f$
1792 *
1793 * @code
1794 *   auto face_worker = [&](const Iterator & cell,
1795 *   const unsigned int &f,
1796 *   const unsigned int &sf,
1797 *   const Iterator & ncell,
1798 *   const unsigned int &nf,
1799 *   const unsigned int &nsf,
1800 *   ScratchData<dim> & scratch_data,
1801 *   CopyData & copy_data) {
1802 *   FEInterfaceValues<dim> &fe_interface_values =
1803 *   scratch_data.fe_interface_values;
1804 *   fe_interface_values.reinit(cell, f, sf, ncell, nf, nsf);
1805 *  
1806 *   copy_data.face_data.emplace_back();
1807 *   CopyData::FaceData &copy_data_face = copy_data.face_data.back();
1808 *  
1809 *   copy_data_face.cell_indices[0] = cell->active_cell_index();
1810 *   copy_data_face.cell_indices[1] = ncell->active_cell_index();
1811 *  
1812 *   const double coeff1 = coefficient.value(cell->center());
1813 *   const double coeff2 = coefficient.value(ncell->center());
1814 *  
1815 *   std::vector<Tensor<1, dim>> grad_u[2];
1816 *  
1817 *   for (unsigned int i = 0; i < 2; ++i)
1818 *   {
1819 *   grad_u[i].resize(fe_interface_values.n_quadrature_points);
1820 *   fe_interface_values.get_fe_face_values(i).get_function_gradients(
1821 *   temp_solution, grad_u[i]);
1822 *   }
1823 *  
1824 *   double jump_norm_square = 0.;
1825 *  
1826 *   for (unsigned int qpoint = 0;
1827 *   qpoint < fe_interface_values.n_quadrature_points;
1828 *   ++qpoint)
1829 *   {
1830 *   const double jump =
1831 *   coeff1 * grad_u[0][qpoint] * fe_interface_values.normal(qpoint) -
1832 *   coeff2 * grad_u[1][qpoint] * fe_interface_values.normal(qpoint);
1833 *  
1834 *   jump_norm_square += jump * jump * fe_interface_values.JxW(qpoint);
1835 *   }
1836 *  
1837 *   const double h = cell->face(f)->measure();
1838 *   copy_data_face.values[0] = 0.5 * h * jump_norm_square;
1839 *   copy_data_face.values[1] = copy_data_face.values[0];
1840 *   };
1841 *  
1842 *   auto copier = [&](const CopyData &copy_data) {
1843 *   if (copy_data.cell_index != numbers::invalid_unsigned_int)
1844 *   estimated_error_square_per_cell[copy_data.cell_index] += copy_data.value;
1845 *  
1846 *   for (auto &cdf : copy_data.face_data)
1847 *   for (unsigned int j = 0; j < 2; ++j)
1848 *   estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
1849 *   };
1850 *  
1851 *   const unsigned int n_gauss_points = degree + 1;
1852 *   ScratchData<dim> scratch_data(mapping,
1853 *   fe,
1854 *   n_gauss_points,
1855 *   update_hessians | update_quadrature_points |
1856 *   update_JxW_values,
1857 *   update_values | update_gradients |
1858 *   update_JxW_values | update_normal_vectors);
1859 *   CopyData copy_data;
1860 *  
1861 * @endcode
1862 *
1863 * We need to assemble each interior face once but we need to make sure that
1864 * both processes assemble the face term between a locally owned and a ghost
1865 * cell. This is achieved by setting the
1866 * MeshWorker::assemble_ghost_faces_both flag. We need to do this, because
1867 * we do not communicate the error estimator contributions here.
1868 *
1869 * @code
1870 *   MeshWorker::mesh_loop(dof_handler.begin_active(),
1871 *   dof_handler.end(),
1872 *   cell_worker,
1873 *   copier,
1874 *   scratch_data,
1875 *   copy_data,
1876 *   MeshWorker::assemble_own_cells |
1877 *   MeshWorker::assemble_ghost_faces_both |
1878 *   MeshWorker::assemble_own_interior_faces_once,
1879 *   /*boundary_worker=*/nullptr,
1880 *   face_worker);
1881 *  
1882 *   const double global_error_estimate =
1883 *   std::sqrt(Utilities::MPI::sum(estimated_error_square_per_cell.l1_norm(),
1884 *   mpi_communicator));
1885 *   pcout << " Global error estimate: " << global_error_estimate
1886 *   << std::endl;
1887 *   }
1888 *  
1889 *  
1890 * @endcode
1891 *
1892 *
1893 * <a name="LaplaceProblemrefine_grid"></a>
1894 * <h4>LaplaceProblem::refine_grid()</h4>
1895 *
1896
1897 *
1898 * We use the cell-wise estimator stored in the vector @p estimate_vector and
1899 * refine a fixed number of cells (chosen here to roughly double the number of
1900 * DoFs in each step).
1901 *
1902 * @code
1903 *   template <int dim, int degree>
1904 *   void LaplaceProblem<dim, degree>::refine_grid()
1905 *   {
1906 *   TimerOutput::Scope timing(computing_timer, "Refine grid");
1907 *  
1908 *   const double refinement_fraction = 1. / (std::pow(2.0, dim) - 1.);
1909 *   parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
1910 *   triangulation, estimated_error_square_per_cell, refinement_fraction, 0.0);
1911 *  
1912 *   triangulation.execute_coarsening_and_refinement();
1913 *   }
1914 *  
1915 *  
1916 * @endcode
1917 *
1918 *
1919 * <a name="LaplaceProblemoutput_results"></a>
1920 * <h4>LaplaceProblem::output_results()</h4>
1921 *
1922
1923 *
1924 * The output_results() function is similar to the ones found in many of the
1925 * tutorials (see @ref step_40 "step-40" for example).
1926 *
1927 * @code
1928 *   template <int dim, int degree>
1929 *   void LaplaceProblem<dim, degree>::output_results(const unsigned int cycle)
1930 *   {
1931 *   TimerOutput::Scope timing(computing_timer, "Output results");
1932 *  
1933 *   VectorType temp_solution;
1934 *   temp_solution.reinit(locally_owned_dofs,
1935 *   locally_relevant_dofs,
1936 *   mpi_communicator);
1937 *   temp_solution = solution;
1938 *  
1939 *   DataOut<dim> data_out;
1940 *   data_out.attach_dof_handler(dof_handler);
1941 *   data_out.add_data_vector(temp_solution, "solution");
1942 *  
1943 *   Vector<float> subdomain(triangulation.n_active_cells());
1944 *   for (unsigned int i = 0; i < subdomain.size(); ++i)
1945 *   subdomain(i) = triangulation.locally_owned_subdomain();
1946 *   data_out.add_data_vector(subdomain, "subdomain");
1947 *  
1948 *   Vector<float> level(triangulation.n_active_cells());
1949 *   for (const auto &cell : triangulation.active_cell_iterators())
1950 *   level(cell->active_cell_index()) = cell->level();
1951 *   data_out.add_data_vector(level, "level");
1952 *  
1953 *   if (estimated_error_square_per_cell.size() > 0)
1954 *   data_out.add_data_vector(estimated_error_square_per_cell,
1955 *   "estimated_error_square_per_cell");
1956 *  
1957 *   data_out.build_patches();
1958 *  
1959 *   const std::string pvtu_filename = data_out.write_vtu_with_pvtu_record(
1960 *   "", "solution", cycle, mpi_communicator, 2 /*n_digits*/, 1 /*n_groups*/);
1961 *  
1962 *   pcout << " Wrote " << pvtu_filename << std::endl;
1963 *   }
1964 *  
1965 *  
1966 * @endcode
1967 *
1968 *
1969 * <a name="LaplaceProblemrun"></a>
1970 * <h4>LaplaceProblem::run()</h4>
1971 *
1972
1973 *
1974 * As in most tutorials, this function calls the various functions defined
1975 * above to set up, assemble, solve, and output the results.
1976 *
1977 * @code
1978 *   template <int dim, int degree>
1979 *   void LaplaceProblem<dim, degree>::run()
1980 *   {
1981 *   for (unsigned int cycle = 0; cycle < settings.n_steps; ++cycle)
1982 *   {
1983 *   pcout << "Cycle " << cycle << ':' << std::endl;
1984 *   if (cycle > 0)
1985 *   refine_grid();
1986 *  
1987 *   pcout << " Number of active cells: "
1988 *   << triangulation.n_global_active_cells();
1989 *  
1990 * @endcode
1991 *
1992 * We only output level cell data for the GMG methods (same with DoF
1993 * data below). Note that the partition efficiency is irrelevant for AMG
1994 * since the level hierarchy is not distributed or used during the
1995 * computation.
1996 *
1997 * @code
1998 *   if (settings.solver == Settings::gmg_mf ||
1999 *   settings.solver == Settings::gmg_mb)
2000 *   pcout << " (" << triangulation.n_global_levels() << " global levels)"
2001 *   << std::endl
2002 *   << " Partition efficiency: "
2003 *   << 1.0 / MGTools::workload_imbalance(triangulation);
2004 *   pcout << std::endl;
2005 *  
2006 *   setup_system();
2007 *  
2008 * @endcode
2009 *
2010 * Only set up the multilevel hierarchy for GMG.
2011 *
2012 * @code
2013 *   if (settings.solver == Settings::gmg_mf ||
2014 *   settings.solver == Settings::gmg_mb)
2015 *   setup_multigrid();
2016 *  
2017 *   pcout << " Number of degrees of freedom: " << dof_handler.n_dofs();
2018 *   if (settings.solver == Settings::gmg_mf ||
2019 *   settings.solver == Settings::gmg_mb)
2020 *   {
2021 *   pcout << " (by level: ";
2022 *   for (unsigned int level = 0; level < triangulation.n_global_levels();
2023 *   ++level)
2024 *   pcout << dof_handler.n_dofs(level)
2025 *   << (level == triangulation.n_global_levels() - 1 ? ")" :
2026 *   ", ");
2027 *   }
2028 *   pcout << std::endl;
2029 *  
2030 * @endcode
2031 *
2032 * For the matrix-free method, we only assemble the right-hand side.
2033 * For both matrix-based methods, we assemble both active matrix and
2034 * right-hand side, and only assemble the multigrid matrices for
2035 * matrix-based GMG.
2036 *
2037 * @code
2038 *   if (settings.solver == Settings::gmg_mf)
2039 *   assemble_rhs();
2040 *   else /*gmg_mb or amg*/
2041 *   {
2042 *   assemble_system();
2043 *   if (settings.solver == Settings::gmg_mb)
2044 *   assemble_multigrid();
2045 *   }
2046 *  
2047 *   solve();
2048 *   estimate();
2049 *  
2050 *   if (settings.output)
2051 *   output_results(cycle);
2052 *  
2053 *   computing_timer.print_summary();
2054 *   computing_timer.reset();
2055 *   }
2056 *   }
2057 *  
2058 *  
2059 * @endcode
2060 *
2061 *
2062 * <a name="Themainfunction"></a>
2063 * <h3>The main() function</h3>
2064 *
2065
2066 *
2067 * This is a similar main function to @ref step_40 "step-40", with the exception that
2068 * we require the user to pass a .prm file as a sole command line
2069 * argument (see @ref step_29 "step-29" and the documentation of the ParameterHandler
2070 * class for a complete discussion of parameter files).
2071 *
2072 * @code
2073 *   int main(int argc, char *argv[])
2074 *   {
2075 *   using namespace dealii;
2076 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
2077 *  
2078 *   Settings settings;
2079 *   if (!settings.try_parse((argc > 1) ? (argv[1]) : ""))
2080 *   return 0;
2081 *  
2082 *   try
2083 *   {
2084 *   constexpr unsigned int fe_degree = 2;
2085 *  
2086 *   switch (settings.dimension)
2087 *   {
2088 *   case 2:
2089 *   {
2090 *   LaplaceProblem<2, fe_degree> test(settings);
2091 *   test.run();
2092 *  
2093 *   break;
2094 *   }
2095 *  
2096 *   case 3:
2097 *   {
2098 *   LaplaceProblem<3, fe_degree> test(settings);
2099 *   test.run();
2100 *  
2101 *   break;
2102 *   }
2103 *  
2104 *   default:
2105 *   Assert(false, ExcMessage("This program only works in 2d and 3d."));
2106 *   }
2107 *   }
2108 *   catch (std::exception &exc)
2109 *   {
2110 *   std::cerr << std::endl
2111 *   << std::endl
2112 *   << "----------------------------------------------------"
2113 *   << std::endl;
2114 *   std::cerr << "Exception on processing: " << std::endl
2115 *   << exc.what() << std::endl
2116 *   << "Aborting!" << std::endl
2117 *   << "----------------------------------------------------"
2118 *   << std::endl;
2119 *   MPI_Abort(MPI_COMM_WORLD, 1);
2120 *   return 1;
2121 *   }
2122 *   catch (...)
2123 *   {
2124 *   std::cerr << std::endl
2125 *   << std::endl
2126 *   << "----------------------------------------------------"
2127 *   << std::endl;
2128 *   std::cerr << "Unknown exception!" << std::endl
2129 *   << "Aborting!" << std::endl
2130 *   << "----------------------------------------------------"
2131 *   << std::endl;
2132 *   MPI_Abort(MPI_COMM_WORLD, 2);
2133 *   return 1;
2134 *   }
2135 *  
2136 *   return 0;
2137 *   }
2138 * @endcode
2139<a name="Results"></a><h1>Results</h1>
2140
2141
2142When you run the program using the following command
2143@code
2144mpirun -np 16 ./step-50 gmg_mf_2d.prm
2145@endcode
2146the screen output should look like the following:
2147@code
2148Cycle 0:
2149 Number of active cells: 12 (2 global levels)
2150 Partition efficiency: 0.1875
2151 Number of degrees of freedom: 65 (by level: 21, 65)
2152 Number of CG iterations: 10
2153 Global error estimate: 0.355373
2154 Wrote solution_00.pvtu
2155
2156
2157+---------------------------------------------+------------+------------+
2158| Total wallclock time elapsed since start | 0.0163s | |
2159| | | |
2160| Section | no. calls | wall time | % of total |
2161+---------------------------------+-----------+------------+------------+
2162| Assemble right-hand side | 1 | 0.000374s | 2.3% |
2163| Estimate | 1 | 0.000724s | 4.4% |
2164| Output results | 1 | 0.00277s | 17% |
2165| Setup | 1 | 0.00225s | 14% |
2166| Setup multigrid | 1 | 0.00181s | 11% |
2167| Solve | 1 | 0.00364s | 22% |
2168| Solve: 1 multigrid V-cycle | 1 | 0.000354s | 2.2% |
2169| Solve: CG | 1 | 0.00151s | 9.3% |
2170| Solve: Preconditioner setup | 1 | 0.00125s | 7.7% |
2171+---------------------------------+-----------+------------+------------+
2172
2173Cycle 1:
2174 Number of active cells: 24 (3 global levels)
2175 Partition efficiency: 0.276786
2176 Number of degrees of freedom: 139 (by level: 21, 65, 99)
2177 Number of CG iterations: 10
2178 Global error estimate: 0.216726
2179 Wrote solution_01.pvtu
2180
2181
2182+---------------------------------------------+------------+------------+
2183| Total wallclock time elapsed since start | 0.0169s | |
2184| | | |
2185| Section | no. calls | wall time | % of total |
2186+---------------------------------+-----------+------------+------------+
2187| Assemble right-hand side | 1 | 0.000309s | 1.8% |
2188| Estimate | 1 | 0.00156s | 9.2% |
2189| Output results | 1 | 0.00222s | 13% |
2190| Refine grid | 1 | 0.00278s | 16% |
2191| Setup | 1 | 0.00196s | 12% |
2192| Setup multigrid | 1 | 0.0023s | 14% |
2193| Solve | 1 | 0.00565s | 33% |
2194| Solve: 1 multigrid V-cycle | 1 | 0.000349s | 2.1% |
2195| Solve: CG | 1 | 0.00285s | 17% |
2196| Solve: Preconditioner setup | 1 | 0.00195s | 12% |
2197+---------------------------------+-----------+------------+------------+
2198
2199Cycle 2:
2200 Number of active cells: 51 (4 global levels)
2201 Partition efficiency: 0.41875
2202 Number of degrees of freedom: 245 (by level: 21, 65, 225, 25)
2203 Number of CG iterations: 11
2204 Global error estimate: 0.112098
2205 Wrote solution_02.pvtu
2206
2207
2208+---------------------------------------------+------------+------------+
2209| Total wallclock time elapsed since start | 0.0183s | |
2210| | | |
2211| Section | no. calls | wall time | % of total |
2212+---------------------------------+-----------+------------+------------+
2213| Assemble right-hand side | 1 | 0.000274s | 1.5% |
2214| Estimate | 1 | 0.00127s | 6.9% |
2215| Output results | 1 | 0.00227s | 12% |
2216| Refine grid | 1 | 0.0024s | 13% |
2217| Setup | 1 | 0.00191s | 10% |
2218| Setup multigrid | 1 | 0.00295s | 16% |
2219| Solve | 1 | 0.00702s | 38% |
2220| Solve: 1 multigrid V-cycle | 1 | 0.000398s | 2.2% |
2221| Solve: CG | 1 | 0.00376s | 21% |
2222| Solve: Preconditioner setup | 1 | 0.00238s | 13% |
2223+---------------------------------+-----------+------------+------------+
2224.
2225.
2226.
2227@endcode
2228Here, the timing of the `solve()` function is split up in 3 parts: setting
2229up the multigrid preconditioner, execution of a single multigrid V-cycle, and
2230the CG solver. The V-cycle that is timed is unnecessary for the overall solve
2231and only meant to give an insight at the different costs for AMG and GMG.
2232Also it should be noted that when using the AMG solver, "Workload imbalance"
2233is not included in the output since the hierarchy of coarse meshes is not
2234required.
2235
2236All results in this section are gathered on Intel Xeon Platinum 8280 (Cascade
2237Lake) nodes which have 56 cores and 192GB per node and support AVX-512 instructions,
2238allowing for vectorization over 8 doubles (vectorization used only in the matrix-free
2239computations). The code is compiled using gcc 7.1.0 with intel-mpi 17.0.3. Trilinos
224012.10.1 is used for the matrix-based GMG/AMG computations.
2241
2242We can then gather a variety of information by calling the program
2243with the input files that are provided in the directory in which
2244@ref step_50 "step-50" is located. Using these, and adjusting the number of mesh
2245refinement steps, we can produce information about how well the
2246program scales.
2247
2248The following table gives weak scaling timings for this program on up to 256M DoFs
2249and 7,168 processors. (Recall that weak scaling keeps the number of
2250degrees of freedom per processor constant while increasing the number of
2251processors; i.e., it considers larger and larger problems.)
2252Here, @f$\mathbb{E}@f$ is the partition efficiency from the
2253 introduction (also equal to 1.0/workload imbalance), "Setup" is a combination
2254of setup, setup multigrid, assemble, and assemble multigrid from the timing blocks,
2255and "Prec" is the preconditioner setup. Ideally all times would stay constant
2256over each problem size for the individual solvers, but since the partition
2257efficiency decreases from 0.371 to 0.161 from largest to smallest problem size,
2258we expect to see an approximately @f$0.371/0.161=2.3@f$ times increase in timings
2259for GMG. This is, in fact, pretty close to what we really get:
2260
2261<table align="center" class="doxtable">
2262<tr>
2263 <th colspan="4"></th>
2264 <th></th>
2265 <th colspan="4">MF-GMG</th>
2266 <th></th>
2267 <th colspan="4">MB-GMG</th>
2268 <th></th>
2269 <th colspan="4">AMG</th>
2270</tr>
2271<tr>
2272 <th align="right">Procs</th>
2273 <th align="right">Cycle</th>
2274 <th align="right">DoFs</th>
2275 <th align="right">@f$\mathbb{E}@f$</th>
2276 <th></th>
2277 <th align="right">Setup</th>
2278 <th align="right">Prec</th>
2279 <th align="right">Solve</th>
2280 <th align="right">Total</th>
2281 <th></th>
2282 <th align="right">Setup</th>
2283 <th align="right">Prec</th>
2284 <th align="right">Solve</th>
2285 <th align="right">Total</th>
2286 <th></th>
2287 <th align="right">Setup</th>
2288 <th align="right">Prec</th>
2289 <th align="right">Solve</th>
2290 <th align="right">Total</th>
2291</tr>
2292<tr>
2293 <td align="right">112</th>
2294 <td align="right">13</th>
2295 <td align="right">4M</th>
2296 <td align="right">0.37</th>
2297 <td></td>
2298 <td align="right">0.742</th>
2299 <td align="right">0.393</th>
2300 <td align="right">0.200</th>
2301 <td align="right">1.335</th>
2302 <td></td>
2303 <td align="right">1.714</th>
2304 <td align="right">2.934</th>
2305 <td align="right">0.716</th>
2306 <td align="right">5.364</th>
2307 <td></td>
2308 <td align="right">1.544</th>
2309 <td align="right">0.456</th>
2310 <td align="right">1.150</th>
2311 <td align="right">3.150</th>
2312</tr>
2313<tr>
2314 <td align="right">448</th>
2315 <td align="right">15</th>
2316 <td align="right">16M</th>
2317 <td align="right">0.29</th>
2318 <td></td>
2319 <td align="right">0.884</th>
2320 <td align="right">0.535</th>
2321 <td align="right">0.253</th>
2322 <td align="right">1.672</th>
2323 <td></td>
2324 <td align="right">1.927</th>
2325 <td align="right">3.776</th>
2326 <td align="right">1.190</th>
2327 <td align="right">6.893</th>
2328 <td></td>
2329 <td align="right">1.544</th>
2330 <td align="right">0.456</th>
2331 <td align="right">1.150</th>
2332 <td align="right">3.150</th>
2333</tr>
2334<tr>
2335 <td align="right">1,792</th>
2336 <td align="right">17</th>
2337 <td align="right">65M</th>
2338 <td align="right">0.22</th>
2339 <td></td>
2340 <td align="right">1.122</th>
2341 <td align="right">0.686</th>
2342 <td align="right">0.309</th>
2343 <td align="right">2.117</th>
2344 <td></td>
2345 <td align="right">2.171</th>
2346 <td align="right">4.862</th>
2347 <td align="right">1.660</th>
2348 <td align="right">8.693</th>
2349 <td></td>
2350 <td align="right">1.654</th>
2351 <td align="right">0.546</th>
2352 <td align="right">1.460</th>
2353 <td align="right">3.660</th>
2354</tr>
2355<tr>
2356 <td align="right">7,168</th>
2357 <td align="right">19</th>
2358 <td align="right">256M</th>
2359 <td align="right">0.16</th>
2360 <td></td>
2361 <td align="right">1.214</th>
2362 <td align="right">0.893</th>
2363 <td align="right">0.521</th>
2364 <td align="right">2.628</th>
2365 <td></td>
2366 <td align="right">2.386</th>
2367 <td align="right">7.260</th>
2368 <td align="right">2.560</th>
2369 <td align="right">12.206</th>
2370 <td></td>
2371 <td align="right">1.844</th>
2372 <td align="right">1.010</th>
2373 <td align="right">1.890</th>
2374 <td align="right">4.744</th>
2375</tr>
2376</table>
2377
2378On the other hand, the algebraic multigrid in the last set of columns
2379is relatively unaffected by the increasing imbalance of the mesh
2380hierarchy (because it doesn't use the mesh hierarchy) and the growth
2381in time is rather driven by other factors that are well documented in
2382the literature (most notably that the algorithmic complexity of
2383some parts of algebraic multigrid methods appears to be @f${\cal O}(N
2384\log N)@f$ instead of @f${\cal O}(N)@f$ for geometric multigrid).
2385
2386The upshort of the table above is that the matrix-free geometric multigrid
2387method appears to be the fastest approach to solving this equation if
2388not by a huge margin. Matrix-based methods, on the other hand, are
2389consistently the worst.
2390
2391The following figure provides strong scaling results for each method, i.e.,
2392we solve the same problem on more and more processors. Specifically,
2393we consider the problems after 16 mesh refinement cycles
2394(32M DoFs) and 19 cycles (256M DoFs), on between 56 to 28,672 processors:
2395
2396<img width="600px" src="https://www.dealii.org/images/steps/developer/step-50-strong-scaling.png" alt="">
2397
2398While the matrix-based GMG solver and AMG scale similarly and have a
2399similar time to solution (at least as long as there is a substantial
2400number of unknowns per processor -- say, several 10,000), the
2401matrix-free GMG solver scales much better and solves the finer problem
2402in roughly the same time as the AMG solver for the coarser mesh with
2403only an eighth of the number of processors. Conversely, it can solve the
2404same problem on the same number of processors in about one eighth the
2405time.
2406
2407
2408<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
2409
2410
2411<a name="Testingconvergenceandhigherorderelements"></a><h4> Testing convergence and higher order elements </h4>
2412
2413
2414The finite element degree is currently hard-coded as 2, see the template
2415arguments of the main class. It is easy to change. To test, it would be
2416interesting to switch to a test problem with a reference solution. This way,
2417you can compare error rates.
2418
2419<a name="Coarsesolver"></a><h4> Coarse solver </h4>
2420
2421
2422A more interesting example would involve a more complicated coarse mesh (see
2423@ref step_49 "step-49" for inspiration). The issue in that case is that the coarsest
2424level of the mesh hierarchy is actually quite large, and one would
2425have to think about ways to solve the coarse level problem
2426efficiently. (This is not an issue for algebraic multigrid methods
2427because they would just continue to build coarser and coarser levels
2428of the matrix, regardless of their geometric origin.)
2429
2430In the program here, we simply solve the coarse level problem with a
2431Conjugate Gradient method without any preconditioner. That is acceptable
2432if the coarse problem is really small -- for example, if the coarse
2433mesh had a single cell, then the coarse mesh problems has a @f$9\times 9@f$
2434matrix in 2d, and a @f$27\times 27@f$ matrix in 3d; for the coarse mesh we
2435use on the @f$L@f$-shaped domain of the current program, these sizes are
2436@f$21\times 21@f$ in 2d and @f$117\times 117@f$ in 3d. But if the coarse mesh
2437consists of hundreds or thousands of cells, this approach will no
2438longer work and might start to dominate the overall run-time of each V-cycle.
2439A common approach is then to solve the coarse mesh problem using an
2440algebraic multigrid preconditioner; this would then, however, require
2441assembling the coarse matrix (even for the matrix-free version) as
2442input to the AMG implementation.
2443 *
2444 *
2445<a name="PlainProg"></a>
2446<h1> The plain program</h1>
2447@include "step-50.cc"
2448*/
void reinit(const IndexSet &local_constraints=IndexSet())
Definition fe_q.h:551
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
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 AssertThrow(cond, exc)
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)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_relevant_level_dofs(const DoFHandler< dim, spacedim > &dof_handler, const unsigned int level)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
@ matrix
Contents is actually a matrix.
PETScWrappers::PreconditionBoomerAMG PreconditionAMG
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:75
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
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
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)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
void free(T *&pointer)
Definition cuda.h:97
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
std::string compress(const std::string &input)
Definition utilities.cc:390
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)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:13826
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
Definition types.h:33
unsigned int boundary_id
Definition types.h:141
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
TasksParallelScheme tasks_parallel_scheme
const TriangulationDescription::Settings settings