step-40.h

Go to the documentation of this file.

 
 *   #endif
 *       LA::MPI::PreconditionAMG preconditioner;
 *       preconditioner.initialize(system_matrix, data);
 *   
 *       solver.solve(system_matrix,
 *                    completely_distributed_solution,
 *                    system_rhs,
 *                    preconditioner);
 *   
 *       pcout << "   Solved in " << solver_control.last_step() << " iterations."
 *             << std::endl;
 *   
 *       constraints.distribute(completely_distributed_solution);
 *   
 *       locally_relevant_solution = completely_distributed_solution;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemrefine_grid"></a> 
 * <h4>LaplaceProblem::refine_grid</h4>
 * 
 
 * 
 * The function that estimates the error and refines the grid is again
 * almost exactly like the one in @ref step_6 "step-6". The only difference is that the
 * function that flags cells to be refined is now in namespace
 * parallel::distributed::GridRefinement -- a namespace that has functions
 * that can communicate between all involved processors and determine global
 * thresholds to use in deciding which cells to refine and which to coarsen.
 *   
 
 * 
 * Note that we didn't have to do anything special about the
 * KellyErrorEstimator class: we just give it a vector with as many elements
 * as the local triangulation has cells (locally owned cells, ghost cells,
 * and artificial ones), but it only fills those entries that correspond to
 * cells that are locally owned.
 * 
 * @code
 *     template <int dim>
 *     void LaplaceProblem<dim>::refine_grid()
 *     {
 *       TimerOutput::Scope t(computing_timer, "refine");
 *   
 *       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *       KellyErrorEstimator<dim>::estimate(
 *         dof_handler,
 *         QGauss<dim - 1>(fe.degree + 1),
 *         std::map<types::boundary_id, const Function<dim> *>(),
 *         locally_relevant_solution,
 *         estimated_error_per_cell);
 *       parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
 *         triangulation, estimated_error_per_cell, 0.3, 0.03);
 *       triangulation.execute_coarsening_and_refinement();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemoutput_results"></a> 
 * <h4>LaplaceProblem::output_results</h4>
 * 
 
 * 
 * Compared to the corresponding function in @ref step_6 "step-6", the one here is
 * a tad more complicated. There are two reasons: the first one is
 * that we do not just want to output the solution but also for each
 * cell which processor owns it (i.e. which "subdomain" it is
 * in). Secondly, as discussed at length in @ref step_17 "step-17" and @ref step_18 "step-18",
 * generating graphical data can be a bottleneck in
 * parallelizing. In those two programs, we simply generate one
 * output file per process. That worked because the
 * parallel::shared::Triangulation cannot be used with large numbers
 * of MPI processes anyway.  But this doesn't scale: Creating a
 * single file per processor will overwhelm the filesystem with a
 * large number of processors.
 *   
 
 * 
 * We here follow a more sophisticated approach that uses
 * high-performance, parallel IO routines using MPI I/O to write to
 * a small, fixed number of visualization files (here 8). We also
 * generate a .pvtu record referencing these .vtu files, which can
 * be opened directly in visualizatin tools like ParaView and VisIt.
 *   
 
 * 
 * To start, the top of the function looks like it usually does. In addition
 * to attaching the solution vector (the one that has entries for all locally
 * relevant, not only the locally owned, elements), we attach a data vector
 * that stores, for each cell, the subdomain the cell belongs to. This is
 * slightly tricky, because of course not every processor knows about every
 * cell. The vector we attach therefore has an entry for every cell that the
 * current processor has in its mesh (locally owned ones, ghost cells, and
 * artificial cells), but the DataOut class will ignore all entries that
 * correspond to cells that are not owned by the current processor. As a
 * consequence, it doesn't actually matter what values we write into these
 * vector entries: we simply fill the entire vector with the number of the
 * current MPI process (i.e. the subdomain_id of the current process); this
 * correctly sets the values we care for, i.e. the entries that correspond
 * to locally owned cells, while providing the wrong value for all other
 * elements -- but these are then ignored anyway.
 * 
 * @code
 *     template <int dim>
 *     void LaplaceProblem<dim>::output_results(const unsigned int cycle) const
 *     {
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(locally_relevant_solution, "u");
 *   
 *       Vector<float> subdomain(triangulation.n_active_cells());
 *       for (unsigned int i = 0; i < subdomain.size(); ++i)
 *         subdomain(i) = triangulation.locally_owned_subdomain();
 *       data_out.add_data_vector(subdomain, "subdomain");
 *   
 *       data_out.build_patches();
 *   
 * @endcode
 * 
 * The final step is to write this data to disk. We write up to 8 VTU files
 * in parallel with the help of MPI-IO. Additionally a PVTU record is
 * generated, which groups the written VTU files.
 * 
 * @code
 *       data_out.write_vtu_with_pvtu_record(
 *         "./", "solution", cycle, mpi_communicator, 2, 8);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="LaplaceProblemrun"></a> 
 * <h4>LaplaceProblem::run</h4>
 * 
 
 * 
 * The function that controls the overall behavior of the program is again
 * like the one in @ref step_6 "step-6". The minor difference are the use of
 * <code>pcout</code> instead of <code>std::cout</code> for output to the
 * console (see also @ref step_17 "step-17").
 *   
 
 * 
 * A functional difference to @ref step_6 "step-6" is the use of a square domain and that
 * we start with a slightly finer mesh (5 global refinement cycles) -- there
 * just isn't much of a point showing a massively %parallel program starting
 * on 4 cells (although admittedly the point is only slightly stronger
 * starting on 1024).
 * 
 * @code
 *     template <int dim>
 *     void LaplaceProblem<dim>::run()
 *     {
 *       pcout << "Running with "
 *   #ifdef USE_PETSC_LA
 *             << "PETSc"
 *   #else
 *             << "Trilinos"
 *   #endif
 *             << " on " << Utilities::MPI::n_mpi_processes(mpi_communicator)
 *             << " MPI rank(s)..." << std::endl;
 *   
 *       const unsigned int n_cycles = 8;
 *       for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
 *         {
 *           pcout << "Cycle " << cycle << ':' << std::endl;
 *   
 *           if (cycle == 0)
 *             {
 *               GridGenerator::hyper_cube(triangulation);
 *               triangulation.refine_global(5);
 *             }
 *           else
 *             refine_grid();
 *   
 *           setup_system();
 *   
 *           pcout << "   Number of active cells:       "
 *                 << triangulation.n_global_active_cells() << std::endl
 *                 << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl;
 *   
 *           assemble_system();
 *           solve();
 *   
 *           {
 *             TimerOutput::Scope t(computing_timer, "output");
 *             output_results(cycle);
 *           }
 *   
 *           computing_timer.print_summary();
 *           computing_timer.reset();
 *   
 *           pcout << std::endl;
 *         }
 *     }
 *   } // namespace Step40
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="main"></a> 
 * <h4>main()</h4>
 * 
 
 * 
 * The final function, <code>main()</code>, again has the same structure as in
 * all other programs, in particular @ref step_6 "step-6". Like the other programs that use
 * MPI, we have to initialize and finalize MPI, which is done using the helper
 * object Utilities::MPI::MPI_InitFinalize. The constructor of that class also
 * initializes libraries that depend on MPI, such as p4est, PETSc, SLEPc, and
 * Zoltan (though the last two are not used in this tutorial). The order here
 * is important: we cannot use any of these libraries until they are
 * initialized, so it does not make sense to do anything before creating an
 * instance of Utilities::MPI::MPI_InitFinalize.
 * 
 
 * 
 * After the solver finishes, the LaplaceProblem destructor will run followed
 * by Utilities::MPI::MPI_InitFinalize::~MPI_InitFinalize(). This order is
 * also important: Utilities::MPI::MPI_InitFinalize::~MPI_InitFinalize() calls
 * <code>PetscFinalize</code> (and finalization functions for other
 * libraries), which will delete any in-use PETSc objects. This must be done
 * after we destruct the Laplace solver to avoid double deletion
 * errors. Fortunately, due to the order of destructor call rules of C++, we
 * do not need to worry about any of this: everything happens in the correct
 * order (i.e., the reverse of the order of construction). The last function
 * called by Utilities::MPI::MPI_InitFinalize::~MPI_InitFinalize() is
 * <code>MPI_Finalize</code>: i.e., once this object is destructed the program
 * should exit since MPI will no longer be available.
 * 
 * @code
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         using namespace dealii;
 *         using namespace Step40;
 *   
 *         Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 *   
 *         LaplaceProblem<2> laplace_problem_2d;
 *         laplace_problem_2d.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
