deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+00:00
\(\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-26.h
Go to the documentation of this file.
1,
680 *   const unsigned int component) const
681 *   {
682 *   (void)component;
683 *   Assert(component == 0, ExcIndexRange(component, 0, 1));
684 *   return 0;
685 *   }
686 *  
687 *  
688 *  
689 * @endcode
690 *
691 *
692 * <a name="step_26-ThecodeHeatEquationcodeimplementation"></a>
693 * <h3>The <code>HeatEquation</code> implementation</h3>
694 *
695
696 *
697 * It is time now for the implementation of the main class. Let's
698 * start with the constructor which selects a linear element, a time
699 * step constant at 1/500 (remember that one period of the source
700 * on the right hand side was set to 0.2 above, so we resolve each
701 * period with 100 time steps) and chooses the Crank Nicolson method
702 * by setting @f$\theta=1/2@f$.
703 *
704 * @code
705 *   template <int dim>
706 *   HeatEquation<dim>::HeatEquation()
707 *   : fe(1)
708 *   , dof_handler(triangulation)
709 *   , time_step(1. / 500)
710 *   , theta(0.5)
711 *   {}
712 *  
713 *  
714 *  
715 * @endcode
716 *
717 *
718 * <a name="step_26-codeHeatEquationsetup_systemcode"></a>
719 * <h4><code>HeatEquation::setup_system</code></h4>
720 *
721
722 *
723 * The next function is the one that sets up the DoFHandler object,
724 * computes the constraints, and sets the linear algebra objects
725 * to their correct sizes. We also compute the mass and Laplace
726 * matrix here by simply calling two functions in the library.
727 *
728
729 *
730 * Note that we do not take the hanging node constraints into account when
731 * assembling the matrices (both functions have an AffineConstraints argument
732 * that defaults to an empty object). This is because we are going to
733 * condense the constraints in run() after combining the matrices for the
734 * current time-step.
735 *
736 * @code
737 *   template <int dim>
738 *   void HeatEquation<dim>::setup_system()
739 *   {
740 *   dof_handler.distribute_dofs(fe);
741 *  
742 *   std::cout << std::endl
743 *   << "===========================================" << std::endl
744 *   << "Number of active cells: " << triangulation.n_active_cells()
745 *   << std::endl
746 *   << "Number of degrees of freedom: " << dof_handler.n_dofs()
747 *   << std::endl
748 *   << std::endl;
749 *  
750 *   constraints.clear();
751 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
752 *   constraints.close();
753 *  
754 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
755 *   DoFTools::make_sparsity_pattern(dof_handler,
756 *   dsp,
757 *   constraints,
758 *   /*keep_constrained_dofs = */ true);
759 *   sparsity_pattern.copy_from(dsp);
760 *  
761 *   mass_matrix.reinit(sparsity_pattern);
762 *   laplace_matrix.reinit(sparsity_pattern);
763 *   system_matrix.reinit(sparsity_pattern);
764 *  
765 *   MatrixCreator::create_mass_matrix(dof_handler,
766 *   QGauss<dim>(fe.degree + 1),
767 *   mass_matrix);
768 *   MatrixCreator::create_laplace_matrix(dof_handler,
769 *   QGauss<dim>(fe.degree + 1),
770 *   laplace_matrix);
771 *  
772 *   solution.reinit(dof_handler.n_dofs());
773 *   old_solution.reinit(dof_handler.n_dofs());
774 *   system_rhs.reinit(dof_handler.n_dofs());
775 *   }
776 *  
777 *  
778 * @endcode
779 *
780 *
781 * <a name="step_26-codeHeatEquationsolve_time_stepcode"></a>
782 * <h4><code>HeatEquation::solve_time_step</code></h4>
783 *
784
785 *
786 * The next function is the one that solves the actual linear system
787 * for a single time step. There is nothing surprising here:
788 *
789 * @code
790 *   template <int dim>
791 *   void HeatEquation<dim>::solve_time_step()
792 *   {
793 *   SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
794 *   SolverCG<Vector<double>> cg(solver_control);
795 *  
796 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
797 *   preconditioner.initialize(system_matrix, 1.0);
798 *  
799 *   cg.solve(system_matrix, solution, system_rhs, preconditioner);
800 *  
801 *   constraints.distribute(solution);
802 *  
803 *   std::cout << " " << solver_control.last_step() << " CG iterations."
804 *   << std::endl;
805 *   }
806 *  
807 *  
808 *  
809 * @endcode
810 *
811 *
812 * <a name="step_26-codeHeatEquationoutput_resultscode"></a>
813 * <h4><code>HeatEquation::output_results</code></h4>
814 *
815
816 *
817 * Neither is there anything new in generating graphical output other than the
818 * fact that we tell the DataOut object what the current time and time step
819 * number is, so that this can be written into the output file :
820 *
821 * @code
822 *   template <int dim>
823 *   void HeatEquation<dim>::output_results() const
824 *   {
825 *   DataOut<dim> data_out;
826 *  
827 *   data_out.attach_dof_handler(dof_handler);
828 *   data_out.add_data_vector(solution, "U");
829 *  
830 *   data_out.build_patches();
831 *  
832 *   data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
833 *  
834 *   const std::string filename =
835 *   "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtk";
836 *   std::ofstream output(filename);
837 *   data_out.write_vtk(output);
838 *   }
839 *  
840 *  
841 * @endcode
842 *
843 *
844 * <a name="step_26-codeHeatEquationrefine_meshcode"></a>
845 * <h4><code>HeatEquation::refine_mesh</code></h4>
846 *
847
848 *
849 * This function is the interesting part of the program. It takes care of
850 * the adaptive mesh refinement. The three tasks
851 * this function performs is to first find out which cells to
852 * refine/coarsen, then to actually do the refinement and eventually
853 * transfer the solution vectors between the two different grids. The first
854 * task is simply achieved by using the well-established Kelly error
855 * estimator on the solution. The second task is to actually do the
856 * remeshing. That involves only basic functions as well, such as the
857 * <code>refine_and_coarsen_fixed_fraction</code> that refines those cells
858 * with the largest estimated error that together make up 60 per cent of the
859 * error, and coarsens those cells with the smallest error that make up for
860 * a combined 40 per cent of the error. Note that for problems such as the
861 * current one where the areas where something is going on are shifting
862 * around, we want to aggressively coarsen so that we can move cells
863 * around to where it is necessary.
864 *
865
866 *
867 * As already discussed in the introduction, too small a mesh leads to
868 * too small a time step, whereas too large a mesh leads to too little
869 * resolution. Consequently, after the first two steps, we have two
870 * loops that limit refinement and coarsening to an allowable range of
871 * cells:
872 *
873 * @code
874 *   template <int dim>
875 *   void HeatEquation<dim>::refine_mesh(const unsigned int min_grid_level,
876 *   const unsigned int max_grid_level)
877 *   {
878 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
879 *  
880 *   KellyErrorEstimator<dim>::estimate(
881 *   dof_handler,
882 *   QGauss<dim - 1>(fe.degree + 1),
883 *   std::map<types::boundary_id, const Function<dim> *>(),
884 *   solution,
885 *   estimated_error_per_cell);
886 *  
887 *   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
888 *   estimated_error_per_cell,
889 *   0.6,
890 *   0.4);
891 *  
892 *   if (triangulation.n_levels() > max_grid_level)
893 *   for (const auto &cell :
894 *   triangulation.active_cell_iterators_on_level(max_grid_level))
895 *   cell->clear_refine_flag();
896 *   for (const auto &cell :
897 *   triangulation.active_cell_iterators_on_level(min_grid_level))
898 *   cell->clear_coarsen_flag();
899 * @endcode
900 *
901 * These two loops above are slightly different but this is easily
902 * explained. In the first loop, instead of calling
903 * <code>triangulation.end()</code> we may as well have called
904 * <code>triangulation.end_active(max_grid_level)</code>. The two
905 * calls should yield the same iterator since iterators are sorted
906 * by level and there should not be any cells on levels higher than
907 * on level <code>max_grid_level</code>. In fact, this very piece
908 * of code makes sure that this is the case.
909 *
910
911 *
912 * As part of mesh refinement we need to transfer the solution vectors
913 * from the old mesh to the new one. To this end we use the
914 * SolutionTransfer class and we have to prepare the solution vectors that
915 * should be transferred to the new grid (we will lose the old grid once
916 * we have done the refinement so the transfer has to happen concurrently
917 * with refinement). At the point where we call this function, we will
918 * have just computed the solution, so we no longer need the old_solution
919 * variable (it will be overwritten by the solution just after the mesh
920 * may have been refined, i.e., at the end of the time step; see below).
921 * In other words, we only need the one solution vector, and we copy it
922 * to a temporary object where it is safe from being reset when we further
923 * down below call <code>setup_system()</code>.
924 *
925
926 *
927 * Consequently, we initialize a SolutionTransfer object by attaching
928 * it to the old DoF handler. We then prepare the triangulation and the
929 * data vector for refinement (in this order).
930 *
931 * @code
932 *   SolutionTransfer<dim> solution_trans(dof_handler);
933 *  
934 *   Vector<double> previous_solution;
935 *   previous_solution = solution;
936 *   triangulation.prepare_coarsening_and_refinement();
937 *   solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
938 *  
939 * @endcode
940 *
941 * Now everything is ready, so do the refinement and recreate the DoF
942 * structure on the new grid, and finally initialize the matrix structures
943 * and the new vectors in the <code>setup_system</code> function. Next, we
944 * actually perform the interpolation of the solution from old to new
945 * grid. The final step is to apply the hanging node constraints to the
946 * solution vector, i.e., to make sure that the values of degrees of
947 * freedom located on hanging nodes are so that the solution is
948 * continuous. This is necessary since SolutionTransfer only operates on
949 * cells locally, without regard to the neighborhood.
950 *
951 * @code
952 *   triangulation.execute_coarsening_and_refinement();
953 *   setup_system();
954 *  
955 *   solution_trans.interpolate(solution);
956 *   constraints.distribute(solution);
957 *   }
958 *  
959 *  
960 *  
961 * @endcode
962 *
963 *
964 * <a name="step_26-codeHeatEquationruncode"></a>
965 * <h4><code>HeatEquation::run</code></h4>
966 *
967
968 *
969 * This is the main driver of the program, where we loop over all
970 * time steps. At the top of the function, we set the number of
971 * initial global mesh refinements and the number of initial cycles of
972 * adaptive mesh refinement by repeating the first time step a few
973 * times. Then we create a mesh, initialize the various objects we will
974 * work with, set a label for where we should start when re-running
975 * the first time step, and interpolate the initial solution onto
976 * out mesh (we choose the zero function here, which of course we could
977 * do in a simpler way by just setting the solution vector to zero). We
978 * also output the initial time step once.
979 *
980
981 *
982 * @note If you're an experienced programmer, you may be surprised
983 * that we use a <code>goto</code> statement in this piece of code!
984 * <code>goto</code> statements are not particularly well liked any
985 * more since Edsgar Dijkstra, one of the greats of computer science,
986 * wrote a letter in 1968 called "Go To Statement considered harmful"
987 * (see <a href="http://en.wikipedia.org/wiki/Considered_harmful">here</a>).
988 * The author of this code subscribes to this notion whole-heartedly:
989 * <code>goto</code> is hard to understand. In fact, deal.II contains
990 * virtually no occurrences: excluding code that was essentially
991 * transcribed from books and not counting duplicated code pieces,
992 * there are 3 locations in about 600,000 lines of code at the time
993 * this note is written; we also use it in 4 tutorial programs, in
994 * exactly the same context as here. Instead of trying to justify
995 * the occurrence here, let's first look at the code and we'll come
996 * back to the issue at the end of function.
997 *
998 * @code
999 *   template <int dim>
1000 *   void HeatEquation<dim>::run()
1001 *   {
1002 *   const unsigned int initial_global_refinement = 2;
1003 *   const unsigned int n_adaptive_pre_refinement_steps = 4;
1004 *  
1006 *   triangulation.refine_global(initial_global_refinement);
1007 *  
1008 *   setup_system();
1009 *  
1010 *   unsigned int pre_refinement_step = 0;
1011 *  
1012 *   Vector<double> tmp;
1013 *   Vector<double> forcing_terms;
1014 *  
1015 *   start_time_iteration:
1016 *  
1017 *   time = 0.0;
1018 *   timestep_number = 0;
1019 *  
1020 *   tmp.reinit(solution.size());
1021 *   forcing_terms.reinit(solution.size());
1022 *  
1023 *  
1024 *   VectorTools::interpolate(dof_handler,
1026 *   old_solution);
1027 *   solution = old_solution;
1028 *  
1029 *   output_results();
1030 *  
1031 * @endcode
1032 *
1033 * Then we start the main loop until the computed time exceeds our
1034 * end time of 0.5. The first task is to build the right hand
1035 * side of the linear system we need to solve in each time step.
1036 * Recall that it contains the term @f$MU^{n-1}-(1-\theta)k_n AU^{n-1}@f$.
1037 * We put these terms into the variable system_rhs, with the
1038 * help of a temporary vector:
1039 *
1040 * @code
1041 *   const double end_time = 0.5;
1042 *   while (time <= end_time)
1043 *   {
1044 *   time += time_step;
1045 *   ++timestep_number;
1046 *  
1047 *   std::cout << "Time step " << timestep_number << " at t=" << time
1048 *   << std::endl;
1049 *  
1050 *   mass_matrix.vmult(system_rhs, old_solution);
1051 *  
1052 *   laplace_matrix.vmult(tmp, old_solution);
1053 *   system_rhs.add(-(1 - theta) * time_step, tmp);
1054 *  
1055 * @endcode
1056 *
1057 * The second piece is to compute the contributions of the source
1058 * terms. This corresponds to the term @f$k_n
1059 * \left[ (1-\theta)F^{n-1} + \theta F^n \right]@f$. The following
1060 * code calls VectorTools::create_right_hand_side to compute the
1061 * vectors @f$F@f$, where we set the time of the right hand side
1062 * (source) function before we evaluate it. The result of this
1063 * all ends up in the forcing_terms variable:
1064 *
1065 * @code
1066 *   RightHandSide<dim> rhs_function;
1067 *   rhs_function.set_time(time);
1068 *   VectorTools::create_right_hand_side(dof_handler,
1069 *   QGauss<dim>(fe.degree + 1),
1070 *   rhs_function,
1071 *   tmp);
1072 *   forcing_terms = tmp;
1073 *   forcing_terms *= time_step * theta;
1074 *  
1075 *   rhs_function.set_time(time - time_step);
1076 *   VectorTools::create_right_hand_side(dof_handler,
1077 *   QGauss<dim>(fe.degree + 1),
1078 *   rhs_function,
1079 *   tmp);
1080 *  
1081 *   forcing_terms.add(time_step * (1 - theta), tmp);
1082 *  
1083 * @endcode
1084 *
1085 * Next, we add the forcing terms to the ones that
1086 * come from the time stepping, and also build the matrix
1087 * @f$M+k_n\theta A@f$ that we have to invert in each time step.
1088 * The final piece of these operations is to eliminate
1089 * hanging node constrained degrees of freedom from the
1090 * linear system:
1091 *
1092 * @code
1093 *   system_rhs += forcing_terms;
1094 *  
1095 *   system_matrix.copy_from(mass_matrix);
1096 *   system_matrix.add(theta * time_step, laplace_matrix);
1097 *  
1098 *   constraints.condense(system_matrix, system_rhs);
1099 *  
1100 * @endcode
1101 *
1102 * There is one more operation we need to do before we
1103 * can solve it: boundary values. To this end, we create
1104 * a boundary value object, set the proper time to the one
1105 * of the current time step, and evaluate it as we have
1106 * done many times before. The result is used to also
1107 * set the correct boundary values in the linear system:
1108 *
1109 * @code
1110 *   {
1111 *   BoundaryValues<dim> boundary_values_function;
1112 *   boundary_values_function.set_time(time);
1113 *  
1114 *   std::map<types::global_dof_index, double> boundary_values;
1116 *   0,
1117 *   boundary_values_function,
1118 *   boundary_values);
1119 *  
1120 *   MatrixTools::apply_boundary_values(boundary_values,
1121 *   system_matrix,
1122 *   solution,
1123 *   system_rhs);
1124 *   }
1125 *  
1126 * @endcode
1127 *
1128 * With this out of the way, all we have to do is solve the
1129 * system, generate graphical data, and...
1130 *
1131 * @code
1132 *   solve_time_step();
1133 *  
1134 *   output_results();
1135 *  
1136 * @endcode
1137 *
1138 * ...take care of mesh refinement. Here, what we want to do is
1139 * (i) refine the requested number of times at the very beginning
1140 * of the solution procedure, after which we jump to the top to
1141 * restart the time iteration, (ii) refine every fifth time
1142 * step after that.
1143 *
1144
1145 *
1146 * The time loop and, indeed, the main part of the program ends
1147 * with starting into the next time step by setting old_solution
1148 * to the solution we have just computed.
1149 *
1150 * @code
1151 *   if ((timestep_number == 1) &&
1152 *   (pre_refinement_step < n_adaptive_pre_refinement_steps))
1153 *   {
1154 *   refine_mesh(initial_global_refinement,
1155 *   initial_global_refinement +
1156 *   n_adaptive_pre_refinement_steps);
1157 *   ++pre_refinement_step;
1158 *  
1159 *   std::cout << std::endl;
1160 *  
1161 *   goto start_time_iteration;
1162 *   }
1163 *   else if ((timestep_number > 0) && (timestep_number % 5 == 0))
1164 *   {
1165 *   refine_mesh(initial_global_refinement,
1166 *   initial_global_refinement +
1167 *   n_adaptive_pre_refinement_steps);
1168 *   tmp.reinit(solution.size());
1169 *   forcing_terms.reinit(solution.size());
1170 *   }
1171 *  
1172 *   old_solution = solution;
1173 *   }
1174 *   }
1175 *   } // namespace Step26
1176 * @endcode
1177 *
1178 * Now that you have seen what the function does, let us come back to the issue
1179 * of the <code>goto</code>. In essence, what the code does is
1180 * something like this:
1181 * <div class=CodeFragmentInTutorialComment>
1182 * @code
1183 * void run ()
1184 * {
1185 * initialize;
1186 * start_time_iteration:
1187 * for (timestep=1...)
1188 * {
1189 * solve timestep;
1190 * if (timestep==1 && not happy with the result)
1191 * {
1192 * adjust some data structures;
1193 * goto start_time_iteration; // simply try again
1194 * }
1195 * postprocess;
1196 * }
1197 * }
1198 * @endcode
1199 * </div>
1200 * Here, the condition "happy with the result" is whether we'd like to keep
1201 * the current mesh or would rather refine the mesh and start over on the
1202 * new mesh. We could of course replace the use of the <code>goto</code>
1203 * by the following:
1204 * <div class=CodeFragmentInTutorialComment>
1205 * @code
1206 * void run ()
1207 * {
1208 * initialize;
1209 * while (true)
1210 * {
1211 * solve timestep;
1212 * if (not happy with the result)
1213 * adjust some data structures;
1214 * else
1215 * break;
1216 * }
1217 * postprocess;
1218 *
1219
1220 * for (timestep=2...)
1221 * {
1222 * solve timestep;
1223 * postprocess;
1224 * }
1225 * }
1226 * @endcode
1227 * </div>
1228 * This has the advantage of getting rid of the <code>goto</code>
1229 * but the disadvantage of having to duplicate the code that implements
1230 * the "solve timestep" and "postprocess" operations in two different
1231 * places. This could be countered by putting these parts of the code
1232 * (sizable chunks in the actual implementation above) into their
1233 * own functions, but a <code>while(true)</code> loop with a
1234 * <code>break</code> statement is not really all that much easier
1235 * to read or understand than a <code>goto</code>.
1236 *
1237
1238 *
1239 * In the end, one might simply agree that <i>in general</i>
1240 * <code>goto</code> statements are a bad idea but be pragmatic and
1241 * state that there may be occasions where they can help avoid code
1242 * duplication and awkward control flow. This may be one of these
1243 * places, and it matches the position Steve McConnell takes in his
1244 * excellent book "Code Complete" @cite CodeComplete about good
1245 * programming practices (see the mention of this book in the
1246 * introduction of @ref step_1 "step-1") that spends a surprising ten pages on the
1247 * question of <code>goto</code> in general.
1248 *
1249
1250 *
1251 *
1252
1253 *
1254 *
1255 * <a name="step_26-Thecodemaincodefunction"></a>
1256 * <h3>The <code>main</code> function</h3>
1257 *
1258
1259 *
1260 * Having made it this far, there is, again, nothing
1261 * much to discuss for the main function of this
1262 * program: it looks like all such functions since @ref step_6 "step-6".
1263 *
1264 * @code
1265 *   int main()
1266 *   {
1267 *   try
1268 *   {
1269 *   using namespace Step26;
1270 *  
1271 *   HeatEquation<2> heat_equation_solver;
1272 *   heat_equation_solver.run();
1273 *   }
1274 *   catch (std::exception &exc)
1275 *   {
1276 *   std::cerr << std::endl
1277 *   << std::endl
1278 *   << "----------------------------------------------------"
1279 *   << std::endl;
1280 *   std::cerr << "Exception on processing: " << std::endl
1281 *   << exc.what() << std::endl
1282 *   << "Aborting!" << std::endl
1283 *   << "----------------------------------------------------"
1284 *   << std::endl;
1285 *  
1286 *   return 1;
1287 *   }
1288 *   catch (...)
1289 *   {
1290 *   std::cerr << std::endl
1291 *   << std::endl
1292 *   << "----------------------------------------------------"
1293 *   << std::endl;
1294 *   std::cerr << "Unknown exception!" << std::endl
1295 *   << "Aborting!" << std::endl
1296 *   << "----------------------------------------------------"
1297 *   << std::endl;
1298 *   return 1;
1299 *   }
1300 *  
1301 *   return 0;
1302 *   }
1303 * @endcode
1304<a name="step_26-Results"></a><h1>Results</h1>
1305
1306
1307As in many of the tutorials, the actual output of the program matters less
1308than how we arrived there. Nonetheless, here it is:
1309@code
1310===========================================
1311Number of active cells: 48
1312Number of degrees of freedom: 65
1313
1314Time step 1 at t=0.002
1315 7 CG iterations.
1316
1317===========================================
1318Number of active cells: 60
1319Number of degrees of freedom: 81
1320
1321
1322Time step 1 at t=0.002
1323 7 CG iterations.
1324
1325===========================================
1326Number of active cells: 105
1327Number of degrees of freedom: 136
1328
1329
1330Time step 1 at t=0.002
1331 7 CG iterations.
1332
1333[...]
1334
1335Time step 249 at t=0.498
1336 13 CG iterations.
1337Time step 250 at t=0.5
1338 14 CG iterations.
1339
1340===========================================
1341Number of active cells: 1803
1342Number of degrees of freedom: 2109
1343@endcode
1344
1345Maybe of more interest is a visualization of the solution and the mesh on which
1346it was computed:
1347
1348<img src="https://www.dealii.org/images/steps/developer/step-26.movie.gif" alt="Animation of the solution of step 26.">
1349
1350The movie shows how the two sources switch on and off and how the mesh reacts
1351to this. It is quite obvious that the mesh as is is probably not the best we
1352could come up with. We'll get back to this in the next section.
1353
1354
1355<a name="step-26-extensions"></a>
1356<a name="step_26-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1357
1358
1359There are at least two areas where one can improve this program significantly:
1360adaptive time stepping and a better choice of the mesh.
1361
1362<a name="step_26-Adaptivetimestepping"></a><h4>Adaptive time stepping</h4>
1363
1364
1365Having chosen an implicit time stepping scheme, we are not bound by any
1366CFL-like condition on the time step. Furthermore, because the time scales on
1367which change happens on a given cell in the heat equation are not bound to the
1368cells diameter (unlike the case with the wave equation, where we had a fixed
1369speed of information transport that couples the temporal and spatial scales),
1370we can choose the time step as we please. Or, better, choose it as we deem
1371necessary for accuracy.
1372
1373Looking at the solution, it is clear that the action does not happen uniformly
1374over time: a lot is changing around the times when we switch on a source, things
1375become less dramatic once a source is on for a little while, and we enter a
1376long phase of decline when both sources are off. During these times, we could
1377surely get away with a larger time step than before without sacrificing too
1378much accuracy.
1379
1380The literature has many suggestions on how to choose the time step size
1381adaptively. Much can be learned, for example, from the way ODE solvers choose
1382their time steps. One can also be inspired by a posteriori error estimators
1383that can, ideally, be written in a way that they consist of a temporal and a
1384spatial contribution to the overall error. If the temporal one is too large,
1385we should choose a smaller time step. Ideas in this direction can be found,
1386for example, in the PhD thesis of a former principal developer of deal.II,
1387Ralf Hartmann, published by the University of Heidelberg, Germany, in 2002
1388(see @cite Har02).
1389
1390
1391<a name="step_26-Bettertimesteppingmethods"></a><h4>Better time stepping methods</h4>
1392
1393
1394We here use one of the simpler time stepping methods, namely the second order
1395in time Crank-Nicolson method. This is surely already better than the even
1396more widely used (and even less accurate) implicit (=backward) Euler method,
1397but many other, more accurate methods such as BDF or
1398Runge-Kutta methods are available and should be used as they do not represent
1399much additional effort. It is not difficult to implement this for the current
1400program, if one wanted; a more systematic treatment is also given in @ref step_86 "step-86".
1401
1402As a general rule, however, one should not be implementing time stepping
1403methods by hand, as we do here, for problems that do not require
1404exploiting special properties of the equation and consequently require
1405specialized time stepping methods. (The heat equation does not fall into
1406this category, and "standard" time stepping methods are all we need here.)
1407Rather, one should use one of the available
1408high-quality libraries for time stepping, for the same reasons as one should
1409not be implementing finite element methods by hand but use deal.II instead.
1410Indeed, deal.II has interfaces to two such time stepping library,
1411[SUNDIALS](https://computing.llnl.gov/projects/sundials) and the
1412[PETSc TS package](https://petsc.org/release/manualpages/TS/TS/), already available.
1413In particular, the SUNDIALS::ARKode class would make for a great starting
1414point for the use of much better (and much more accurate) time steppers,
1415as would PETScWrappers::TimeStepper;
1416the methods one would then get also implement the automatic time step
1417control mentioned above. To see how this works, take a look at
1418@ref step_86 "step-86".
1419
1420
1421<a name="step_26-Betterrefinementcriteria"></a><h4>Better refinement criteria</h4>
1422
1423
1424If you look at the meshes in the movie above, it is clear that they are not
1425particularly well suited to the task at hand. In fact, they look rather
1426random.
1427
1428There are two factors at play. First, there are some islands where cells
1429have been refined but that are surrounded by non-refined cells (and there
1430are probably also a few occasional coarsened islands). These are not terrible,
1431as they most of the time do not affect the approximation quality of the mesh,
1432but they also don't help because so many of their additional degrees of
1433freedom are in fact constrained by hanging node constraints. That said,
1434this is easy to fix: the Triangulation class takes an argument to its
1435constructor indicating a level of "mesh smoothing". Passing one of many
1436possible flags, this instructs the triangulation to refine some additional
1437cells, or not to refine some cells, so that the resulting mesh does not have
1438these artifacts.
1439
1440The second problem is more severe: the mesh appears to lag the solution.
1441The underlying reason is that we only adapt the mesh once every fifth
1442time step, and only allow for a single refinement in these cases. Whenever a
1443source switches on, the solution had been very smooth in this area before and
1444the mesh was consequently rather coarse. This implies that the next time step
1445when we refine the mesh, we will get one refinement level more in this area,
1446and five time steps later another level, etc. But this is not enough: first,
1447we should refine immediately when a source switches on (after all, in the
1448current context we at least know what the right hand side is), and we should
1449allow for more than one refinement level. Of course, all of this can be done
1450using deal.II, it just requires a bit of algorithmic thinking in how to make
1451this work!
1452
1453
1454<a name="step_26-Positivitypreservation"></a><h4>Positivity preservation</h4>
1455
1456
1457To increase the accuracy and resolution of your simulation in time, one
1458typically decreases the time step size @f$k_n@f$. If you start playing around
1459with the time step in this particular example, you will notice that the
1460solution becomes partly negative, if @f$k_n@f$ is below a certain threshold.
1461This is not what we would expect to happen (in nature).
1462
1463To get an idea of this behavior mathematically, let us consider a general,
1464fully discrete problem:
1465@f{align*}{
1466 A u^{n} = B u^{n-1}.
1467@f}
1468The general form of the @f$i@f$th equation then reads:
1469@f{align*}{
1470 a_{ii} u^{n}_i &= b_{ii} u^{n-1}_i +
1471 \sum\limits_{j \in S_i} \left( b_{ij} u^{n-1}_j - a_{ij} u^{n}_j \right),
1472@f}
1473where @f$S_i@f$ is the set of degrees of freedom that DoF @f$i@f$ couples with (i.e.,
1474for which either the matrix @f$A@f$ or matrix @f$B@f$ has a nonzero entry at position
1475@f$(i,j)@f$). If all coefficients
1476fulfill the following conditions:
1477@f{align*}{
1478 a_{ii} &> 0, & b_{ii} &\geq 0, & a_{ij} &\leq 0, & b_{ij} &\geq 0,
1479 &
1480 \forall j &\in S_i,
1481@f}
1482all solutions @f$u^{n}@f$ keep their sign from the previous ones @f$u^{n-1}@f$, and
1483consequently from the initial values @f$u^0@f$. See e.g.
1484<a href="http://bookstore.siam.org/cs14/">Kuzmin, H&auml;m&auml;l&auml;inen</a>
1485for more information on positivity preservation.
1486
1487Depending on the PDE to solve and the time integration scheme used, one is
1488able to deduce conditions for the time step @f$k_n@f$. For the heat equation with
1489the Crank-Nicolson scheme,
1490<a href="https://doi.org/10.2478/cmam-2010-0025">Schatz et. al.</a> have
1491translated it to the following ones:
1492@f{align*}{
1493 (1 - \theta) k a_{ii} &\leq m_{ii},\qquad \forall i,
1494 &
1495 \theta k \left| a_{ij} \right| &\geq m_{ij},\qquad j \neq i,
1496@f}
1497where @f$M = m_{ij}@f$ denotes the @ref GlossMassMatrix "mass matrix" and @f$A = a_{ij}@f$ the stiffness
1498matrix with @f$a_{ij} \leq 0@f$ for @f$j \neq i@f$, respectively. With
1499@f$a_{ij} \leq 0@f$, we can formulate bounds for the global time step @f$k@f$ as
1500follows:
1501@f{align*}{
1502 k_{\text{max}} &= \frac{ 1 }{ 1 - \theta }
1503 \min\left( \frac{ m_{ii} }{ a_{ii} } \right),~ \forall i,
1504 &
1505 k_{\text{min}} &= \frac{ 1 }{ \theta }
1506 \max\left( \frac{ m_{ij} }{ \left|a_{ij}\right| } \right),~ j \neq i.
1507@f}
1508In other words, the time step is constrained by <i>both a lower
1509and upper bound</i> in case of a Crank-Nicolson scheme. These bounds should be
1510considered along with the CFL condition to ensure significance of the performed
1511simulations.
1512
1513Being unable to make the time step as small as we want to get more
1514accuracy without losing the positivity property is annoying. It raises
1515the question of whether we can at least <i>compute</i> the minimal time step
1516we can choose to ensure positivity preservation in this particular tutorial.
1517Indeed, we can use
1518the SparseMatrix objects for both mass and stiffness that are created via
1519the MatrixCreator functions. Iterating through each entry via SparseMatrixIterators
1520lets us check for diagonal and off-diagonal entries to set a proper time step
1521dynamically. For quadratic matrices, the diagonal element is stored as the
1522first member of a row (see SparseMatrix documentation). An exemplary code
1523snippet on how to grab the entries of interest from the <code>mass_matrix</code>
1524is shown below.
1525
1526@code
1527Assert (mass_matrix.m() == mass_matrix.n(), ExcNotQuadratic());
1528const unsigned int num_rows = mass_matrix.m();
1529double mass_matrix_min_diag = std::numeric_limits<double>::max(),
1530 mass_matrix_max_offdiag = 0.;
1531
1532SparseMatrixIterators::Iterator<double,true> row_it (&mass_matrix, 0);
1533
1534for(unsigned int m = 0; m<num_rows; ++m)
1535{
1536 // check the diagonal element
1537 row_it = mass_matrix.begin(m);
1538 mass_matrix_min_diag = std::min(row_it->value(), mass_matrix_min_diag);
1539 ++row_it;
1540
1541 // check the off-diagonal elements
1542 for(; row_it != mass_matrix.end(m); ++row_it)
1543 mass_matrix_max_offdiag = std::max(row_it->value(), mass_matrix_max_offdiag);
1544}
1545@endcode
1546
1547Using the information so computed, we can bound the time step via the formulas
1548above.
1549 *
1550 *
1551<a name="step_26-PlainProg"></a>
1552<h1> The plain program</h1>
1553@include "step-26.cc"
1554*/
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void random(DoFHandler< dim, spacedim > &dof_handler)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
double diameter(const Triangulation< dim, spacedim > &tria)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)
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={})
void create_right_hand_side(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, const Function< spacedim, typename VectorType::value_type > &rhs, VectorType &rhs_vector, const AffineConstraints< typename VectorType::value_type > &constraints=AffineConstraints< typename VectorType::value_type >())
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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)