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