When you run the program, on a single processor or with your local MPI
installation on a few, you should get output like this:
@code
Cycle 0:
   Number of active cells:       1024
   Number of degrees of freedom: 4225
   Solved in 10 iterations.
 
 
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |     0.176s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| assembly                        |         1 |    0.0209s |        12% |
| output                          |         1 |    0.0189s |        11% |
| setup                           |         1 |    0.0299s |        17% |
| solve                           |         1 |    0.0419s |        24% |
+---------------------------------+-----------+------------+------------+
 
 
Cycle 1:
   Number of active cells:       1954
   Number of degrees of freedom: 8399
   Solved in 10 iterations.
 
 
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |     0.327s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| assembly                        |         1 |    0.0368s |        11% |
| output                          |         1 |    0.0208s |       6.4% |
| refine                          |         1 |     0.157s |        48% |
| setup                           |         1 |    0.0452s |        14% |
| solve                           |         1 |    0.0668s |        20% |
+---------------------------------+-----------+------------+------------+
 
 
Cycle 2:
   Number of active cells:       3664
   Number of degrees of freedom: 16183
   Solved in 11 iterations.
 
...
@endcode
 
The exact numbers differ, depending on how many processors we use;
this is due to the fact that the preconditioner depends on the
partitioning of the problem, the solution then differs in the last few
digits, and consequently the mesh refinement differs slightly.
The primary thing to notice here, though, is that the number of
iterations does not increase with the size of the problem. This
guarantees that we can efficiently solve even the largest problems.
 
