553 *
const double m = 0.5;
554 *
const double c1 = 0.;
555 *
const double c2 = 0.;
556 *
const double factor =
558 *
double result = 1.;
559 *
for (
unsigned int d = 0;
d < dim; ++
d)
570 * <a name=
"SineGordonProblemclass"></a>
571 * <h3>SineGordonProblem
class</h3>
575 * This is the main
class that builds on the class in @ref step_25 "step-25". However, we
577 * the geometry data. Also, we use a distributed
triangulation in
this
582 *
class SineGordonProblem
585 * SineGordonProblem();
591 *
void make_grid_and_dofs();
592 *
void output_results(
const unsigned int timestep_number);
594 * #ifdef DEAL_II_WITH_P4EST
609 *
const unsigned int n_global_refinements;
610 *
double time, time_step;
611 *
const double final_time;
612 *
const double cfl_number;
613 *
const unsigned int output_timestep_skip;
620 * <a name=
"SineGordonProblemSineGordonProblem"></a>
621 * <h4>SineGordonProblem::SineGordonProblem</h4>
625 * This is the constructor of the SineGordonProblem
class. The time interval
626 * and time step size are defined here. Moreover, we use the degree of the
627 * finite element that we defined at the top of the program to initialize a
628 *
FE_Q finite element based on Gauss-Lobatto support points. These points
629 * are convenient because in conjunction with a
QGaussLobatto quadrature
631 * compromising accuracy too much (note that the integration is inexact,
632 * though), see also the discussion in the introduction. Note that
FE_Q
633 * selects the Gauss-Lobatto nodal points by
default due to their improved
634 * conditioning versus equidistant points. To make things more
explicit, we
635 * state the selection of the nodal points nonetheless.
639 * SineGordonProblem<dim>::SineGordonProblem()
642 * #ifdef DEAL_II_WITH_P4EST
648 * , n_global_refinements(10 - 2 * dim)
652 * , cfl_number(.1 / fe_degree)
653 * , output_timestep_skip(200)
659 * <a name=
"SineGordonProblemmake_grid_and_dofs"></a>
660 * <h4>SineGordonProblem::make_grid_and_dofs</h4>
664 * As in @ref step_25
"step-25" this functions sets up a cube grid in <code>dim</code>
665 * dimensions of extent @f$[-15,15]@f$. We
refine the mesh more in the
center of
666 * the domain since the solution is concentrated there. We
first refine all
667 * cells whose
center is within a radius of 11, and then
refine once more
668 *
for a radius 6. This simple ad hoc refinement could be done better by
669 * adapting the mesh to the solution
using error estimators during the time
670 * stepping as done in other example programs, and
using
676 *
void SineGordonProblem<dim>::make_grid_and_dofs()
684 *
for (; cell != end_cell; ++cell)
685 *
if (cell->is_locally_owned())
686 *
if (cell->center().norm() < 11)
687 * cell->set_refine_flag();
692 *
for (; cell != end_cell; ++cell)
693 *
if (cell->is_locally_owned())
694 *
if (cell->center().norm() < 6)
695 * cell->set_refine_flag();
699 * pcout <<
" Number of global active cells: "
700 * #ifdef DEAL_II_WITH_P4EST
707 * dof_handler.distribute_dofs(fe);
709 * pcout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
715 * We generate hanging node constraints
for ensuring continuity of the
716 * solution. As in @ref step_40
"step-40", we need to equip the constraint
matrix with
717 * the
IndexSet of locally relevant degrees of freedom to avoid it to
718 * consume too much memory
for big problems. Next, the <code>
MatrixFree
719 * </code>
object for the problem is
set up. Note that we specify a
720 * particular scheme
for shared-memory parallelization (hence
one would
721 * use multithreading
for intra-node parallelism and not MPI; we here
722 * choose the standard option —
if we wanted to disable shared
723 * memory parallelization even in
case where there is more than
one TBB
724 * thread available in the program, we would choose
726 * instead of
using the
default QGauss quadrature argument, we supply a
728 * behavior. Finally, three solution vectors are initialized.
MatrixFree
729 * expects a particular layout of ghost indices (as it handles index
730 * access in MPI-local
numbers that need to match between the vector and
731 *
MatrixFree), so we just ask it to initialize the vectors to be sure the
732 * ghost exchange is properly handled.
736 * constraints.clear();
737 * constraints.reinit(locally_relevant_dofs);
739 * constraints.close();
745 * matrix_free_data.
reinit(dof_handler,
750 * matrix_free_data.initialize_dof_vector(solution);
751 * old_solution.reinit(solution);
752 * old_old_solution.reinit(solution);
760 * <a name=
"SineGordonProblemoutput_results"></a>
761 * <h4>SineGordonProblem::output_results</h4>
765 * This
function prints the
norm of the solution and writes the solution
766 * vector to a file. The
norm is standard (except
for the fact that we need
767 * to accumulate the norms over all processors
for the
parallel grid which
769 *
second is similar to what we did in @ref step_40
"step-40" or @ref step_37
"step-37". Note that we can
770 * use the same vector
for output as the
one used during computations: The
772 * all locally owned cells (
this is what is needed in the local evaluations,
773 * too), including ghost vector entries on these cells. This is the only
775 * as well as in
DataOut. The only action to take at this
point is to make
776 * sure that the vector updates its ghost values before we read from
777 * them. This is a feature present only in the
779 * and Trilinos, on the other hand, need to be copied to special vectors
780 * including ghost values (see the relevant section in @ref step_40 "step-40"). If we also
781 * wanted to access all degrees of freedom on ghost cells (
e.g. when
782 * computing error estimators that use the jump of solution over cell
783 * boundaries), we would need more information and create a vector
784 * initialized with locally relevant dofs just as in @ref step_40 "step-40". Observe also
785 * that we need to distribute constraints for output - they are not filled
786 * during computations (rather, they are interpolated on the fly in the
792 * SineGordonProblem<dim>::output_results(const
unsigned int timestep_number)
794 * constraints.distribute(solution);
797 * solution.update_ghost_values();
804 *
const double solution_norm =
809 * pcout <<
" Time:" << std::setw(8) << std::setprecision(3) << time
810 * <<
", solution norm: " << std::setprecision(5) << std::setw(7)
811 * << solution_norm << std::endl;
820 *
"./",
"solution", timestep_number, MPI_COMM_WORLD, 3);
827 * <a name=
"SineGordonProblemrun"></a>
832 * This
function is called by the main
function and steps into the
833 * subroutines of the
class.
837 * After printing some information about the
parallel setup, the
first
838 * action is to
set up the grid and the cell
operator. Then, the time step
839 * is computed from the CFL number given in the constructor and the finest
840 * mesh size. The finest mesh size is computed as the
diameter of the last
842 * the mesh. This is only possible
for meshes where all elements on a
level
843 * have the same size, otherwise,
one needs to
loop over all cells. Note
844 * that we need to query all the processors
for their finest cell since
845 * not all processors might hold a region where the mesh is at the finest
846 *
level. Then, we readjust the time step a little to hit the
final time
854 * pcout <<
"Number of MPI ranks: "
856 * pcout <<
"Number of threads on each rank: "
859 *
const unsigned int n_vect_bits = 8 *
sizeof(
double) * n_vect_doubles;
860 * pcout <<
"Vectorization over " << n_vect_doubles
861 * <<
" doubles = " << n_vect_bits <<
" bits ("
866 * make_grid_and_dofs();
868 *
const double local_min_cell_diameter =
870 *
const double global_min_cell_diameter =
872 * time_step = cfl_number * global_min_cell_diameter;
873 * time_step = (final_time - time) / (
int((final_time - time) / time_step));
874 * pcout <<
" Time step size: " << time_step
875 * <<
", finest cell: " << global_min_cell_diameter << std::endl
880 * Next the
initial value is
set. Since we have a two-step time stepping
881 * method, we also need a
value of the solution at time-time_step. For
882 * accurate results,
one would need to compute
this from the time
883 * derivative of the solution at
initial time, but here we ignore
this
889 * We then go on by writing the
initial state to file and collecting
890 * the two starting solutions in a <tt>std::vector</tt> of pointers that
892 * an instance of the <code> SineGordonOperation class </code> based on
893 * the finite element degree specified at the top of this file is
set up.
897 * InitialCondition<dim>(1, time),
900 * InitialCondition<dim>(1, time - time_step),
905 * previous_solutions({&old_solution, &old_old_solution});
907 * SineGordonOperation<dim, fe_degree> sine_gordon_op(matrix_free_data,
912 * Now
loop over the time steps. In each iteration, we
shift the solution
914 * `SineGordonOperator`
class. Then, we write the solution to a file. We
915 * clock the wall times
for the computational time needed as wall as the
916 * time needed to create the output and report the
numbers when the time
917 * stepping is finished.
921 * Note how
this shift is implemented: We simply
call the
swap method on
922 * the two vectors which swaps only some pointers without the need to
copy
923 * data around, a relatively expensive operation within an
explicit time
924 * stepping method. Let us see what happens in more detail: First, we
925 * exchange <code>old_solution</code> with <code>old_old_solution</code>,
926 * which means that <code>old_old_solution</code> gets
927 * <code>old_solution</code>, which is what we expect. Similarly,
928 * <code>old_solution</code> gets the content from <code>solution</code>
929 * in the next step. After
this, <code>solution</code> holds
930 * <code>old_old_solution</code>, but that will be overwritten during
this
934 *
unsigned int timestep_number = 1;
938 *
double output_time = 0;
939 *
for (time += time_step; time <= final_time;
940 * time += time_step, ++timestep_number)
943 * old_old_solution.swap(old_solution);
944 * old_solution.swap(solution);
945 * sine_gordon_op.apply(solution, previous_solutions);
949 *
if (timestep_number % output_timestep_skip == 0)
950 * output_results(timestep_number / output_timestep_skip);
955 * output_results(timestep_number / output_timestep_skip + 1);
959 * <<
" Performed " << timestep_number <<
" time steps." << std::endl;
961 * pcout <<
" Average wallclock time per time step: "
962 * << wtime / timestep_number <<
"s" << std::endl;
964 * pcout <<
" Spent " << output_time <<
"s on output and " << wtime
965 * <<
"s on computations." << std::endl;
974 * <a name=
"Thecodemaincodefunction"></a>
975 * <h3>The <code>main</code>
function</h3>
979 * As in @ref step_40
"step-40", we initialize MPI at the start of the program. Since we will
980 * in
general mix MPI parallelization with threads, we also
set the third
981 * argument in MPI_InitFinalize that controls the number of threads to an
982 *
invalid number, which means that the TBB library chooses the number of
983 * threads automatically, typically to the number of available cores in the
984 * system. As an alternative, you can also
set this number manually
if you
985 * want to
set a specific number of threads (
e.g. when MPI-only is required).
988 *
int main(
int argc,
char **argv)
990 *
using namespace Step48;
998 * SineGordonProblem<dimension> sg_problem;
1001 *
catch (std::exception &exc)
1003 * std::cerr << std::endl
1005 * <<
"----------------------------------------------------"
1007 * std::cerr <<
"Exception on processing: " << std::endl
1008 * << exc.what() << std::endl
1009 * <<
"Aborting!" << std::endl
1010 * <<
"----------------------------------------------------"
1017 * std::cerr << std::endl
1019 * <<
"----------------------------------------------------"
1021 * std::cerr <<
"Unknown exception!" << std::endl
1022 * <<
"Aborting!" << std::endl
1023 * <<
"----------------------------------------------------"
1031 <a name=
"Results"></a><h1>Results</h1>
1034 <a name=
"Comparisonwithasparsematrix"></a><h3>Comparison with a sparse
matrix</h3>
1037 In order to demonstrate the gain in
using the
MatrixFree class instead of
1038 the standard <code>deal.II</code> assembly routines for evaluating the
1039 information from old time steps, we study a simple serial
run of the code on a
1040 nonadaptive mesh. Since much time is spent on evaluating the sine function, we
1041 do not only show the
numbers of the full sine-Gordon equation but also for the
1042 wave equation (the sine-term skipped from the sine-Gordon equation). We use
1043 both
second and fourth order elements. The results are summarized in the
1046 <table align=
"center" class=
"doxtable">
1049 <th colspan=
"3">wave equation</th>
1050 <th colspan=
"2">sine-Gordon</th>
1061 <td>2D, @f$\mathcal{Q}_2@f$</td>
1062 <td align=
"right"> 0.0106</td>
1063 <td align=
"right"> 0.00971</td>
1064 <td align=
"right"> 0.109</td>
1065 <td align=
"right"> 0.0243</td>
1066 <td align=
"right"> 0.124</td>
1069 <td>2D, @f$\mathcal{Q}_4@f$</td>
1070 <td align=
"right"> 0.0328</td>
1071 <td align=
"right"> 0.0706</td>
1072 <td align=
"right"> 0.528</td>
1073 <td align=
"right"> 0.0714</td>
1074 <td align=
"right"> 0.502</td>
1077 <td>3D, @f$\mathcal{Q}_2@f$</td>
1078 <td align=
"right"> 0.0151</td>
1079 <td align=
"right"> 0.0320</td>
1080 <td align=
"right"> 0.331</td>
1081 <td align=
"right"> 0.0376</td>
1082 <td align=
"right"> 0.364</td>
1085 <td>3D, @f$\mathcal{Q}_4@f$</td>
1086 <td align=
"right"> 0.0918</td>
1087 <td align=
"right"> 0.844</td>
1088 <td align=
"right"> 6.83</td>
1089 <td align=
"right"> 0.194</td>
1090 <td align=
"right"> 6.95</td>
1094 It is apparent that the
matrix-
free code outperforms the standard assembly
1095 routines in deal.II by far. In 3D and
for fourth order elements,
one operator
1096 evaluation is also almost ten times as fast as a sparse
matrix-vector
1099 <a name=
"Parallelrunin2Dand3D"></a><h3>Parallel
run in 2D and 3D</h3>
1102 We start with the program output obtained on a workstation with 12 cores / 24
1103 threads (
one Intel Xeon E5-2687W v4 CPU running at 3.2 GHz, hyperthreading
1104 enabled), running the program in release mode:
1107 Number of MPI ranks: 1
1108 Number of threads on each rank: 24
1109 Vectorization over 4 doubles = 256 bits (AVX)
1111 Number of global active cells: 15412
1112 Number of degrees of freedom: 249065
1113 Time step size: 0.00292997, finest cell: 0.117188
1115 Time: -10, solution
norm: 9.5599
1116 Time: -9.41, solution
norm: 17.678
1117 Time: -8.83, solution
norm: 23.504
1118 Time: -8.24, solution
norm: 27.5
1119 Time: -7.66, solution
norm: 29.513
1120 Time: -7.07, solution
norm: 29.364
1121 Time: -6.48, solution
norm: 27.23
1122 Time: -5.9, solution
norm: 23.527
1123 Time: -5.31, solution
norm: 18.439
1124 Time: -4.73, solution
norm: 11.935
1125 Time: -4.14, solution
norm: 5.5284
1126 Time: -3.55, solution
norm: 8.0354
1127 Time: -2.97, solution
norm: 14.707
1128 Time: -2.38, solution
norm: 20
1129 Time: -1.8, solution
norm: 22.834
1130 Time: -1.21, solution
norm: 22.771
1131 Time: -0.624, solution
norm: 20.488
1132 Time: -0.0381, solution
norm: 16.697
1133 Time: 0.548, solution
norm: 11.221
1134 Time: 1.13, solution
norm: 5.3912
1135 Time: 1.72, solution
norm: 8.4528
1136 Time: 2.31, solution
norm: 14.335
1137 Time: 2.89, solution
norm: 18.555
1138 Time: 3.48, solution
norm: 20.894
1139 Time: 4.06, solution
norm: 21.305
1140 Time: 4.65, solution
norm: 19.903
1141 Time: 5.24, solution
norm: 16.864
1142 Time: 5.82, solution
norm: 12.223
1143 Time: 6.41, solution
norm: 6.758
1144 Time: 6.99, solution
norm: 7.2423
1145 Time: 7.58, solution
norm: 12.888
1146 Time: 8.17, solution
norm: 17.273
1147 Time: 8.75, solution
norm: 19.654
1148 Time: 9.34, solution
norm: 19.838
1149 Time: 9.92, solution
norm: 17.964
1150 Time: 10, solution
norm: 17.595
1152 Performed 6826 time steps.
1153 Average wallclock time per time step: 0.0013453s
1154 Spent 14.976s on output and 9.1831s on computations.
1157 In 3D, the respective output looks like
1160 Number of MPI ranks: 1
1161 Number of threads on each rank: 24
1162 Vectorization over 4 doubles = 256 bits (AVX)
1164 Number of global active cells: 17592
1165 Number of degrees of freedom: 1193881
1166 Time step size: 0.0117233, finest cell: 0.46875
1168 Time: -10, solution
norm: 29.558
1169 Time: -7.66, solution
norm: 129.13
1170 Time: -5.31, solution
norm: 67.753
1171 Time: -2.97, solution
norm: 79.245
1172 Time: -0.621, solution
norm: 123.52
1173 Time: 1.72, solution
norm: 43.525
1174 Time: 4.07, solution
norm: 93.285
1175 Time: 6.41, solution
norm: 97.722
1176 Time: 8.76, solution
norm: 36.734
1177 Time: 10, solution
norm: 94.115
1179 Performed 1706 time steps.
1180 Average wallclock time per time step: 0.0084542s
1181 Spent 16.766s on output and 14.423s on computations.
1184 It takes 0.008 seconds
for one time step with more than a million
1185 degrees of freedom (note that we would need many processors to reach such
1186 numbers when solving linear systems).
1188 If we replace the thread-parallelization by a pure MPI parallelization, the
1189 timings change into:
1191 $ mpirun -n 24 ./step-48
1192 Number of MPI ranks: 24
1193 Number of threads on each rank: 1
1194 Vectorization over 4 doubles = 256 bits (AVX)
1196 Performed 1706 time steps.
1197 Average wallclock time per time step: 0.0051747s
1198 Spent 2.0535s on output and 8.828s on computations.
1201 We observe a dramatic speedup
for the output (which makes sense, given that
1202 most code of the output is not parallelized via threads, whereas it is
for
1203 MPI), but less than the theoretical factor of 12 we would expect from the
1204 parallelism. More interestingly, the computations also get faster when
1205 switching from the threads-only variant to the MPI-only variant. This is a
1207 2019). The main reason is that the decisions regarding work on conflicting
1208 cell batches made to enable execution in
parallel are overly pessimistic:
1209 While they ensure that no work on neighboring cells is done on different
1210 threads at the same time,
this conservative setting implies that data from
1211 neighboring cells is also evicted from caches by the time neighbors get
1212 touched. Furthermore, the current scheme is not able to provide a constant
1213 load
for all 24 threads
for the given mesh with 17,592 cells.
1215 The current program allows to also mix MPI parallelization with thread
1216 parallelization. This is most beneficial when running programs on clusters
1217 with multiple nodes,
using MPI
for the inter-node parallelization and threads
1218 for the intra-node parallelization. On the workstation used above, we can
run
1219 threads in the hyperthreading region (i.e.,
using 2 threads
for each of the 12
1220 MPI ranks). An important setting
for mixing MPI with threads is to ensure
1221 proper binning of tasks to CPUs. On many clusters the placing is either
1222 automatically via the `mpirun/mpiexec` environment, or there can be manual
1223 settings. Here, we simply report the
run times the plain version of the
1224 program (noting that things could be improved towards the timings of the
1225 MPI-only program when proper pinning is done):
1227 $ mpirun -n 12 ./step-48
1228 Number of MPI ranks: 12
1229 Number of threads on each rank: 2
1230 Vectorization over 4 doubles = 256 bits (AVX)
1232 Performed 1706 time steps.
1233 Average wallclock time per time step: 0.0056651s
1234 Spent 2.5175s on output and 9.6646s on computations.
1239 <a name=
"Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1242 There are several things in this program that could be improved to make it
1243 even more efficient (besides improved boundary conditions and physical
1244 stuff as discussed in @ref step_25
"step-25"):
1246 <ul> <li> <
b>Faster evaluation of sine terms:</
b> As becomes obvious
1247 from the comparison of the plain wave equation and the sine-Gordon
1248 equation above, the evaluation of the sine terms dominates the total
1249 time for the finite element operator application. There are a few
1250 reasons for this: Firstly, the deal.II sine computation of a
1252 the operator application). This could be cured by handing the sine
1253 computation to a library with vectorized sine computations like
1254 Intel
's math kernel library (MKL). By using the function
1255 <code>vdSin</code> in MKL, the program uses half the computing time
1256 in 2D and 40 percent less time in 3D. On the other hand, the sine
1257 computation is structurally much more complicated than the simple
1258 arithmetic operations like additions and multiplications in the rest
1259 of the local operation.
1261 <li> <b>Higher order time stepping:</b> While the implementation allows for
1262 arbitrary order in the spatial part (by adjusting the degree of the finite
1263 element), the time stepping scheme is a standard second-order leap-frog
1264 scheme. Since solutions in wave propagation problems are usually very
1265 smooth, the error is likely dominated by the time stepping part. Of course,
1266 this could be cured by using smaller time steps (at a fixed spatial
1267 resolution), but it would be more efficient to use higher order time
1268 stepping as well. While it would be straight-forward to do so for a
1269 first-order system (use some Runge–Kutta scheme of higher order,
1270 probably combined with adaptive time step selection like the <a
1271 href="http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method">Dormand–Prince
1272 method</a>), it is more challenging for the second-order formulation. At
1273 least in the finite difference community, people usually use the PDE to find
1274 spatial correction terms that improve the temporal error.
1279 <a name="PlainProg"></a>
1280 <h1> The plain program</h1>
1281 @include "step-48.cc"