When run on a sufficiently large number of machines (say a few
thousand), this program can relatively easily solve problems with well
over one billion unknowns in less than a minute. On the other hand,
such big problems can no longer be visualized, so we also ran the
program on only 16 processors. Here are a mesh, along with its
partitioning onto the 16 processors, and the corresponding solution:
 
<table width="100%">
<tr>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.mesh.png" alt="">
</td>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.solution.png" alt="">
</td>
</tr>
</table>
 
The mesh on the left has a mere 7,069 cells. This is of course a
problem we would easily have been able to solve already on a single
processor using @ref step_6 "step-6", but the point of the program was to show how
to write a program that scales to many more machines. For example,
here are two graphs that show how the run time of a large number of parts
of the program scales on problems with around 52 and 375 million degrees of
freedom if we take more and more processors (these and the next couple of
graphs are taken from an earlier version of the
@ref distributed_paper "Distributed Computing paper"; updated graphs showing
data of runs on even larger numbers of processors, and a lot
more interpretation can be found in the final version of the paper):
 
<table width="100%">
<tr>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.strong2.png" alt="">
</td>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.strong.png" alt="">
</td>
</tr>
</table>
 
As can clearly be seen, the program scales nicely to very large
numbers of processors.
(For a discussion of what we consider "scalable" programs, see
@ref GlossParallelScaling "this glossary entry".)
The curves, in particular the linear solver, become a
bit wobbly at the right end of the graphs since each processor has too little
to do to offset the cost of communication (the part of the whole problem each
processor has to solve in the above two examples is only 13,000 and 90,000
degrees of freedom when 4,096 processors are used; a good rule of thumb is that
parallel programs work well if each processor has at least 100,000 unknowns).
 
While the strong scaling graphs above show that we can solve a problem of
fixed size faster and faster if we take more and more processors, the more
interesting question may be how big problems can become so that they can still
be solved within a reasonable time on a machine of a particular size. We show
this in the following two graphs for 256 and 4096 processors:
 
<table width="100%">
<tr>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.256.png" alt="">
</td>
<td>
  <img src="https://www.dealii.org/images/steps/developer/step-40.4096.png" alt="">
</td>
</tr>
</table>
 
What these graphs show is that all parts of the program scale linearly with
the number of degrees of freedom. This time, lines are wobbly at the left as
the size of local problems is too small. For more discussions of these results
we refer to the @ref distributed_paper "Distributed Computing paper".
 
So how large are the largest problems one can solve? At the time of writing
this problem, the
limiting factor is that the program uses the BoomerAMG algebraic
multigrid method from the <a
href="http://acts.nersc.gov/hypre/" target="_top">Hypre package</a> as
a preconditioner, which unfortunately uses signed 32-bit integers to
index the elements of a %distributed matrix. This limits the size of
problems to @f$2^{31}-1=2,147,483,647@f$ degrees of freedom. From the graphs
above it is obvious that the scalability would extend beyond this
number, and one could expect that given more than the 4,096 machines
shown above would also further reduce the compute time. That said, one
can certainly expect that this limit will eventually be lifted by the
hypre developers.
 
On the other hand, this does not mean that deal.II cannot solve bigger
problems. Indeed, @ref step_37 "step-37" shows how one can solve problems that are not
just a little, but very substantially larger than anything we have shown
here.
 
 
 
<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
In a sense, this program is the ultimate solver for the Laplace
equation: it can essentially solve the equation to whatever accuracy
you want, if only you have enough processors available. Since the
Laplace equation by itself is not terribly interesting at this level
of accuracy, the more interesting possibilities for extension
therefore concern not so much this program but what comes beyond
it. For example, several of the other programs in this tutorial have
significant run times, especially in 3d. It would therefore be
interesting to use the techniques explained here to extend other
programs to support parallel distributed computations. We have done
this for @ref step_31 "step-31" in the @ref step_32 "step-32" tutorial program, but the same would
apply to, for example, @ref step_23 "step-23" and @ref step_25 "step-25" for hyperbolic time
dependent problems, @ref step_33 "step-33" for gas dynamics, or @ref step_35 "step-35" for the
Navier-Stokes equations.
 
Maybe equally interesting is the problem of postprocessing. As
mentioned above, we only show pictures of the solution and the mesh
for 16 processors because 4,096 processors solving 1 billion unknowns
would produce graphical output on the order of several 10
gigabyte. Currently, no program is able to visualize this amount of
data in any reasonable way unless it also runs on at least several
hundred processors. There are, however, approaches where visualization
programs directly communicate with solvers on each processor with each
visualization process rendering the part of the scene computed by the
solver on this processor. Implementing such an interface would allow
to quickly visualize things that are otherwise not amenable to
graphical display.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-40.cc"
*/