step-32.h

Go to the documentation of this file.

 
 *       constexpr double kappa                 = 1e-6;    /* m^2 / s    */
 *       constexpr double reference_density     = 3300;    /* kg / m^3   */
 *       constexpr double reference_temperature = 293;     /* K          */
 *       constexpr double expansion_coefficient = 2e-5;    /* 1/K        */
 *       constexpr double specific_heat         = 1250;    /* J / K / kg */
 *       constexpr double radiogenic_heating    = 7.4e-12; /* W / kg     */
 *   
 *   
 *       constexpr double R0 = 6371000. - 2890000.; /* m          */
 *       constexpr double R1 = 6371000. - 35000.;   /* m          */
 *   
 *       constexpr double T0 = 4000 + 273; /* K          */
 *       constexpr double T1 = 700 + 273;  /* K          */
 *   
 *   
 * @endcode
 * 
 * The next set of definitions are for functions that encode the density
 * as a function of temperature, the gravity vector, and the initial
 * values for the temperature. Again, all of these (along with the values
 * they compute) are discussed in the introduction:
 * 
 * @code
 *       double density(const double temperature)
 *       {
 *         return (
 *           reference_density *
 *           (1 - expansion_coefficient * (temperature - reference_temperature)));
 *       }
 *   
 *   
 *       template <int dim>
 *       Tensor<1, dim> gravity_vector(const Point<dim> &p)
 *       {
 *         const double r = p.norm();
 *         return -(1.245e-6 * r + 7.714e13 / r / r) * p / r;
 *       }
 *   
 *   
 *   
 *       template <int dim>
 *       class TemperatureInitialValues : public Function<dim>
 *       {
 *       public:
 *         TemperatureInitialValues()
 *           : Function<dim>(1)
 *         {}
 *   
 *         virtual double value(const Point<dim>  &p,
 *                              const unsigned int component = 0) const override;
 *   
 *         virtual void vector_value(const Point<dim> &p,
 *                                   Vector<double>   &value) const override;
 *       };
 *   
 *   
 *   
 *       template <int dim>
 *       double TemperatureInitialValues<dim>::value(const Point<dim> &p,
 *                                                   const unsigned int) const
 *       {
 *         const double r = p.norm();
 *         const double h = R1 - R0;
 *   
 *         const double s = (r - R0) / h;
 *         const double q =
 *           (dim == 3) ? std::max(0.0, cos(numbers::PI * abs(p[2] / R1))) : 1.0;
 *         const double phi = std::atan2(p[0], p[1]);
 *         const double tau = s + 0.2 * s * (1 - s) * std::sin(6 * phi) * q;
 *   
 *         return T0 * (1.0 - tau) + T1 * tau;
 *       }
 *   
 *   
 *       template <int dim>
 *       void
 *       TemperatureInitialValues<dim>::vector_value(const Point<dim> &p,
 *                                                   Vector<double>   &values) const
 *       {
 *         for (unsigned int c = 0; c < this->n_components; ++c)
 *           values(c) = TemperatureInitialValues<dim>::value(p, c);
 *       }
 *   
 *   
 * @endcode
 * 
 * As mentioned in the introduction we need to rescale the pressure to
 * avoid the relative ill-conditioning of the momentum and mass
 * conservation equations. The scaling factor is @f$\frac{\eta}{L}@f$ where
 * @f$L@f$ was a typical length scale. By experimenting it turns out that a
 * good length scale is the diameter of plumes, which is around 10 km:
 * 
 * @code
 *       constexpr double pressure_scaling = eta / 10000;
 *   
 * @endcode
 * 
 * The final number in this namespace is a constant that denotes the
 * number of seconds per (average, tropical) year. We use this only when
 * generating screen output: internally, all computations of this program
 * happen in SI units (kilogram, meter, seconds) but writing geological
 * times in seconds yields numbers that one can't relate to reality, and
 * so we convert to years using the factor defined here:
 * 
 * @code
 *       const double year_in_seconds = 60 * 60 * 24 * 365.2425;
 *   
 *     } // namespace EquationData
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-PreconditioningtheStokessystem"></a> 
 * <h3>Preconditioning the Stokes system</h3>
 * 
 
 * 
 * This namespace implements the preconditioner. As discussed in the
 * introduction, this preconditioner differs in a number of key portions
 * from the one used in @ref step_31 "step-31". Specifically, it is a right preconditioner,
 * implementing the matrix
 * @f{align*}{
 * \left(\begin{array}{cc}A^{-1} & -A^{-1}B^TS^{-1}
 * \\0 & S^{-1}
 * \end{array}\right)
 * @f}
 * where the two inverse matrix operations
 * are approximated by linear solvers or, if the right flag is given to the
 * constructor of this class, by a single AMG V-cycle for the velocity
 * block. The three code blocks of the <code>vmult</code> function implement
 * the multiplications with the three blocks of this preconditioner matrix
 * and should be self explanatory if you have read through @ref step_31 "step-31" or the
 * discussion of composing solvers in @ref step_20 "step-20".
 * 
 * @code
 *     namespace LinearSolvers
 *     {
 *       template <class PreconditionerTypeA, class PreconditionerTypeMp>
 *       class BlockSchurPreconditioner : public EnableObserverPointer
 *       {
 *       public:
 *         BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
 *                                  const TrilinosWrappers::BlockSparseMatrix &Spre,
 *                                  const PreconditionerTypeMp &Mppreconditioner,
 *                                  const PreconditionerTypeA  &Apreconditioner,
 *                                  const bool                  do_solve_A)
 *           : stokes_matrix(&S)
 *           , stokes_preconditioner_matrix(&Spre)
 *           , mp_preconditioner(Mppreconditioner)
 *           , a_preconditioner(Apreconditioner)
 *           , do_solve_A(do_solve_A)
 *         {}
 *   
 *         void vmult(TrilinosWrappers::MPI::BlockVector       &dst,
 *                    const TrilinosWrappers::MPI::BlockVector &src) const
 *         {
 *           TrilinosWrappers::MPI::Vector utmp(src.block(0));
 *   
 *           {
 *             SolverControl solver_control(5000, 1e-6 * src.block(1).l2_norm());
 *   
 *             SolverCG<TrilinosWrappers::MPI::Vector> solver(solver_control);
 *   
 *             solver.solve(stokes_preconditioner_matrix->block(1, 1),
 *                          dst.block(1),
 *                          src.block(1),
 *                          mp_preconditioner);
 *   
 *             dst.block(1) *= -1.0;
 *           }
 *   
 *           {
 *             stokes_matrix->block(0, 1).vmult(utmp, dst.block(1));
 *             utmp *= -1.0;
 *             utmp.add(src.block(0));
 *           }
 *   
 *           if (do_solve_A == true)
 *             {
 *               SolverControl solver_control(5000, utmp.l2_norm() * 1e-2);
 *               TrilinosWrappers::SolverCG solver(solver_control);
 *               solver.solve(stokes_matrix->block(0, 0),
 *                            dst.block(0),
 *                            utmp,
 *                            a_preconditioner);
 *             }
 *           else
 *             a_preconditioner.vmult(dst.block(0), utmp);
 *         }
 *   
 *       private:
 *         const ObserverPointer<const TrilinosWrappers::BlockSparseMatrix>
 *           stokes_matrix;
 *         const ObserverPointer<const TrilinosWrappers::BlockSparseMatrix>
 *                                     stokes_preconditioner_matrix;
 *         const PreconditionerTypeMp &mp_preconditioner;
 *         const PreconditionerTypeA  &a_preconditioner;
 *         const bool                  do_solve_A;
 *       };
 *     } // namespace LinearSolvers
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-Definitionofassemblydatastructures"></a> 
 * <h3>Definition of assembly data structures</h3>
 *   
 
 * 
 * As described in the introduction, we will use the WorkStream mechanism
 * discussed in the @ref threads topic to parallelize operations among the
 * processors of a single machine. The WorkStream class requires that data
 * is passed around in two kinds of data structures, one for scratch data
 * and one to pass data from the assembly function to the function that
 * copies local contributions into global objects.
 *   
 
 * 
 * The following namespace (and the two sub-namespaces) contains a
 * collection of data structures that serve this purpose, one pair for each
 * of the four operations discussed in the introduction that we will want to
 * parallelize. Each assembly routine gets two sets of data: a Scratch array
 * that collects all the classes and arrays that are used for the
 * calculation of the cell contribution, and a CopyData array that keeps
 * local matrices and vectors which will be written into the global
 * matrix. Whereas CopyData is a container for the final data that is
 * written into the global matrices and vector (and, thus, absolutely
 * necessary), the Scratch arrays are merely there for performance reasons
 * &mdash; it would be much more expensive to set up a FEValues object on
 * each cell, than creating it only once and updating some derivative data.
 *   
 
 * 
 * @ref step_31 "step-31" had four assembly routines: One for the preconditioner matrix of
 * the Stokes system, one for the Stokes matrix and right hand side, one for
 * the temperature matrices and one for the right hand side of the
 * temperature equation. We here organize the scratch arrays and CopyData
 * objects for each of those four assembly components using a
 * <code>struct</code> environment (since we consider these as temporary
 * objects we pass around, rather than classes that implement functionality
 * of their own, though this is a more subjective point of view to
 * distinguish between <code>struct</code>s and <code>class</code>es).
 *   
 
 * 
 * Regarding the Scratch objects, each struct is equipped with a constructor
 * that creates an @ref FEValues object using the @ref FiniteElement,
 * Quadrature, @ref Mapping (which describes the interpolation of curved
 * boundaries), and @ref UpdateFlags instances. Moreover, we manually
 * implement a copy constructor (since the FEValues class is not copyable by
 * itself), and provide some additional vector fields that are used to hold
 * intermediate data during the computation of local contributions.
 *   
 
 * 
 * Let us start with the scratch arrays and, specifically, the one used for
 * assembly of the Stokes preconditioner:
 * 
 * @code
 *     namespace Assembly
 *     {
 *       namespace Scratch
 *       {
 *         template <int dim>
 *         struct StokesPreconditioner
 *         {
 *           StokesPreconditioner(const FiniteElement<dim> &stokes_fe,
 *                                const Quadrature<dim>    &stokes_quadrature,
 *                                const Mapping<dim>       &mapping,
 *                                const UpdateFlags         update_flags);
 *   
 *           StokesPreconditioner(const StokesPreconditioner &data);
 *   
 *   
 *           FEValues<dim> stokes_fe_values;
 *   
 *           std::vector<Tensor<2, dim>> grad_phi_u;
 *           std::vector<double>         phi_p;
 *         };
 *   
 *         template <int dim>
 *         StokesPreconditioner<dim>::StokesPreconditioner(
 *           const FiniteElement<dim> &stokes_fe,
 *           const Quadrature<dim>    &stokes_quadrature,
 *           const Mapping<dim>       &mapping,
 *           const UpdateFlags         update_flags)
 *           : stokes_fe_values(mapping, stokes_fe, stokes_quadrature, update_flags)
 *           , grad_phi_u(stokes_fe.n_dofs_per_cell())
 *           , phi_p(stokes_fe.n_dofs_per_cell())
 *         {}
 *   
 *   
 *   
 *         template <int dim>
 *         StokesPreconditioner<dim>::StokesPreconditioner(
 *           const StokesPreconditioner &scratch)
 *           : stokes_fe_values(scratch.stokes_fe_values.get_mapping(),
 *                              scratch.stokes_fe_values.get_fe(),
 *                              scratch.stokes_fe_values.get_quadrature(),
 *                              scratch.stokes_fe_values.get_update_flags())
 *           , grad_phi_u(scratch.grad_phi_u)
 *           , phi_p(scratch.phi_p)
 *         {}
 *   
 *   
 *   
 * @endcode
 * 
 * The next one is the scratch object used for the assembly of the full
 * Stokes system. Observe that we derive the StokesSystem scratch class
 * from the StokesPreconditioner class above. We do this because all the
 * objects that are necessary for the assembly of the preconditioner are
 * also needed for the actual matrix system and right hand side, plus
 * some extra data. This makes the program more compact. Note also that
 * the assembly of the Stokes system and the temperature right hand side
 * further down requires data from temperature and velocity,
 * respectively, so we actually need two FEValues objects for those two
 * cases.
 * 
 * @code
 *         template <int dim>
 *         struct StokesSystem : public StokesPreconditioner<dim>
 *         {
 *           StokesSystem(const FiniteElement<dim> &stokes_fe,
 *                        const Mapping<dim>       &mapping,
 *                        const Quadrature<dim>    &stokes_quadrature,
 *                        const UpdateFlags         stokes_update_flags,
 *                        const FiniteElement<dim> &temperature_fe,
 *                        const UpdateFlags         temperature_update_flags);
 *   
 *           StokesSystem(const StokesSystem<dim> &data);
 *   
 *   
 *           FEValues<dim> temperature_fe_values;
 *   
 *           std::vector<Tensor<1, dim>>          phi_u;
 *           std::vector<SymmetricTensor<2, dim>> grads_phi_u;
 *           std::vector<double>                  div_phi_u;
 *   
 *           std::vector<double> old_temperature_values;
 *         };
 *   
 *   
 *         template <int dim>
 *         StokesSystem<dim>::StokesSystem(
 *           const FiniteElement<dim> &stokes_fe,
 *           const Mapping<dim>       &mapping,
 *           const Quadrature<dim>    &stokes_quadrature,
 *           const UpdateFlags         stokes_update_flags,
 *           const FiniteElement<dim> &temperature_fe,
 *           const UpdateFlags         temperature_update_flags)
 *           : StokesPreconditioner<dim>(stokes_fe,
 *                                       stokes_quadrature,
 *                                       mapping,
 *                                       stokes_update_flags)
 *           , temperature_fe_values(mapping,
 *                                   temperature_fe,
 *                                   stokes_quadrature,
 *                                   temperature_update_flags)
 *           , phi_u(stokes_fe.n_dofs_per_cell())
 *           , grads_phi_u(stokes_fe.n_dofs_per_cell())
 *           , div_phi_u(stokes_fe.n_dofs_per_cell())
 *           , old_temperature_values(stokes_quadrature.size())
 *         {}
 *   
 *   
 *         template <int dim>
 *         StokesSystem<dim>::StokesSystem(const StokesSystem<dim> &scratch)
 *           : StokesPreconditioner<dim>(scratch)
 *           , temperature_fe_values(
 *               scratch.temperature_fe_values.get_mapping(),
 *               scratch.temperature_fe_values.get_fe(),
 *               scratch.temperature_fe_values.get_quadrature(),
 *               scratch.temperature_fe_values.get_update_flags())
 *           , phi_u(scratch.phi_u)
 *           , grads_phi_u(scratch.grads_phi_u)
 *           , div_phi_u(scratch.div_phi_u)
 *           , old_temperature_values(scratch.old_temperature_values)
 *         {}
 *   
 *   
 * @endcode
 * 
 * After defining the objects used in the assembly of the Stokes system,
 * we do the same for the assembly of the matrices necessary for the
 * temperature system. The general structure is very similar:
 * 
 * @code
 *         template <int dim>
 *         struct TemperatureMatrix
 *         {
 *           TemperatureMatrix(const FiniteElement<dim> &temperature_fe,
 *                             const Mapping<dim>       &mapping,
 *                             const Quadrature<dim>    &temperature_quadrature);
 *   
 *           TemperatureMatrix(const TemperatureMatrix &data);
 *   
 *   
 *           FEValues<dim> temperature_fe_values;
 *   
 *           std::vector<double>         phi_T;
 *           std::vector<Tensor<1, dim>> grad_phi_T;
 *         };
 *   
 *   
 *         template <int dim>
 *         TemperatureMatrix<dim>::TemperatureMatrix(
 *           const FiniteElement<dim> &temperature_fe,
 *           const Mapping<dim>       &mapping,
 *           const Quadrature<dim>    &temperature_quadrature)
 *           : temperature_fe_values(mapping,
 *                                   temperature_fe,
 *                                   temperature_quadrature,
 *                                   update_values | update_gradients |
 *                                     update_JxW_values)
 *           , phi_T(temperature_fe.n_dofs_per_cell())
 *           , grad_phi_T(temperature_fe.n_dofs_per_cell())
 *         {}
 *   
 *   
 *         template <int dim>
 *         TemperatureMatrix<dim>::TemperatureMatrix(
 *           const TemperatureMatrix &scratch)
 *           : temperature_fe_values(
 *               scratch.temperature_fe_values.get_mapping(),
 *               scratch.temperature_fe_values.get_fe(),
 *               scratch.temperature_fe_values.get_quadrature(),
 *               scratch.temperature_fe_values.get_update_flags())
 *           , phi_T(scratch.phi_T)
 *           , grad_phi_T(scratch.grad_phi_T)
 *         {}
 *   
 *   
 * @endcode
 * 
 * The final scratch object is used in the assembly of the right hand
 * side of the temperature system. This object is significantly larger
 * than the ones above because a lot more quantities enter the
 * computation of the right hand side of the temperature equation. In
 * particular, the temperature values and gradients of the previous two
 * time steps need to be evaluated at the quadrature points, as well as
 * the velocities and the strain rates (i.e. the symmetric gradients of
 * the velocity) that enter the right hand side as friction heating
 * terms. Despite the number of terms, the following should be rather
 * self explanatory:
 * 
 * @code
 *         template <int dim>
 *         struct TemperatureRHS
 *         {
 *           TemperatureRHS(const FiniteElement<dim> &temperature_fe,
 *                          const FiniteElement<dim> &stokes_fe,
 *                          const Mapping<dim>       &mapping,
 *                          const Quadrature<dim>    &quadrature);
 *   
 *           TemperatureRHS(const TemperatureRHS &data);
 *   
 *   
 *           FEValues<dim> temperature_fe_values;
 *           FEValues<dim> stokes_fe_values;
 *   
 *           std::vector<double>         phi_T;
 *           std::vector<Tensor<1, dim>> grad_phi_T;
 *   
 *           std::vector<Tensor<1, dim>> old_velocity_values;
 *           std::vector<Tensor<1, dim>> old_old_velocity_values;
 *   
 *           std::vector<SymmetricTensor<2, dim>> old_strain_rates;
 *           std::vector<SymmetricTensor<2, dim>> old_old_strain_rates;
 *   
 *           std::vector<double>         old_temperature_values;
 *           std::vector<double>         old_old_temperature_values;
 *           std::vector<Tensor<1, dim>> old_temperature_grads;
 *           std::vector<Tensor<1, dim>> old_old_temperature_grads;
 *           std::vector<double>         old_temperature_laplacians;
 *           std::vector<double>         old_old_temperature_laplacians;
 *         };
 *   
 *   
 *         template <int dim>
 *         TemperatureRHS<dim>::TemperatureRHS(
 *           const FiniteElement<dim> &temperature_fe,
 *           const FiniteElement<dim> &stokes_fe,
 *           const Mapping<dim>       &mapping,
 *           const Quadrature<dim>    &quadrature)
 *           : temperature_fe_values(mapping,
 *                                   temperature_fe,
 *                                   quadrature,
 *                                   update_values | update_gradients |
 *                                     update_hessians | update_quadrature_points |
 *                                     update_JxW_values)
 *           , stokes_fe_values(mapping,
 *                              stokes_fe,
 *                              quadrature,
 *                              update_values | update_gradients)
 *           , phi_T(temperature_fe.n_dofs_per_cell())
 *           , grad_phi_T(temperature_fe.n_dofs_per_cell())
 *           ,
 *   
 *           old_velocity_values(quadrature.size())
 *           , old_old_velocity_values(quadrature.size())
 *           , old_strain_rates(quadrature.size())
 *           , old_old_strain_rates(quadrature.size())
 *           ,
 *   
 *           old_temperature_values(quadrature.size())
 *           , old_old_temperature_values(quadrature.size())
 *           , old_temperature_grads(quadrature.size())
 *           , old_old_temperature_grads(quadrature.size())
 *           , old_temperature_laplacians(quadrature.size())
 *           , old_old_temperature_laplacians(quadrature.size())
 *         {}
 *   
 *   
 *         template <int dim>
 *         TemperatureRHS<dim>::TemperatureRHS(const TemperatureRHS &scratch)
 *           : temperature_fe_values(
 *               scratch.temperature_fe_values.get_mapping(),
 *               scratch.temperature_fe_values.get_fe(),
 *               scratch.temperature_fe_values.get_quadrature(),
 *               scratch.temperature_fe_values.get_update_flags())
 *           , stokes_fe_values(scratch.stokes_fe_values.get_mapping(),
 *                              scratch.stokes_fe_values.get_fe(),
 *                              scratch.stokes_fe_values.get_quadrature(),
 *                              scratch.stokes_fe_values.get_update_flags())
 *           , phi_T(scratch.phi_T)
 *           , grad_phi_T(scratch.grad_phi_T)
 *           ,
 *   
 *           old_velocity_values(scratch.old_velocity_values)
 *           , old_old_velocity_values(scratch.old_old_velocity_values)
 *           , old_strain_rates(scratch.old_strain_rates)
 *           , old_old_strain_rates(scratch.old_old_strain_rates)
 *           ,
 *   
 *           old_temperature_values(scratch.old_temperature_values)
 *           , old_old_temperature_values(scratch.old_old_temperature_values)
 *           , old_temperature_grads(scratch.old_temperature_grads)
 *           , old_old_temperature_grads(scratch.old_old_temperature_grads)
 *           , old_temperature_laplacians(scratch.old_temperature_laplacians)
 *           , old_old_temperature_laplacians(scratch.old_old_temperature_laplacians)
 *         {}
 *       } // namespace Scratch
 *   
 *   
 * @endcode
 * 
 * The CopyData objects are even simpler than the Scratch objects as all
 * they have to do is to store the results of local computations until
 * they can be copied into the global matrix or vector objects. These
 * structures therefore only need to provide a constructor, a copy
 * operation, and some arrays for local matrix, local vectors and the
 * relation between local and global degrees of freedom (a.k.a.
 * <code>local_dof_indices</code>). Again, we have one such structure for
 * each of the four operations we will parallelize using the WorkStream
 * class:
 * 
 * @code
 *       namespace CopyData
 *       {
 *         template <int dim>
 *         struct StokesPreconditioner
 *         {
 *           StokesPreconditioner(const FiniteElement<dim> &stokes_fe);
 *           StokesPreconditioner(const StokesPreconditioner &data);
 *           StokesPreconditioner &operator=(const StokesPreconditioner &) = default;
 *   
 *           FullMatrix<double>                   local_matrix;
 *           std::vector<types::global_dof_index> local_dof_indices;
 *         };
 *   
 *         template <int dim>
 *         StokesPreconditioner<dim>::StokesPreconditioner(
 *           const FiniteElement<dim> &stokes_fe)
 *           : local_matrix(stokes_fe.n_dofs_per_cell(), stokes_fe.n_dofs_per_cell())
 *           , local_dof_indices(stokes_fe.n_dofs_per_cell())
 *         {}
 *   
 *         template <int dim>
 *         StokesPreconditioner<dim>::StokesPreconditioner(
 *           const StokesPreconditioner &data)
 *           : local_matrix(data.local_matrix)
 *           , local_dof_indices(data.local_dof_indices)
 *         {}
 *   
 *   
 *   
 *         template <int dim>
 *         struct StokesSystem : public StokesPreconditioner<dim>
 *         {
 *           StokesSystem(const FiniteElement<dim> &stokes_fe);
 *   
 *           Vector<double> local_rhs;
 *         };
 *   
 *         template <int dim>
 *         StokesSystem<dim>::StokesSystem(const FiniteElement<dim> &stokes_fe)
 *           : StokesPreconditioner<dim>(stokes_fe)
 *           , local_rhs(stokes_fe.n_dofs_per_cell())
 *         {}
 *   
 *   
 *   
 *         template <int dim>
 *         struct TemperatureMatrix
 *         {
 *           TemperatureMatrix(const FiniteElement<dim> &temperature_fe);
 *   
 *           FullMatrix<double>                   local_mass_matrix;
 *           FullMatrix<double>                   local_stiffness_matrix;
 *           std::vector<types::global_dof_index> local_dof_indices;
 *         };
 *   
 *         template <int dim>
 *         TemperatureMatrix<dim>::TemperatureMatrix(
 *           const FiniteElement<dim> &temperature_fe)
 *           : local_mass_matrix(temperature_fe.n_dofs_per_cell(),
 *                               temperature_fe.n_dofs_per_cell())
 *           , local_stiffness_matrix(temperature_fe.n_dofs_per_cell(),
 *                                    temperature_fe.n_dofs_per_cell())
 *           , local_dof_indices(temperature_fe.n_dofs_per_cell())
 *         {}
 *   
 *   
 *   
 *         template <int dim>
 *         struct TemperatureRHS
 *         {
 *           TemperatureRHS(const FiniteElement<dim> &temperature_fe);
 *   
 *           Vector<double>                       local_rhs;
 *           std::vector<types::global_dof_index> local_dof_indices;
 *           FullMatrix<double>                   matrix_for_bc;
 *         };
 *   
 *         template <int dim>
 *         TemperatureRHS<dim>::TemperatureRHS(
 *           const FiniteElement<dim> &temperature_fe)
 *           : local_rhs(temperature_fe.n_dofs_per_cell())
 *           , local_dof_indices(temperature_fe.n_dofs_per_cell())
 *           , matrix_for_bc(temperature_fe.n_dofs_per_cell(),
 *                           temperature_fe.n_dofs_per_cell())
 *         {}
 *       } // namespace CopyData
 *     }   // namespace Assembly
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-ThecodeBoussinesqFlowProblemcodeclasstemplate"></a> 
 * <h3>The <code>BoussinesqFlowProblem</code> class template</h3>
 *   
 
 * 
 * This is the declaration of the main class. It is very similar to @ref step_31 "step-31"
 * but there are a number differences we will comment on below.
 *   
 
 * 
 * The top of the class is essentially the same as in @ref step_31 "step-31", listing the
 * public methods and a set of private functions that do the heavy
 * lifting. Compared to @ref step_31 "step-31" there are only two additions to this
 * section: the function <code>get_cfl_number()</code> that computes the
 * maximum CFL number over all cells which we then compute the global time
 * step from, and the function <code>get_entropy_variation()</code> that is
 * used in the computation of the entropy stabilization. It is akin to the
 * <code>get_extrapolated_temperature_range()</code> we have used in @ref step_31 "step-31"
 * for this purpose, but works on the entropy instead of the temperature
 * instead.
 * 
 * @code
 *     template <int dim>
 *     class BoussinesqFlowProblem
 *     {
 *     public:
 *       struct Parameters;
 *       BoussinesqFlowProblem(Parameters &parameters);
 *       void run();
 *   
 *     private:
 *       void   setup_dofs();
 *       void   assemble_stokes_preconditioner();
 *       void   build_stokes_preconditioner();
 *       void   assemble_stokes_system();
 *       void   assemble_temperature_matrix();
 *       void   assemble_temperature_system(const double maximal_velocity);
 *       double get_maximal_velocity() const;
 *       double get_cfl_number() const;
 *       double get_entropy_variation(const double average_temperature) const;
 *       std::pair<double, double> get_extrapolated_temperature_range() const;
 *       void                      solve();
 *       void                      output_results();
 *       void                      refine_mesh(const unsigned int max_grid_level);
 *   
 *       double compute_viscosity(
 *         const std::vector<double>         &old_temperature,
 *         const std::vector<double>         &old_old_temperature,
 *         const std::vector<Tensor<1, dim>> &old_temperature_grads,
 *         const std::vector<Tensor<1, dim>> &old_old_temperature_grads,
 *         const std::vector<double>         &old_temperature_laplacians,
 *         const std::vector<double>         &old_old_temperature_laplacians,
 *         const std::vector<Tensor<1, dim>> &old_velocity_values,
 *         const std::vector<Tensor<1, dim>> &old_old_velocity_values,
 *         const std::vector<SymmetricTensor<2, dim>> &old_strain_rates,
 *         const std::vector<SymmetricTensor<2, dim>> &old_old_strain_rates,
 *         const double                                global_u_infty,
 *         const double                                global_T_variation,
 *         const double                                average_temperature,
 *         const double                                global_entropy_variation,
 *         const double                                cell_diameter) const;
 *   
 *     public:
 * @endcode
 * 
 * The first significant new component is the definition of a struct for
 * the parameters according to the discussion in the introduction. This
 * structure is initialized by reading from a parameter file during
 * construction of this object.
 * 
 * @code
 *       struct Parameters
 *       {
 *         Parameters(const std::string &parameter_filename);
 *   
 *         static void declare_parameters(ParameterHandler &prm);
 *         void        parse_parameters(ParameterHandler &prm);
 *   
 *         double end_time;
 *   
 *         unsigned int initial_global_refinement;
 *         unsigned int initial_adaptive_refinement;
 *   
 *         bool         generate_graphical_output;
 *         unsigned int graphical_output_interval;
 *   
 *         unsigned int adaptive_refinement_interval;
 *   
 *         double stabilization_alpha;
 *         double stabilization_c_R;
 *         double stabilization_beta;
 *   
 *         unsigned int stokes_velocity_degree;
 *         bool         use_locally_conservative_discretization;
 *   
 *         unsigned int temperature_degree;
 *       };
 *   
 *     private:
 *       Parameters &parameters;
 *   
 * @endcode
 * 
 * The <code>pcout</code> (for <i>%parallel <code>std::cout</code></i>)
 * object is used to simplify writing output: each MPI process can use
 * this to generate output as usual, but since each of these processes
 * will (hopefully) produce the same output it will just be replicated
 * many times over; with the ConditionalOStream class, only the output
 * generated by one MPI process will actually be printed to screen,
 * whereas the output by all the other threads will simply be forgotten.
 * 
 * @code
 *       ConditionalOStream pcout;
 *   
 * @endcode
 * 
 * The following member variables will then again be similar to those in
 * @ref step_31 "step-31" (and to other tutorial programs). As mentioned in the
 * introduction, we fully distribute computations, so we will have to use
 * the parallel::distributed::Triangulation class (see @ref step_40 "step-40") but the
 * remainder of these variables is rather standard with two exceptions:
 *     
 
 * 
 * - The <code>mapping</code> variable is used to denote a higher-order
 * polynomial mapping. As mentioned in the introduction, we use this
 * mapping when forming integrals through quadrature for all cells.
 *     
 
 * 
 * - In a bit of naming confusion, you will notice below that some of the
 * variables from namespace TrilinosWrappers are taken from namespace
 * TrilinosWrappers::MPI (such as the right hand side vectors) whereas
 * others are not (such as the various matrices). This is due to legacy
 * reasons. We will frequently have to query velocities
 * and temperatures at arbitrary quadrature points; consequently, rather
 * than importing ghost information of a vector whenever we need access
 * to degrees of freedom that are relevant locally but owned by another
 * processor, we solve linear systems in %parallel but then immediately
 * initialize a vector including ghost entries of the solution for further
 * processing. The various <code>*_solution</code> vectors are therefore
 * filled immediately after solving their respective linear system in
 * %parallel and will always contain values for all
 * @ref GlossLocallyRelevantDof "locally relevant degrees of freedom";
 * the fully distributed vectors that we obtain from the solution process
 * and that only ever contain the
 * @ref GlossLocallyOwnedDof "locally owned degrees of freedom" are
 * destroyed immediately after the solution process and after we have
 * copied the relevant values into the member variable vectors.
 * 
 * @code
 *       parallel::distributed::Triangulation<dim> triangulation;
 *       double                                    global_Omega_diameter;
 *   
 *       const MappingQ<dim> mapping;
 *   
 *       const FESystem<dim>       stokes_fe;
 *       DoFHandler<dim>           stokes_dof_handler;
 *       AffineConstraints<double> stokes_constraints;
 *   
 *       TrilinosWrappers::BlockSparseMatrix stokes_matrix;
 *       TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
 *   
 *       TrilinosWrappers::MPI::BlockVector stokes_solution;
 *       TrilinosWrappers::MPI::BlockVector old_stokes_solution;
 *       TrilinosWrappers::MPI::BlockVector stokes_rhs;
 *   
 *   
 *       const FE_Q<dim>           temperature_fe;
 *       DoFHandler<dim>           temperature_dof_handler;
 *       AffineConstraints<double> temperature_constraints;
 *   
 *       TrilinosWrappers::SparseMatrix temperature_mass_matrix;
 *       TrilinosWrappers::SparseMatrix temperature_stiffness_matrix;
 *       TrilinosWrappers::SparseMatrix temperature_matrix;
 *   
 *       TrilinosWrappers::MPI::Vector temperature_solution;
 *       TrilinosWrappers::MPI::Vector old_temperature_solution;
 *       TrilinosWrappers::MPI::Vector old_old_temperature_solution;
 *       TrilinosWrappers::MPI::Vector temperature_rhs;
 *   
 *   
 *       double       time_step;
 *       double       old_time_step;
 *       unsigned int timestep_number;
 *   
 *       std::shared_ptr<TrilinosWrappers::PreconditionAMG>    Amg_preconditioner;
 *       std::shared_ptr<TrilinosWrappers::PreconditionJacobi> Mp_preconditioner;
 *       std::shared_ptr<TrilinosWrappers::PreconditionJacobi> T_preconditioner;
 *   
 *       bool rebuild_stokes_matrix;
 *       bool rebuild_stokes_preconditioner;
 *       bool rebuild_temperature_matrices;
 *       bool rebuild_temperature_preconditioner;
 *   
 * @endcode
 * 
 * The next member variable, <code>computing_timer</code> is used to
 * conveniently account for compute time spent in certain "sections" of
 * the code that are repeatedly entered. For example, we will enter (and
 * leave) sections for Stokes matrix assembly and would like to accumulate
 * the run time spent in this section over all time steps. Every so many
 * time steps as well as at the end of the program (through the destructor
 * of the TimerOutput class) we will then produce a nice summary of the
 * times spent in the different sections into which we categorize the
 * run-time of this program.
 * 
 * @code
 *       TimerOutput computing_timer;
 *   
 * @endcode
 * 
 * After these member variables we have a number of auxiliary functions
 * that have been broken out of the ones listed above. Specifically, there
 * are first three functions that we call from <code>setup_dofs</code> and
 * then the ones that do the assembling of linear systems:
 * 
 * @code
 *       void setup_stokes_matrix(
 *         const std::vector<IndexSet> &stokes_partitioning,
 *         const std::vector<IndexSet> &stokes_relevant_partitioning);
 *       void setup_stokes_preconditioner(
 *         const std::vector<IndexSet> &stokes_partitioning,
 *         const std::vector<IndexSet> &stokes_relevant_partitioning);
 *       void setup_temperature_matrices(
 *         const IndexSet &temperature_partitioning,
 *         const IndexSet &temperature_relevant_partitioning);
 *   
 *   
 * @endcode
 * 
 * Following the @ref MTWorkStream "task-based parallelization" paradigm,
 * we split all the assembly routines into two parts: a first part that
 * can do all the calculations on a certain cell without taking care of
 * other threads, and a second part (which is writing the local data into
 * the global matrices and vectors) which can be entered by only one
 * thread at a time. In order to implement that, we provide functions for
 * each of those two steps for all the four assembly routines that we use
 * in this program. The following eight functions do exactly this:
 * 
 * @code
 *       void local_assemble_stokes_preconditioner(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         Assembly::Scratch::StokesPreconditioner<dim>         &scratch,
 *         Assembly::CopyData::StokesPreconditioner<dim>        &data);
 *   
 *       void copy_local_to_global_stokes_preconditioner(
 *         const Assembly::CopyData::StokesPreconditioner<dim> &data);
 *   
 *   
 *       void local_assemble_stokes_system(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         Assembly::Scratch::StokesSystem<dim>                 &scratch,
 *         Assembly::CopyData::StokesSystem<dim>                &data);
 *   
 *       void copy_local_to_global_stokes_system(
 *         const Assembly::CopyData::StokesSystem<dim> &data);
 *   
 *   
 *       void local_assemble_temperature_matrix(
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         Assembly::Scratch::TemperatureMatrix<dim>            &scratch,
 *         Assembly::CopyData::TemperatureMatrix<dim>           &data);
 *   
 *       void copy_local_to_global_temperature_matrix(
 *         const Assembly::CopyData::TemperatureMatrix<dim> &data);
 *   
 *   
 *   
 *       void local_assemble_temperature_rhs(
 *         const std::pair<double, double> global_T_range,
 *         const double                    global_max_velocity,
 *         const double                    global_entropy_variation,
 *         const typename DoFHandler<dim>::active_cell_iterator &cell,
 *         Assembly::Scratch::TemperatureRHS<dim>               &scratch,
 *         Assembly::CopyData::TemperatureRHS<dim>              &data);
 *   
 *       void copy_local_to_global_temperature_rhs(
 *         const Assembly::CopyData::TemperatureRHS<dim> &data);
 *   
 * @endcode
 * 
 * Finally, we forward declare a member class that we will define later on
 * and that will be used to compute a number of quantities from our
 * solution vectors that we'd like to put into the output files for
 * visualization.
 * 
 * @code
 *       class Postprocessor;
 *     };
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemclassimplementation"></a> 
 * <h3>BoussinesqFlowProblem class implementation</h3>
 * 
 
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemParameters"></a> 
 * <h4>BoussinesqFlowProblem::Parameters</h4>
 *   
 
 * 
 * Here comes the definition of the parameters for the Stokes problem. We
 * allow to set the end time for the simulation, the level of refinements
 * (both global and adaptive, which in the sum specify what maximum level
 * the cells are allowed to have), and the interval between refinements in
 * the time stepping.
 *   
 
 * 
 * Then, we let the user specify constants for the stabilization parameters
 * (as discussed in the introduction), the polynomial degree for the Stokes
 * velocity space, whether to use the locally conservative discretization
 * based on FE_DGP elements for the pressure or not (FE_Q elements for
 * pressure), and the polynomial degree for the temperature interpolation.
 *   
 
 * 
 * The constructor checks for a valid input file (if not, a file with
 * default parameters for the quantities is written), and eventually parses
 * the parameters.
 * 
 * @code
 *     template <int dim>
 *     BoussinesqFlowProblem<dim>::Parameters::Parameters(
 *       const std::string &parameter_filename)
 *       : end_time(1e8)
 *       , initial_global_refinement(2)
 *       , initial_adaptive_refinement(2)
 *       , adaptive_refinement_interval(10)
 *       , stabilization_alpha(2)
 *       , stabilization_c_R(0.11)
 *       , stabilization_beta(0.078)
 *       , stokes_velocity_degree(2)
 *       , use_locally_conservative_discretization(true)
 *       , temperature_degree(2)
 *     {
 *       ParameterHandler prm;
 *       BoussinesqFlowProblem<dim>::Parameters::declare_parameters(prm);
 *   
 *       std::ifstream parameter_file(parameter_filename);
 *   
 *       if (!parameter_file)
 *         {
 *           parameter_file.close();
 *   
 *           std::ofstream parameter_out(parameter_filename);
 *           prm.print_parameters(parameter_out, ParameterHandler::PRM);
 *   
 *           AssertThrow(
 *             false,
 *             ExcMessage(
 *               "Input parameter file <" + parameter_filename +
 *               "> not found. Creating a template file of the same name."));
 *         }
 *   
 *       prm.parse_input(parameter_file);
 *       parse_parameters(prm);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Next we have a function that declares the parameters that we expect in
 * the input file, together with their data types, default values and a
 * description:
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::Parameters::declare_parameters(
 *       ParameterHandler &prm)
 *     {
 *       prm.declare_entry("End time",
 *                         "1e8",
 *                         Patterns::Double(0),
 *                         "The end time of the simulation in years.");
 *       prm.declare_entry("Initial global refinement",
 *                         "2",
 *                         Patterns::Integer(0),
 *                         "The number of global refinement steps performed on "
 *                         "the initial coarse mesh, before the problem is first "
 *                         "solved there.");
 *       prm.declare_entry("Initial adaptive refinement",
 *                         "2",
 *                         Patterns::Integer(0),
 *                         "The number of adaptive refinement steps performed after "
 *                         "initial global refinement.");
 *       prm.declare_entry("Time steps between mesh refinement",
 *                         "10",
 *                         Patterns::Integer(1),
 *                         "The number of time steps after which the mesh is to be "
 *                         "adapted based on computed error indicators.");
 *       prm.declare_entry("Generate graphical output",
 *                         "false",
 *                         Patterns::Bool(),
 *                         "Whether graphical output is to be generated or not. "
 *                         "You may not want to get graphical output if the number "
 *                         "of processors is large.");
 *       prm.declare_entry("Time steps between graphical output",
 *                         "50",
 *                         Patterns::Integer(1),
 *                         "The number of time steps between each generation of "
 *                         "graphical output files.");
 *   
 *       prm.enter_subsection("Stabilization parameters");
 *       {
 *         prm.declare_entry("alpha",
 *                           "2",
 *                           Patterns::Double(1, 2),
 *                           "The exponent in the entropy viscosity stabilization.");
 *         prm.declare_entry("c_R",
 *                           "0.11",
 *                           Patterns::Double(0),
 *                           "The c_R factor in the entropy viscosity "
 *                           "stabilization.");
 *         prm.declare_entry("beta",
 *                           "0.078",
 *                           Patterns::Double(0),
 *                           "The beta factor in the artificial viscosity "
 *                           "stabilization. An appropriate value for 2d is 0.052 "
 *                           "and 0.078 for 3d.");
 *       }
 *       prm.leave_subsection();
 *   
 *       prm.enter_subsection("Discretization");
 *       {
 *         prm.declare_entry(
 *           "Stokes velocity polynomial degree",
 *           "2",
 *           Patterns::Integer(1),
 *           "The polynomial degree to use for the velocity variables "
 *           "in the Stokes system.");
 *         prm.declare_entry(
 *           "Temperature polynomial degree",
 *           "2",
 *           Patterns::Integer(1),
 *           "The polynomial degree to use for the temperature variable.");
 *         prm.declare_entry(
 *           "Use locally conservative discretization",
 *           "true",
 *           Patterns::Bool(),
 *           "Whether to use a Stokes discretization that is locally "
 *           "conservative at the expense of a larger number of degrees "
 *           "of freedom, or to go with a cheaper discretization "
 *           "that does not locally conserve mass (although it is "
 *           "globally conservative.");
 *       }
 *       prm.leave_subsection();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * And then we need a function that reads the contents of the
 * ParameterHandler object we get by reading the input file and puts the
 * results into variables that store the values of the parameters we have
 * previously declared:
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::Parameters::parse_parameters(
 *       ParameterHandler &prm)
 *     {
 *       end_time                  = prm.get_double("End time");
 *       initial_global_refinement = prm.get_integer("Initial global refinement");
 *       initial_adaptive_refinement =
 *         prm.get_integer("Initial adaptive refinement");
 *   
 *       adaptive_refinement_interval =
 *         prm.get_integer("Time steps between mesh refinement");
 *   
 *       generate_graphical_output = prm.get_bool("Generate graphical output");
 *       graphical_output_interval =
 *         prm.get_integer("Time steps between graphical output");
 *   
 *       prm.enter_subsection("Stabilization parameters");
 *       {
 *         stabilization_alpha = prm.get_double("alpha");
 *         stabilization_c_R   = prm.get_double("c_R");
 *         stabilization_beta  = prm.get_double("beta");
 *       }
 *       prm.leave_subsection();
 *   
 *       prm.enter_subsection("Discretization");
 *       {
 *         stokes_velocity_degree =
 *           prm.get_integer("Stokes velocity polynomial degree");
 *         temperature_degree = prm.get_integer("Temperature polynomial degree");
 *         use_locally_conservative_discretization =
 *           prm.get_bool("Use locally conservative discretization");
 *       }
 *       prm.leave_subsection();
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemBoussinesqFlowProblem"></a> 
 * <h4>BoussinesqFlowProblem::BoussinesqFlowProblem</h4>
 *   
 
 * 
 * The constructor of the problem is very similar to the constructor in
 * @ref step_31 "step-31". What is different is the %parallel communication: Trilinos uses
 * a message passing interface (MPI) for data distribution. When entering
 * the BoussinesqFlowProblem class, we have to decide how the parallelization
 * is to be done. We choose a rather simple strategy and let all processors
 * that are running the program work together, specified by the communicator
 * <code>MPI_COMM_WORLD</code>. Next, we create the output stream (as we
 * already did in @ref step_18 "step-18") that only generates output on the first MPI
 * process and is completely forgetful on all others. The implementation of
 * this idea is to check the process number when <code>pcout</code> gets a
 * true argument, and it uses the <code>std::cout</code> stream for
 * output. If we are one processor five, for instance, then we will give a
 * <code>false</code> argument to <code>pcout</code>, which means that the
 * output of that processor will not be printed. With the exception of the
 * mapping object (for which we use polynomials of degree 4) all but the
 * final member variable are exactly the same as in @ref step_31 "step-31".
 *   
 
 * 
 * This final object, the TimerOutput object, is then told to restrict
 * output to the <code>pcout</code> stream (processor 0), and then we
 * specify that we want to get a summary table at the end of the program
 * which shows us wallclock times (as opposed to CPU times). We will
 * manually also request intermediate summaries every so many time steps in
 * the <code>run()</code> function below.
 * 
 * @code
 *     template <int dim>
 *     BoussinesqFlowProblem<dim>::BoussinesqFlowProblem(Parameters &parameters_)
 *       : parameters(parameters_)
 *       , pcout(std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0))
 *       ,
 *   
 *       triangulation(MPI_COMM_WORLD,
 *                     typename Triangulation<dim>::MeshSmoothing(
 *                       Triangulation<dim>::smoothing_on_refinement |
 *                       Triangulation<dim>::smoothing_on_coarsening))
 *       ,
 *   
 *       global_Omega_diameter(0.)
 *       ,
 *   
 *       mapping(4)
 *       ,
 *   
 *       stokes_fe(FE_Q<dim>(parameters.stokes_velocity_degree) ^ dim,
 *                 (parameters.use_locally_conservative_discretization ?
 *                    static_cast<const FiniteElement<dim> &>(
 *                      FE_DGP<dim>(parameters.stokes_velocity_degree - 1)) :
 *                    static_cast<const FiniteElement<dim> &>(
 *                      FE_Q<dim>(parameters.stokes_velocity_degree - 1))))
 *       ,
 *   
 *       stokes_dof_handler(triangulation)
 *       ,
 *   
 *       temperature_fe(parameters.temperature_degree)
 *       , temperature_dof_handler(triangulation)
 *       ,
 *   
 *       time_step(0)
 *       , old_time_step(0)
 *       , timestep_number(0)
 *       , rebuild_stokes_matrix(true)
 *       , rebuild_stokes_preconditioner(true)
 *       , rebuild_temperature_matrices(true)
 *       , rebuild_temperature_preconditioner(true)
 *       ,
 *   
 *       computing_timer(MPI_COMM_WORLD,
 *                       pcout,
 *                       TimerOutput::summary,
 *                       TimerOutput::wall_times)
 *     {}
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-TheBoussinesqFlowProblemhelperfunctions"></a> 
 * <h4>The BoussinesqFlowProblem helper functions</h4>
 * 
 * <a name="step_32-BoussinesqFlowProblemget_maximal_velocity"></a> 
 * <h5>BoussinesqFlowProblem::get_maximal_velocity</h5>
 * 
 
 * 
 * Except for two small details, the function to compute the global maximum
 * of the velocity is the same as in @ref step_31 "step-31". The first detail is actually
 * common to all functions that implement loops over all cells in the
 * triangulation: When operating in %parallel, each processor can only work
 * on a chunk of cells since each processor only has a certain part of the
 * entire triangulation. This chunk of cells that we want to work on is
 * identified via a so-called <code>subdomain_id</code>, as we also did in
 * @ref step_18 "step-18". All we need to change is hence to perform the cell-related
 * operations only on cells that are owned by the current process (as
 * opposed to ghost or artificial cells), i.e. for which the subdomain id
 * equals the number of the process ID. Since this is a commonly used
 * operation, there is a shortcut for this operation: we can ask whether the
 * cell is owned by the current processor using
 * <code>cell-@>is_locally_owned()</code>.
 *   
 
 * 
 * The second difference is the way we calculate the maximum value. Before,
 * we could simply have a <code>double</code> variable that we checked
 * against on each quadrature point for each cell. Now, we have to be a bit
 * more careful since each processor only operates on a subset of
 * cells. What we do is to first let each processor calculate the maximum
 * among its cells, and then do a global communication operation
 * <code>Utilities::MPI::max</code> that computes the maximum value among
 * all the maximum values of the individual processors. MPI provides such a
 * call, but it's even simpler to use the respective function in namespace
 * Utilities::MPI using the MPI communicator object since that will do the
 * right thing even if we work without MPI and on a single machine only. The
 * call to <code>Utilities::MPI::max</code> needs two arguments, namely the
 * local maximum (input) and the MPI communicator, which is MPI_COMM_WORLD
 * in this example.
 * 
 * @code
 *     template <int dim>
 *     double BoussinesqFlowProblem<dim>::get_maximal_velocity() const
 *     {
 *       const QIterated<dim> quadrature_formula(QTrapezoid<1>(),
 *                                               parameters.stokes_velocity_degree);
 *       const unsigned int   n_q_points = quadrature_formula.size();
 *   
 *       FEValues<dim>               fe_values(mapping,
 *                               stokes_fe,
 *                               quadrature_formula,
 *                               update_values);
 *       std::vector<Tensor<1, dim>> velocity_values(n_q_points);
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *   
 *       double max_local_velocity = 0;
 *   
 *       for (const auto &cell : stokes_dof_handler.active_cell_iterators())
 *         if (cell->is_locally_owned())
 *           {
 *             fe_values.reinit(cell);
 *             fe_values[velocities].get_function_values(stokes_solution,
 *                                                       velocity_values);
 *   
 *             for (unsigned int q = 0; q < n_q_points; ++q)
 *               max_local_velocity =
 *                 std::max(max_local_velocity, velocity_values[q].norm());
 *           }
 *   
 *       return Utilities::MPI::max(max_local_velocity, MPI_COMM_WORLD);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemget_cfl_number"></a> 
 * <h5>BoussinesqFlowProblem::get_cfl_number</h5>
 * 
 
 * 
 * The next function does something similar, but we now compute the CFL
 * number, i.e., maximal velocity on a cell divided by the cell
 * diameter. This number is necessary to determine the time step size, as we
 * use a semi-explicit time stepping scheme for the temperature equation
 * (see @ref step_31 "step-31" for a discussion). We compute it in the same way as above:
 * Compute the local maximum over all locally owned cells, then exchange it
 * via MPI to find the global maximum.
 * 
 * @code
 *     template <int dim>
 *     double BoussinesqFlowProblem<dim>::get_cfl_number() const
 *     {
 *       const QIterated<dim> quadrature_formula(QTrapezoid<1>(),
 *                                               parameters.stokes_velocity_degree);
 *       const unsigned int   n_q_points = quadrature_formula.size();
 *   
 *       FEValues<dim>               fe_values(mapping,
 *                               stokes_fe,
 *                               quadrature_formula,
 *                               update_values);
 *       std::vector<Tensor<1, dim>> velocity_values(n_q_points);
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *   
 *       double max_local_cfl = 0;
 *   
 *       for (const auto &cell : stokes_dof_handler.active_cell_iterators())
 *         if (cell->is_locally_owned())
 *           {
 *             fe_values.reinit(cell);
 *             fe_values[velocities].get_function_values(stokes_solution,
 *                                                       velocity_values);
 *   
 *             double max_local_velocity = 1e-10;
 *             for (unsigned int q = 0; q < n_q_points; ++q)
 *               max_local_velocity =
 *                 std::max(max_local_velocity, velocity_values[q].norm());
 *             max_local_cfl =
 *               std::max(max_local_cfl, max_local_velocity / cell->diameter());
 *           }
 *   
 *       return Utilities::MPI::max(max_local_cfl, MPI_COMM_WORLD);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemget_entropy_variation"></a> 
 * <h5>BoussinesqFlowProblem::get_entropy_variation</h5>
 * 
 
 * 
 * Next comes the computation of the global entropy variation
 * @f$\|E(T)-\bar{E}(T)\|_\infty@f$ where the entropy @f$E@f$ is defined as
 * discussed in the introduction.  This is needed for the evaluation of the
 * stabilization in the temperature equation as explained in the
 * introduction. The entropy variation is actually only needed if we use
 * @f$\alpha=2@f$ as a power in the residual computation. The infinity norm is
 * computed by the maxima over quadrature points, as usual in discrete
 * computations.
 *   
 
 * 
 * In order to compute this quantity, we first have to find the
 * space-average @f$\bar{E}(T)@f$ and then evaluate the maximum. However, that
 * means that we would need to perform two loops. We can avoid the overhead
 * by noting that @f$\|E(T)-\bar{E}(T)\|_\infty =
 * \max\big(E_{\textrm{max}}(T)-\bar{E}(T),
 * \bar{E}(T)-E_{\textrm{min}}(T)\big)@f$, i.e., the maximum out of the
 * deviation from the average entropy in positive and negative
 * directions. The four quantities we need for the latter formula (maximum
 * entropy, minimum entropy, average entropy, area) can all be evaluated in
 * the same loop over all cells, so we choose this simpler variant.
 * 
 * @code
 *     template <int dim>
 *     double BoussinesqFlowProblem<dim>::get_entropy_variation(
 *       const double average_temperature) const
 *     {
 *       if (parameters.stabilization_alpha != 2)
 *         return 1.;
 *   
 *       const QGauss<dim>  quadrature_formula(parameters.temperature_degree + 1);
 *       const unsigned int n_q_points = quadrature_formula.size();
 *   
 *       FEValues<dim>       fe_values(temperature_fe,
 *                               quadrature_formula,
 *                               update_values | update_JxW_values);
 *       std::vector<double> old_temperature_values(n_q_points);
 *       std::vector<double> old_old_temperature_values(n_q_points);
 *   
 * @endcode
 * 
 * In the two functions above we computed the maximum of numbers that were
 * all non-negative, so we knew that zero was certainly a lower bound. On
 * the other hand, here we need to find the maximum deviation from the
 * average value, i.e., we will need to know the maximal and minimal
 * values of the entropy for which we don't a priori know the sign.
 *     
 
 * 
 * To compute it, we can therefore start with the largest and smallest
 * possible values we can store in a double precision number: The minimum
 * is initialized with a bigger and the maximum with a smaller number than
 * any one that is going to appear. We are then guaranteed that these
 * numbers will be overwritten in the loop on the first cell or, if this
 * processor does not own any cells, in the communication step at the
 * latest. The following loop then computes the minimum and maximum local
 * entropy as well as keeps track of the area/volume of the part of the
 * domain we locally own and the integral over the entropy on it:
 * 
 * @code
 *       double min_entropy = std::numeric_limits<double>::max(),
 *              max_entropy = -std::numeric_limits<double>::max(), area = 0,
 *              entropy_integrated = 0;
 *   
 *       for (const auto &cell : temperature_dof_handler.active_cell_iterators())
 *         if (cell->is_locally_owned())
 *           {
 *             fe_values.reinit(cell);
 *             fe_values.get_function_values(old_temperature_solution,
 *                                           old_temperature_values);
 *             fe_values.get_function_values(old_old_temperature_solution,
 *                                           old_old_temperature_values);
 *             for (unsigned int q = 0; q < n_q_points; ++q)
 *               {
 *                 const double T =
 *                   (old_temperature_values[q] + old_old_temperature_values[q]) / 2;
 *                 const double entropy =
 *                   ((T - average_temperature) * (T - average_temperature));
 *   
 *                 min_entropy = std::min(min_entropy, entropy);
 *                 max_entropy = std::max(max_entropy, entropy);
 *                 area += fe_values.JxW(q);
 *                 entropy_integrated += fe_values.JxW(q) * entropy;
 *               }
 *           }
 *   
 * @endcode
 * 
 * Now we only need to exchange data between processors: we need to sum
 * the two integrals (<code>area</code>, <code>entropy_integrated</code>),
 * and get the extrema for maximum and minimum. We could do this through
 * four different data exchanges, but we can it with two:
 * Utilities::MPI::sum also exists in a variant that takes an array of
 * values that are all to be summed up. And we can also utilize the
 * Utilities::MPI::max function by realizing that forming the minimum over
 * the minimal entropies equals forming the negative of the maximum over
 * the negative of the minimal entropies; this maximum can then be
 * combined with forming the maximum over the maximal entropies.
 * 
 * @code
 *       const double local_sums[2]   = {entropy_integrated, area},
 *                    local_maxima[2] = {-min_entropy, max_entropy};
 *       double global_sums[2], global_maxima[2];
 *   
 *       Utilities::MPI::sum(local_sums, MPI_COMM_WORLD, global_sums);
 *       Utilities::MPI::max(local_maxima, MPI_COMM_WORLD, global_maxima);
 *   
 * @endcode
 * 
 * Having computed everything this way, we can then compute the average
 * entropy and find the @f$L^\infty@f$ norm by taking the larger of the
 * deviation of the maximum or minimum from the average:
 * 
 * @code
 *       const double average_entropy = global_sums[0] / global_sums[1];
 *       const double entropy_diff    = std::max(global_maxima[1] - average_entropy,
 *                                            average_entropy - (-global_maxima[0]));
 *       return entropy_diff;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemget_extrapolated_temperature_range"></a> 
 * <h5>BoussinesqFlowProblem::get_extrapolated_temperature_range</h5>
 * 
 
 * 
 * The next function computes the minimal and maximal value of the
 * extrapolated temperature over the entire domain. Again, this is only a
 * slightly modified version of the respective function in @ref step_31 "step-31". As in
 * the function above, we collect local minima and maxima and then compute
 * the global extrema using the same trick as above.
 *   
 
 * 
 * As already discussed in @ref step_31 "step-31", the function needs to distinguish
 * between the first and all following time steps because it uses a higher
 * order temperature extrapolation scheme when at least two previous time
 * steps are available.
 * 
 * @code
 *     template <int dim>
 *     std::pair<double, double>
 *     BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range() const
 *     {
 *       const QIterated<dim> quadrature_formula(QTrapezoid<1>(),
 *                                               parameters.temperature_degree);
 *       const unsigned int   n_q_points = quadrature_formula.size();
 *   
 *       FEValues<dim>       fe_values(mapping,
 *                               temperature_fe,
 *                               quadrature_formula,
 *                               update_values);
 *       std::vector<double> old_temperature_values(n_q_points);
 *       std::vector<double> old_old_temperature_values(n_q_points);
 *   
 *       double min_local_temperature = std::numeric_limits<double>::max(),
 *              max_local_temperature = -std::numeric_limits<double>::max();
 *   
 *       if (timestep_number != 0)
 *         {
 *           for (const auto &cell : temperature_dof_handler.active_cell_iterators())
 *             if (cell->is_locally_owned())
 *               {
 *                 fe_values.reinit(cell);
 *                 fe_values.get_function_values(old_temperature_solution,
 *                                               old_temperature_values);
 *                 fe_values.get_function_values(old_old_temperature_solution,
 *                                               old_old_temperature_values);
 *   
 *                 for (unsigned int q = 0; q < n_q_points; ++q)
 *                   {
 *                     const double temperature =
 *                       (1. + time_step / old_time_step) *
 *                         old_temperature_values[q] -
 *                       time_step / old_time_step * old_old_temperature_values[q];
 *   
 *                     min_local_temperature =
 *                       std::min(min_local_temperature, temperature);
 *                     max_local_temperature =
 *                       std::max(max_local_temperature, temperature);
 *                   }
 *               }
 *         }
 *       else
 *         {
 *           for (const auto &cell : temperature_dof_handler.active_cell_iterators())
 *             if (cell->is_locally_owned())
 *               {
 *                 fe_values.reinit(cell);
 *                 fe_values.get_function_values(old_temperature_solution,
 *                                               old_temperature_values);
 *   
 *                 for (unsigned int q = 0; q < n_q_points; ++q)
 *                   {
 *                     const double temperature = old_temperature_values[q];
 *   
 *                     min_local_temperature =
 *                       std::min(min_local_temperature, temperature);
 *                     max_local_temperature =
 *                       std::max(max_local_temperature, temperature);
 *                   }
 *               }
 *         }
 *   
 *       double local_extrema[2] = {-min_local_temperature, max_local_temperature};
 *       double global_extrema[2];
 *       Utilities::MPI::max(local_extrema, MPI_COMM_WORLD, global_extrema);
 *   
 *       return std::make_pair(-global_extrema[0], global_extrema[1]);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemcompute_viscosity"></a> 
 * <h5>BoussinesqFlowProblem::compute_viscosity</h5>
 * 
 
 * 
 * The function that calculates the viscosity is purely local and so needs
 * no communication at all. It is mostly the same as in @ref step_31 "step-31" but with an
 * updated formulation of the viscosity if @f$\alpha=2@f$ is chosen:
 * 
 * @code
 *     template <int dim>
 *     double BoussinesqFlowProblem<dim>::compute_viscosity(
 *       const std::vector<double>                  &old_temperature,
 *       const std::vector<double>                  &old_old_temperature,
 *       const std::vector<Tensor<1, dim>>          &old_temperature_grads,
 *       const std::vector<Tensor<1, dim>>          &old_old_temperature_grads,
 *       const std::vector<double>                  &old_temperature_laplacians,
 *       const std::vector<double>                  &old_old_temperature_laplacians,
 *       const std::vector<Tensor<1, dim>>          &old_velocity_values,
 *       const std::vector<Tensor<1, dim>>          &old_old_velocity_values,
 *       const std::vector<SymmetricTensor<2, dim>> &old_strain_rates,
 *       const std::vector<SymmetricTensor<2, dim>> &old_old_strain_rates,
 *       const double                                global_u_infty,
 *       const double                                global_T_variation,
 *       const double                                average_temperature,
 *       const double                                global_entropy_variation,
 *       const double                                cell_diameter) const
 *     {
 *       if (global_u_infty == 0)
 *         return 5e-3 * cell_diameter;
 *   
 *       const unsigned int n_q_points = old_temperature.size();
 *   
 *       double max_residual = 0;
 *       double max_velocity = 0;
 *   
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           const Tensor<1, dim> u =
 *             (old_velocity_values[q] + old_old_velocity_values[q]) / 2;
 *   
 *           const SymmetricTensor<2, dim> strain_rate =
 *             (old_strain_rates[q] + old_old_strain_rates[q]) / 2;
 *   
 *           const double T = (old_temperature[q] + old_old_temperature[q]) / 2;
 *           const double dT_dt =
 *             (old_temperature[q] - old_old_temperature[q]) / old_time_step;
 *           const double u_grad_T =
 *             u * (old_temperature_grads[q] + old_old_temperature_grads[q]) / 2;
 *   
 *           const double kappa_Delta_T =
 *             EquationData::kappa *
 *             (old_temperature_laplacians[q] + old_old_temperature_laplacians[q]) /
 *             2;
 *           const double gamma =
 *             ((EquationData::radiogenic_heating * EquationData::density(T) +
 *               2 * EquationData::eta * strain_rate * strain_rate) /
 *              (EquationData::density(T) * EquationData::specific_heat));
 *   
 *           double residual = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma);
 *           if (parameters.stabilization_alpha == 2)
 *             residual *= std::abs(T - average_temperature);
 *   
 *           max_residual = std::max(residual, max_residual);
 *           max_velocity = std::max(std::sqrt(u * u), max_velocity);
 *         }
 *   
 *       const double max_viscosity =
 *         (parameters.stabilization_beta * max_velocity * cell_diameter);
 *       if (timestep_number == 0)
 *         return max_viscosity;
 *       else
 *         {
 *           Assert(old_time_step > 0, ExcInternalError());
 *   
 *           double entropy_viscosity;
 *           if (parameters.stabilization_alpha == 2)
 *             entropy_viscosity =
 *               (parameters.stabilization_c_R * cell_diameter * cell_diameter *
 *                max_residual / global_entropy_variation);
 *           else
 *             entropy_viscosity =
 *               (parameters.stabilization_c_R * cell_diameter *
 *                global_Omega_diameter * max_velocity * max_residual /
 *                (global_u_infty * global_T_variation));
 *   
 *           return std::min(max_viscosity, entropy_viscosity);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-TheBoussinesqFlowProblemsetupfunctions"></a> 
 * <h4>The BoussinesqFlowProblem setup functions</h4>
 * 
 
 * 
 * The following three functions set up the Stokes matrix, the matrix used
 * for the Stokes preconditioner, and the temperature matrix. The code is
 * mostly the same as in @ref step_31 "step-31", but it has been broken out into three
 * functions of their own for simplicity.
 *   
 
 * 
 * The main functional difference between the code here and that in @ref step_31 "step-31"
 * is that the matrices we want to set up are distributed across multiple
 * processors. Since we still want to build up the sparsity pattern first
 * for efficiency reasons, we could continue to build the <i>entire</i>
 * sparsity pattern as a BlockDynamicSparsityPattern, as we did in
 * @ref step_31 "step-31". However, that would be inefficient: every processor would build
 * the same sparsity pattern, but only initialize a small part of the matrix
 * using it. It also violates the principle that every processor should only
 * work on those cells it owns (and, if necessary the layer of ghost cells
 * around it).
 *   
 
 * 
 * Rather, we use an object of type TrilinosWrappers::BlockSparsityPattern,
 * which is (obviously) a wrapper around a sparsity pattern object provided
 * by Trilinos. The advantage is that the Trilinos sparsity pattern class
 * can communicate across multiple processors: if this processor fills in
 * all the nonzero entries that result from the cells it owns, and every
 * other processor does so as well, then at the end after some MPI
 * communication initiated by the <code>compress()</code> call, we will have
 * the globally assembled sparsity pattern available with which the global
 * matrix can be initialized.
 *   
 
 * 
 * There is one important aspect when initializing Trilinos sparsity
 * patterns in parallel: In addition to specifying the locally owned rows
 * and columns of the matrices via the @p stokes_partitioning index set, we
 * also supply information about all the rows we are possibly going to write
 * into when assembling on a certain processor. The set of locally relevant
 * rows contains all such rows (possibly also a few unnecessary ones, but it
 * is difficult to find the exact row indices before actually getting
 * indices on all cells and resolving constraints). This additional
 * information allows to exactly determine the structure for the
 * off-processor data found during assembly. While Trilinos matrices are
 * able to collect this information on the fly as well (when initializing
 * them from some other reinit method), it is less efficient and leads to
 * problems when assembling matrices with multiple threads. In this program,
 * we pessimistically assume that only one processor at a time can write
 * into the matrix while assembly (whereas the computation is parallel),
 * which is fine for Trilinos matrices. In practice, one can do better by
 * hinting WorkStream at cells that do not share vertices, allowing for
 * parallelism among those cells (see the graph coloring algorithms and
 * WorkStream with colored iterators argument). However, that only works
 * when only one MPI processor is present because Trilinos' internal data
 * structures for accumulating off-processor data on the fly are not thread
 * safe. With the initialization presented here, there is no such problem
 * and one could safely introduce graph coloring for this algorithm.
 *   
 
 * 
 * The only other change we need to make is to tell the
 * DoFTools::make_sparsity_pattern() function that it is only supposed to
 * work on a subset of cells, namely the ones whose
 * <code>subdomain_id</code> equals the number of the current processor, and
 * to ignore all other cells.
 *   
 
 * 
 * This strategy is replicated across all three of the following functions.
 *   
 
 * 
 * Note that Trilinos matrices store the information contained in the
 * sparsity patterns, so we can safely release the <code>sp</code> variable
 * once the matrix has been given the sparsity structure.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::setup_stokes_matrix(
 *       const std::vector<IndexSet> &stokes_partitioning,
 *       const std::vector<IndexSet> &stokes_relevant_partitioning)
 *     {
 *       stokes_matrix.clear();
 *   
 *       TrilinosWrappers::BlockSparsityPattern sp(stokes_partitioning,
 *                                                 stokes_partitioning,
 *                                                 stokes_relevant_partitioning,
 *                                                 MPI_COMM_WORLD);
 *   
 *       Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
 *       for (unsigned int c = 0; c < dim + 1; ++c)
 *         for (unsigned int d = 0; d < dim + 1; ++d)
 *           if (!((c == dim) && (d == dim)))
 *             coupling[c][d] = DoFTools::always;
 *           else
 *             coupling[c][d] = DoFTools::none;
 *   
 *       DoFTools::make_sparsity_pattern(stokes_dof_handler,
 *                                       coupling,
 *                                       sp,
 *                                       stokes_constraints,
 *                                       false,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         MPI_COMM_WORLD));
 *       sp.compress();
 *   
 *       stokes_matrix.reinit(sp);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::setup_stokes_preconditioner(
 *       const std::vector<IndexSet> &stokes_partitioning,
 *       const std::vector<IndexSet> &stokes_relevant_partitioning)
 *     {
 *       Amg_preconditioner.reset();
 *       Mp_preconditioner.reset();
 *   
 *       stokes_preconditioner_matrix.clear();
 *   
 *       TrilinosWrappers::BlockSparsityPattern sp(stokes_partitioning,
 *                                                 stokes_partitioning,
 *                                                 stokes_relevant_partitioning,
 *                                                 MPI_COMM_WORLD);
 *   
 *       Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
 *       for (unsigned int c = 0; c < dim + 1; ++c)
 *         for (unsigned int d = 0; d < dim + 1; ++d)
 *           if (c == d)
 *             coupling[c][d] = DoFTools::always;
 *           else
 *             coupling[c][d] = DoFTools::none;
 *   
 *       DoFTools::make_sparsity_pattern(stokes_dof_handler,
 *                                       coupling,
 *                                       sp,
 *                                       stokes_constraints,
 *                                       false,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         MPI_COMM_WORLD));
 *       sp.compress();
 *   
 *       stokes_preconditioner_matrix.reinit(sp);
 *     }
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::setup_temperature_matrices(
 *       const IndexSet &temperature_partitioner,
 *       const IndexSet &temperature_relevant_partitioner)
 *     {
 *       T_preconditioner.reset();
 *       temperature_mass_matrix.clear();
 *       temperature_stiffness_matrix.clear();
 *       temperature_matrix.clear();
 *   
 *       TrilinosWrappers::SparsityPattern sp(temperature_partitioner,
 *                                            temperature_partitioner,
 *                                            temperature_relevant_partitioner,
 *                                            MPI_COMM_WORLD);
 *       DoFTools::make_sparsity_pattern(temperature_dof_handler,
 *                                       sp,
 *                                       temperature_constraints,
 *                                       false,
 *                                       Utilities::MPI::this_mpi_process(
 *                                         MPI_COMM_WORLD));
 *       sp.compress();
 *   
 *       temperature_matrix.reinit(sp);
 *       temperature_mass_matrix.reinit(sp);
 *       temperature_stiffness_matrix.reinit(sp);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The remainder of the setup function (after splitting out the three
 * functions above) mostly has to deal with the things we need to do for
 * parallelization across processors. Because setting all of this up is a
 * significant compute time expense of the program, we put everything we do
 * here into a timer group so that we can get summary information about the
 * fraction of time spent in this part of the program at its end.
 *   
 
 * 
 * At the top as usual we enumerate degrees of freedom and sort them by
 * component/block, followed by writing their numbers to the screen from
 * processor zero. The DoFHandler::distributed_dofs() function, when applied
 * to a parallel::distributed::Triangulation object, sorts degrees of
 * freedom in such a way that all degrees of freedom associated with
 * subdomain zero come before all those associated with subdomain one,
 * etc. For the Stokes part, this entails, however, that velocities and
 * pressures become intermixed, but this is trivially solved by sorting
 * again by blocks; it is worth noting that this latter operation leaves the
 * relative ordering of all velocities and pressures alone, i.e. within the
 * velocity block we will still have all those associated with subdomain
 * zero before all velocities associated with subdomain one, etc. This is
 * important since we store each of the blocks of this matrix distributed
 * across all processors and want this to be done in such a way that each
 * processor stores that part of the matrix that is roughly equal to the
 * degrees of freedom located on those cells that it will actually work on.
 *   
 
 * 
 * When printing the numbers of degrees of freedom, note that these numbers
 * are going to be large if we use many processors. Consequently, we let the
 * stream put a comma separator in between every three digits. The state of
 * the stream, using the locale, is saved from before to after this
 * operation. While slightly opaque, the code works because the default
 * locale (which we get using the constructor call
 * <code>std::locale("")</code>) implies printing numbers with a comma
 * separator for every third digit (i.e., thousands, millions, billions).
 *   
 
 * 
 * In this function as well as many below, we measure how much time
 * we spend here and collect that in a section called "Setup dof
 * systems" across function invocations. This is done using an
 * TimerOutput::Scope object that gets a timer going in the section
 * with above name of the `computing_timer` object upon construction
 * of the local variable; the timer is stopped again when the
 * destructor of the `timing_section` variable is called.  This, of
 * course, happens either at the end of the function, or if we leave
 * the function through a `return` statement or when an exception is
 * thrown somewhere -- in other words, whenever we leave this
 * function in any way. The use of such "scope" objects therefore
 * makes sure that we do not have to manually add code that tells
 * the timer to stop at every location where this function may be
 * left.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::setup_dofs()
 *     {
 *       TimerOutput::Scope timing_section(computing_timer, "Setup dof systems");
 *   
 *       stokes_dof_handler.distribute_dofs(stokes_fe);
 *   
 *       std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
 *       stokes_sub_blocks[dim] = 1;
 *       DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);
 *   
 *       temperature_dof_handler.distribute_dofs(temperature_fe);
 *   
 *       const std::vector<types::global_dof_index> stokes_dofs_per_block =
 *         DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks);
 *   
 *       const types::global_dof_index n_u = stokes_dofs_per_block[0],
 *                                     n_p = stokes_dofs_per_block[1],
 *                                     n_T = temperature_dof_handler.n_dofs();
 *   
 *       std::locale s = pcout.get_stream().getloc();
 *       pcout.get_stream().imbue(std::locale(""));
 *       pcout << "Number of active cells: " << triangulation.n_global_active_cells()
 *             << " (on " << triangulation.n_levels() << " levels)" << std::endl
 *             << "Number of degrees of freedom: " << n_u + n_p + n_T << " (" << n_u
 *             << '+' << n_p << '+' << n_T << ')' << std::endl
 *             << std::endl;
 *       pcout.get_stream().imbue(s);
 *   
 *   
 * @endcode
 * 
 * After this, we have to set up the various partitioners (of type
 * <code>IndexSet</code>, see the introduction) that describe which parts
 * of each matrix or vector will be stored where, then call the functions
 * that actually set up the matrices, and at the end also resize the
 * various vectors we keep around in this program.
 * 
 
 * 
 * 
 * @code
 *       const IndexSet &stokes_locally_owned_index_set =
 *         stokes_dof_handler.locally_owned_dofs();
 *       const IndexSet stokes_locally_relevant_set =
 *         DoFTools::extract_locally_relevant_dofs(stokes_dof_handler);
 *   
 *       std::vector<IndexSet> stokes_partitioning;
 *       stokes_partitioning.push_back(
 *         stokes_locally_owned_index_set.get_view(0, n_u));
 *       stokes_partitioning.push_back(
 *         stokes_locally_owned_index_set.get_view(n_u, n_u + n_p));
 *   
 *       std::vector<IndexSet> stokes_relevant_partitioning;
 *       stokes_relevant_partitioning.push_back(
 *         stokes_locally_relevant_set.get_view(0, n_u));
 *       stokes_relevant_partitioning.push_back(
 *         stokes_locally_relevant_set.get_view(n_u, n_u + n_p));
 *   
 *       const IndexSet temperature_partitioning =
 *         temperature_dof_handler.locally_owned_dofs();
 *       const IndexSet temperature_relevant_partitioning =
 *         DoFTools::extract_locally_relevant_dofs(temperature_dof_handler);
 *   
 * @endcode
 * 
 * Following this, we can compute constraints for the solution vectors,
 * including hanging node constraints and homogeneous and inhomogeneous
 * boundary values for the Stokes and temperature fields. Note that as for
 * everything else, the constraint objects can not hold <i>all</i>
 * constraints on every processor. Rather, each processor needs to store
 * only those that are actually necessary for correctness given that it
 * only assembles linear systems on cells it owns. As discussed in the
 * @ref distributed_paper "this paper", the set of constraints we need to
 * know about is exactly the set of constraints on all locally relevant
 * degrees of freedom, so this is what we use to initialize the constraint
 * objects.
 * 
 * @code
 *       {
 *         stokes_constraints.clear();
 *         stokes_constraints.reinit(stokes_locally_owned_index_set,
 *                                   stokes_locally_relevant_set);
 *   
 *         DoFTools::make_hanging_node_constraints(stokes_dof_handler,
 *                                                 stokes_constraints);
 *   
 *         const FEValuesExtractors::Vector velocity_components(0);
 *         VectorTools::interpolate_boundary_values(
 *           stokes_dof_handler,
 *           0,
 *           Functions::ZeroFunction<dim>(dim + 1),
 *           stokes_constraints,
 *           stokes_fe.component_mask(velocity_components));
 *   
 *         std::set<types::boundary_id> no_normal_flux_boundaries;
 *         no_normal_flux_boundaries.insert(1);
 *         VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler,
 *                                                         0,
 *                                                         no_normal_flux_boundaries,
 *                                                         stokes_constraints,
 *                                                         mapping);
 *         stokes_constraints.close();
 *       }
 *       {
 *         temperature_constraints.clear();
 *         temperature_constraints.reinit(temperature_partitioning,
 *                                        temperature_relevant_partitioning);
 *   
 *         DoFTools::make_hanging_node_constraints(temperature_dof_handler,
 *                                                 temperature_constraints);
 *         VectorTools::interpolate_boundary_values(
 *           temperature_dof_handler,
 *           0,
 *           EquationData::TemperatureInitialValues<dim>(),
 *           temperature_constraints);
 *         VectorTools::interpolate_boundary_values(
 *           temperature_dof_handler,
 *           1,
 *           EquationData::TemperatureInitialValues<dim>(),
 *           temperature_constraints);
 *         temperature_constraints.close();
 *       }
 *   
 * @endcode
 * 
 * All this done, we can then initialize the various matrix and vector
 * objects to their proper sizes. At the end, we also record that all
 * matrices and preconditioners have to be re-computed at the beginning of
 * the next time step. Note how we initialize the vectors for the Stokes
 * and temperature right hand sides: These are writable vectors (last
 * boolean argument set to @p true) that have the correct one-to-one
 * partitioning of locally owned elements but are still given the relevant
 * partitioning for means of figuring out the vector entries that are
 * going to be set right away. As for matrices, this allows for writing
 * local contributions into the vector with multiple threads (always
 * assuming that the same vector entry is not accessed by multiple threads
 * at the same time). The other vectors only allow for read access of
 * individual elements, including ghosts, but are not suitable for
 * solvers.
 * 
 * @code
 *       setup_stokes_matrix(stokes_partitioning, stokes_relevant_partitioning);
 *       setup_stokes_preconditioner(stokes_partitioning,
 *                                   stokes_relevant_partitioning);
 *       setup_temperature_matrices(temperature_partitioning,
 *                                  temperature_relevant_partitioning);
 *   
 *       stokes_rhs.reinit(stokes_partitioning,
 *                         stokes_relevant_partitioning,
 *                         MPI_COMM_WORLD,
 *                         true);
 *       stokes_solution.reinit(stokes_relevant_partitioning, MPI_COMM_WORLD);
 *       old_stokes_solution.reinit(stokes_solution);
 *   
 *       temperature_rhs.reinit(temperature_partitioning,
 *                              temperature_relevant_partitioning,
 *                              MPI_COMM_WORLD,
 *                              true);
 *       temperature_solution.reinit(temperature_relevant_partitioning,
 *                                   MPI_COMM_WORLD);
 *       old_temperature_solution.reinit(temperature_solution);
 *       old_old_temperature_solution.reinit(temperature_solution);
 *   
 *       rebuild_stokes_matrix              = true;
 *       rebuild_stokes_preconditioner      = true;
 *       rebuild_temperature_matrices       = true;
 *       rebuild_temperature_preconditioner = true;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-TheBoussinesqFlowProblemassemblyfunctions"></a> 
 * <h4>The BoussinesqFlowProblem assembly functions</h4>
 *   
 
 * 
 * Following the discussion in the introduction and in the @ref threads
 * topic, we split the assembly functions into different parts:
 *   
 
 * 
 * <ul> <li> The local calculations of matrices and right hand sides, given
 * a certain cell as input (these functions are named
 * <code>local_assemble_*</code> below). The resulting function is, in other
 * words, essentially the body of the loop over all cells in @ref step_31 "step-31". Note,
 * however, that these functions store the result from the local
 * calculations in variables of classes from the CopyData namespace.
 *   
 
 * 
 * <li>These objects are then given to the second step which writes the
 * local data into the global data structures (these functions are named
 * <code>copy_local_to_global_*</code> below). These functions are pretty
 * trivial.
 *   
 
 * 
 * <li>These two subfunctions are then used in the respective assembly
 * routine (called <code>assemble_*</code> below), where a WorkStream object
 * is set up and runs over all the cells that belong to the processor's
 * subdomain.  </ul>
 * 
 
 * 
 * 
 * <a name="step_32-Stokespreconditionerassembly"></a> 
 * <h5>Stokes preconditioner assembly</h5>
 *   
 
 * 
 * Let us start with the functions that builds the Stokes
 * preconditioner. The first two of these are pretty trivial, given the
 * discussion above. Note in particular that the main point in using the
 * scratch data object is that we want to avoid allocating any objects on
 * the free space each time we visit a new cell. As a consequence, the
 * assembly function below only has automatic local variables, and
 * everything else is accessed through the scratch data object, which is
 * allocated only once before we start the loop over all cells:
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::local_assemble_stokes_preconditioner(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       Assembly::Scratch::StokesPreconditioner<dim>         &scratch,
 *       Assembly::CopyData::StokesPreconditioner<dim>        &data)
 *     {
 *       const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
 *       const unsigned int n_q_points =
 *         scratch.stokes_fe_values.n_quadrature_points;
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *       const FEValuesExtractors::Scalar pressure(dim);
 *   
 *       scratch.stokes_fe_values.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 *   
 *       data.local_matrix = 0;
 *   
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               scratch.grad_phi_u[k] =
 *                 scratch.stokes_fe_values[velocities].gradient(k, q);
 *               scratch.phi_p[k] = scratch.stokes_fe_values[pressure].value(k, q);
 *             }
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *               data.local_matrix(i, j) +=
 *                 (EquationData::eta *
 *                    scalar_product(scratch.grad_phi_u[i], scratch.grad_phi_u[j]) +
 *                  (1. / EquationData::eta) * EquationData::pressure_scaling *
 *                    EquationData::pressure_scaling *
 *                    (scratch.phi_p[i] * scratch.phi_p[j])) *
 *                 scratch.stokes_fe_values.JxW(q);
 *         }
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::copy_local_to_global_stokes_preconditioner(
 *       const Assembly::CopyData::StokesPreconditioner<dim> &data)
 *     {
 *       stokes_constraints.distribute_local_to_global(data.local_matrix,
 *                                                     data.local_dof_indices,
 *                                                     stokes_preconditioner_matrix);
 *     }
 *   
 *   
 * @endcode
 * 
 * Now for the function that actually puts things together, using the
 * WorkStream functions.  WorkStream::run needs a start and end iterator to
 * enumerate the cells it is supposed to work on. Typically, one would use
 * DoFHandler::begin_active() and DoFHandler::end() for that but here we
 * actually only want the subset of cells that in fact are owned by the
 * current processor. This is where the FilteredIterator class comes into
 * play: you give it a range of cells and it provides an iterator that only
 * iterates over that subset of cells that satisfy a certain predicate (a
 * predicate is a function of one argument that either returns true or
 * false). The predicate we use here is IteratorFilters::LocallyOwnedCell,
 * i.e., it returns true exactly if the cell is owned by the current
 * processor. The resulting iterator range is then exactly what we need.
 *   
 
 * 
 * With this obstacle out of the way, we call the WorkStream::run
 * function with this set of cells, scratch and copy objects, and
 * with pointers to two functions: the local assembly and
 * copy-local-to-global function. These functions need to have very
 * specific signatures: three arguments in the first and one
 * argument in the latter case (see the documentation of the
 * WorkStream::run function for the meaning of these arguments).
 * Note how we use a lambda functions to
 * create a function object that satisfies this requirement. It uses
 * function arguments for the local assembly function that specify
 * cell, scratch data, and copy data, as well as function argument
 * for the copy function that expects the
 * data to be written into the global matrix (also see the discussion in
 * @ref step_13 "step-13"'s <code>assemble_linear_system()</code> function). On the other
 * hand, the implicit zeroth argument of member functions (namely
 * the <code>this</code> pointer of the object on which that member
 * function is to operate on) is <i>bound</i> to the
 * <code>this</code> pointer of the current function and is captured. The
 * WorkStream::run function, as a consequence, does not need to know
 * anything about the object these functions work on.
 *   
 
 * 
 * When the WorkStream is executed, it will create several local assembly
 * routines of the first kind for several cells and let some available
 * processors work on them. The function that needs to be synchronized,
 * i.e., the write operation into the global matrix, however, is executed by
 * only one thread at a time in the prescribed order. Of course, this only
 * holds for the parallelization on a single MPI process. Different MPI
 * processes will have their own WorkStream objects and do that work
 * completely independently (and in different memory spaces). In a
 * distributed calculation, some data will accumulate at degrees of freedom
 * that are not owned by the respective processor. It would be inefficient
 * to send data around every time we encounter such a dof. What happens
 * instead is that the Trilinos sparse matrix will keep that data and send
 * it to the owner at the end of assembly, by calling the
 * <code>compress()</code> command.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner()
 *     {
 *       stokes_preconditioner_matrix = 0;
 *   
 *       const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree + 1);
 *   
 *       using CellFilter =
 *         FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>;
 *   
 *       auto worker =
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                Assembly::Scratch::StokesPreconditioner<dim>         &scratch,
 *                Assembly::CopyData::StokesPreconditioner<dim>        &data) {
 *           this->local_assemble_stokes_preconditioner(cell, scratch, data);
 *         };
 *   
 *       auto copier =
 *         [this](const Assembly::CopyData::StokesPreconditioner<dim> &data) {
 *           this->copy_local_to_global_stokes_preconditioner(data);
 *         };
 *   
 *       WorkStream::run(CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                                  stokes_dof_handler.begin_active()),
 *                       CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                                  stokes_dof_handler.end()),
 *                       worker,
 *                       copier,
 *                       Assembly::Scratch::StokesPreconditioner<dim>(
 *                         stokes_fe,
 *                         quadrature_formula,
 *                         mapping,
 *                         update_JxW_values | update_values | update_gradients),
 *                       Assembly::CopyData::StokesPreconditioner<dim>(stokes_fe));
 *   
 *       stokes_preconditioner_matrix.compress(VectorOperation::add);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The final function in this block initiates assembly of the Stokes
 * preconditioner matrix and then in fact builds the Stokes
 * preconditioner. It is mostly the same as in the serial case. The only
 * difference to @ref step_31 "step-31" is that we use a Jacobi preconditioner for the
 * pressure mass matrix instead of IC, as discussed in the introduction.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::build_stokes_preconditioner()
 *     {
 *       if (rebuild_stokes_preconditioner == false)
 *         return;
 *   
 *       TimerOutput::Scope timer_section(computing_timer,
 *                                        "   Build Stokes preconditioner");
 *       pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
 *   
 *       assemble_stokes_preconditioner();
 *   
 *       const FEValuesExtractors::Vector     velocity_components(0);
 *       const std::vector<std::vector<bool>> constant_modes =
 *         DoFTools::extract_constant_modes(
 *           stokes_dof_handler, stokes_fe.component_mask(velocity_components));
 *   
 *       Mp_preconditioner =
 *         std::make_shared<TrilinosWrappers::PreconditionJacobi>();
 *       Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionAMG>();
 *   
 *       TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
 *       Amg_data.constant_modes        = constant_modes;
 *       Amg_data.elliptic              = true;
 *       Amg_data.higher_order_elements = true;
 *       Amg_data.smoother_sweeps       = 2;
 *       Amg_data.aggregation_threshold = 0.02;
 *   
 *       Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1));
 *       Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0),
 *                                      Amg_data);
 *   
 *       rebuild_stokes_preconditioner = false;
 *   
 *       pcout << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-Stokessystemassembly"></a> 
 * <h5>Stokes system assembly</h5>
 * 
 
 * 
 * The next three functions implement the assembly of the Stokes system,
 * again split up into a part performing local calculations, one for writing
 * the local data into the global matrix and vector, and one for actually
 * running the loop over all cells with the help of the WorkStream
 * class. Note that the assembly of the Stokes matrix needs only to be done
 * in case we have changed the mesh. Otherwise, just the
 * (temperature-dependent) right hand side needs to be calculated
 * here. Since we are working with distributed matrices and vectors, we have
 * to call the respective <code>compress()</code> functions in the end of
 * the assembly in order to send non-local data to the owner process.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::local_assemble_stokes_system(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       Assembly::Scratch::StokesSystem<dim>                 &scratch,
 *       Assembly::CopyData::StokesSystem<dim>                &data)
 *     {
 *       const unsigned int dofs_per_cell =
 *         scratch.stokes_fe_values.get_fe().n_dofs_per_cell();
 *       const unsigned int n_q_points =
 *         scratch.stokes_fe_values.n_quadrature_points;
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *       const FEValuesExtractors::Scalar pressure(dim);
 *   
 *       scratch.stokes_fe_values.reinit(cell);
 *   
 *       const typename DoFHandler<dim>::active_cell_iterator temperature_cell =
 *         cell->as_dof_handler_iterator(temperature_dof_handler);
 *       scratch.temperature_fe_values.reinit(temperature_cell);
 *   
 *       if (rebuild_stokes_matrix)
 *         data.local_matrix = 0;
 *       data.local_rhs = 0;
 *   
 *       scratch.temperature_fe_values.get_function_values(
 *         old_temperature_solution, scratch.old_temperature_values);
 *   
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           const double old_temperature = scratch.old_temperature_values[q];
 *   
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value(k, q);
 *               if (rebuild_stokes_matrix)
 *                 {
 *                   scratch.grads_phi_u[k] =
 *                     scratch.stokes_fe_values[velocities].symmetric_gradient(k, q);
 *                   scratch.div_phi_u[k] =
 *                     scratch.stokes_fe_values[velocities].divergence(k, q);
 *                   scratch.phi_p[k] =
 *                     scratch.stokes_fe_values[pressure].value(k, q);
 *                 }
 *             }
 *   
 *           if (rebuild_stokes_matrix == true)
 *             for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 data.local_matrix(i, j) +=
 *                   (EquationData::eta * 2 *
 *                      (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]) -
 *                    (EquationData::pressure_scaling * scratch.div_phi_u[i] *
 *                     scratch.phi_p[j]) -
 *                    (EquationData::pressure_scaling * scratch.phi_p[i] *
 *                     scratch.div_phi_u[j])) *
 *                   scratch.stokes_fe_values.JxW(q);
 *   
 *           const Tensor<1, dim> gravity = EquationData::gravity_vector(
 *             scratch.stokes_fe_values.quadrature_point(q));
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             data.local_rhs(i) += (EquationData::density(old_temperature) *
 *                                   gravity * scratch.phi_u[i]) *
 *                                  scratch.stokes_fe_values.JxW(q);
 *         }
 *   
 *       cell->get_dof_indices(data.local_dof_indices);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::copy_local_to_global_stokes_system(
 *       const Assembly::CopyData::StokesSystem<dim> &data)
 *     {
 *       if (rebuild_stokes_matrix == true)
 *         stokes_constraints.distribute_local_to_global(data.local_matrix,
 *                                                       data.local_rhs,
 *                                                       data.local_dof_indices,
 *                                                       stokes_matrix,
 *                                                       stokes_rhs);
 *       else
 *         stokes_constraints.distribute_local_to_global(data.local_rhs,
 *                                                       data.local_dof_indices,
 *                                                       stokes_rhs);
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::assemble_stokes_system()
 *     {
 *       TimerOutput::Scope timer_section(computing_timer,
 *                                        "   Assemble Stokes system");
 *   
 *       if (rebuild_stokes_matrix == true)
 *         stokes_matrix = 0;
 *   
 *       stokes_rhs = 0;
 *   
 *       const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree + 1);
 *   
 *       using CellFilter =
 *         FilteredIterator<typename DoFHandler<2>::active_cell_iterator>;
 *   
 *       WorkStream::run(
 *         CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                    stokes_dof_handler.begin_active()),
 *         CellFilter(IteratorFilters::LocallyOwnedCell(), stokes_dof_handler.end()),
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                Assembly::Scratch::StokesSystem<dim>                 &scratch,
 *                Assembly::CopyData::StokesSystem<dim>                &data) {
 *           this->local_assemble_stokes_system(cell, scratch, data);
 *         },
 *         [this](const Assembly::CopyData::StokesSystem<dim> &data) {
 *           this->copy_local_to_global_stokes_system(data);
 *         },
 *         Assembly::Scratch::StokesSystem<dim>(
 *           stokes_fe,
 *           mapping,
 *           quadrature_formula,
 *           (update_values | update_quadrature_points | update_JxW_values |
 *            (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0))),
 *           temperature_fe,
 *           update_values),
 *         Assembly::CopyData::StokesSystem<dim>(stokes_fe));
 *   
 *       if (rebuild_stokes_matrix == true)
 *         stokes_matrix.compress(VectorOperation::add);
 *       stokes_rhs.compress(VectorOperation::add);
 *   
 *       rebuild_stokes_matrix = false;
 *   
 *       pcout << std::endl;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-Temperaturematrixassembly"></a> 
 * <h5>Temperature matrix assembly</h5>
 * 
 
 * 
 * The task to be performed by the next three functions is to calculate a
 * mass matrix and a Laplace matrix on the temperature system. These will be
 * combined in order to yield the semi-implicit time stepping matrix that
 * consists of the mass matrix plus a time step-dependent weight factor
 * times the Laplace matrix. This function is again essentially the body of
 * the loop over all cells from @ref step_31 "step-31".
 *   
 
 * 
 * The two following functions perform similar services as the ones above.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::local_assemble_temperature_matrix(
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       Assembly::Scratch::TemperatureMatrix<dim>            &scratch,
 *       Assembly::CopyData::TemperatureMatrix<dim>           &data)
 *     {
 *       const unsigned int dofs_per_cell =
 *         scratch.temperature_fe_values.get_fe().n_dofs_per_cell();
 *       const unsigned int n_q_points =
 *         scratch.temperature_fe_values.n_quadrature_points;
 *   
 *       scratch.temperature_fe_values.reinit(cell);
 *       cell->get_dof_indices(data.local_dof_indices);
 *   
 *       data.local_mass_matrix      = 0;
 *       data.local_stiffness_matrix = 0;
 *   
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               scratch.grad_phi_T[k] =
 *                 scratch.temperature_fe_values.shape_grad(k, q);
 *               scratch.phi_T[k] = scratch.temperature_fe_values.shape_value(k, q);
 *             }
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *               {
 *                 data.local_mass_matrix(i, j) +=
 *                   (scratch.phi_T[i] * scratch.phi_T[j] *
 *                    scratch.temperature_fe_values.JxW(q));
 *                 data.local_stiffness_matrix(i, j) +=
 *                   (EquationData::kappa * scratch.grad_phi_T[i] *
 *                    scratch.grad_phi_T[j] * scratch.temperature_fe_values.JxW(q));
 *               }
 *         }
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::copy_local_to_global_temperature_matrix(
 *       const Assembly::CopyData::TemperatureMatrix<dim> &data)
 *     {
 *       temperature_constraints.distribute_local_to_global(data.local_mass_matrix,
 *                                                          data.local_dof_indices,
 *                                                          temperature_mass_matrix);
 *       temperature_constraints.distribute_local_to_global(
 *         data.local_stiffness_matrix,
 *         data.local_dof_indices,
 *         temperature_stiffness_matrix);
 *     }
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::assemble_temperature_matrix()
 *     {
 *       if (rebuild_temperature_matrices == false)
 *         return;
 *   
 *       TimerOutput::Scope timer_section(computing_timer,
 *                                        "   Assemble temperature matrices");
 *       temperature_mass_matrix      = 0;
 *       temperature_stiffness_matrix = 0;
 *   
 *       const QGauss<dim> quadrature_formula(parameters.temperature_degree + 2);
 *   
 *       using CellFilter =
 *         FilteredIterator<typename DoFHandler<2>::active_cell_iterator>;
 *   
 *       WorkStream::run(
 *         CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                    temperature_dof_handler.begin_active()),
 *         CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                    temperature_dof_handler.end()),
 *         [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
 *                Assembly::Scratch::TemperatureMatrix<dim>            &scratch,
 *                Assembly::CopyData::TemperatureMatrix<dim>           &data) {
 *           this->local_assemble_temperature_matrix(cell, scratch, data);
 *         },
 *         [this](const Assembly::CopyData::TemperatureMatrix<dim> &data) {
 *           this->copy_local_to_global_temperature_matrix(data);
 *         },
 *         Assembly::Scratch::TemperatureMatrix<dim>(temperature_fe,
 *                                                   mapping,
 *                                                   quadrature_formula),
 *         Assembly::CopyData::TemperatureMatrix<dim>(temperature_fe));
 *   
 *       temperature_mass_matrix.compress(VectorOperation::add);
 *       temperature_stiffness_matrix.compress(VectorOperation::add);
 *   
 *       rebuild_temperature_matrices       = false;
 *       rebuild_temperature_preconditioner = true;
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-Temperaturerighthandsideassembly"></a> 
 * <h5>Temperature right hand side assembly</h5>
 * 
 
 * 
 * This is the last assembly function. It calculates the right hand side of
 * the temperature system, which includes the convection and the
 * stabilization terms. It includes a lot of evaluations of old solutions at
 * the quadrature points (which are necessary for calculating the artificial
 * viscosity of stabilization), but is otherwise similar to the other
 * assembly functions. Notice, once again, how we resolve the dilemma of
 * having inhomogeneous boundary conditions, by just making a right hand
 * side at this point (compare the comments for the <code>project()</code>
 * function above): We create some matrix columns with exactly the values
 * that would be entered for the temperature @ref GlossStiffnessMatrix "stiffness matrix", in case we
 * have inhomogeneously constrained dofs. That will account for the correct
 * balance of the right hand side vector with the matrix system of
 * temperature.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::local_assemble_temperature_rhs(
 *       const std::pair<double, double> global_T_range,
 *       const double                    global_max_velocity,
 *       const double                    global_entropy_variation,
 *       const typename DoFHandler<dim>::active_cell_iterator &cell,
 *       Assembly::Scratch::TemperatureRHS<dim>               &scratch,
 *       Assembly::CopyData::TemperatureRHS<dim>              &data)
 *     {
 *       const bool use_bdf2_scheme = (timestep_number != 0);
 *   
 *       const unsigned int dofs_per_cell =
 *         scratch.temperature_fe_values.get_fe().n_dofs_per_cell();
 *       const unsigned int n_q_points =
 *         scratch.temperature_fe_values.n_quadrature_points;
 *   
 *       const FEValuesExtractors::Vector velocities(0);
 *   
 *       data.local_rhs     = 0;
 *       data.matrix_for_bc = 0;
 *       cell->get_dof_indices(data.local_dof_indices);
 *   
 *       scratch.temperature_fe_values.reinit(cell);
 *   
 *       typename DoFHandler<dim>::active_cell_iterator stokes_cell =
 *         cell->as_dof_handler_iterator(stokes_dof_handler);
 *       scratch.stokes_fe_values.reinit(stokes_cell);
 *   
 *       scratch.temperature_fe_values.get_function_values(
 *         old_temperature_solution, scratch.old_temperature_values);
 *       scratch.temperature_fe_values.get_function_values(
 *         old_old_temperature_solution, scratch.old_old_temperature_values);
 *   
 *       scratch.temperature_fe_values.get_function_gradients(
 *         old_temperature_solution, scratch.old_temperature_grads);
 *       scratch.temperature_fe_values.get_function_gradients(
 *         old_old_temperature_solution, scratch.old_old_temperature_grads);
 *   
 *       scratch.temperature_fe_values.get_function_laplacians(
 *         old_temperature_solution, scratch.old_temperature_laplacians);
 *       scratch.temperature_fe_values.get_function_laplacians(
 *         old_old_temperature_solution, scratch.old_old_temperature_laplacians);
 *   
 *       scratch.stokes_fe_values[velocities].get_function_values(
 *         stokes_solution, scratch.old_velocity_values);
 *       scratch.stokes_fe_values[velocities].get_function_values(
 *         old_stokes_solution, scratch.old_old_velocity_values);
 *       scratch.stokes_fe_values[velocities].get_function_symmetric_gradients(
 *         stokes_solution, scratch.old_strain_rates);
 *       scratch.stokes_fe_values[velocities].get_function_symmetric_gradients(
 *         old_stokes_solution, scratch.old_old_strain_rates);
 *   
 *       const double nu =
 *         compute_viscosity(scratch.old_temperature_values,
 *                           scratch.old_old_temperature_values,
 *                           scratch.old_temperature_grads,
 *                           scratch.old_old_temperature_grads,
 *                           scratch.old_temperature_laplacians,
 *                           scratch.old_old_temperature_laplacians,
 *                           scratch.old_velocity_values,
 *                           scratch.old_old_velocity_values,
 *                           scratch.old_strain_rates,
 *                           scratch.old_old_strain_rates,
 *                           global_max_velocity,
 *                           global_T_range.second - global_T_range.first,
 *                           0.5 * (global_T_range.second + global_T_range.first),
 *                           global_entropy_variation,
 *                           cell->diameter());
 *   
 *       for (unsigned int q = 0; q < n_q_points; ++q)
 *         {
 *           for (unsigned int k = 0; k < dofs_per_cell; ++k)
 *             {
 *               scratch.phi_T[k] = scratch.temperature_fe_values.shape_value(k, q);
 *               scratch.grad_phi_T[k] =
 *                 scratch.temperature_fe_values.shape_grad(k, q);
 *             }
 *   
 *   
 *           const double T_term_for_rhs =
 *             (use_bdf2_scheme ?
 *                (scratch.old_temperature_values[q] *
 *                   (1 + time_step / old_time_step) -
 *                 scratch.old_old_temperature_values[q] * (time_step * time_step) /
 *                   (old_time_step * (time_step + old_time_step))) :
 *                scratch.old_temperature_values[q]);
 *   
 *           const double ext_T =
 *             (use_bdf2_scheme ? (scratch.old_temperature_values[q] *
 *                                   (1 + time_step / old_time_step) -
 *                                 scratch.old_old_temperature_values[q] *
 *                                   time_step / old_time_step) :
 *                                scratch.old_temperature_values[q]);
 *   
 *           const Tensor<1, dim> ext_grad_T =
 *             (use_bdf2_scheme ? (scratch.old_temperature_grads[q] *
 *                                   (1 + time_step / old_time_step) -
 *                                 scratch.old_old_temperature_grads[q] * time_step /
 *                                   old_time_step) :
 *                                scratch.old_temperature_grads[q]);
 *   
 *           const Tensor<1, dim> extrapolated_u =
 *             (use_bdf2_scheme ?
 *                (scratch.old_velocity_values[q] * (1 + time_step / old_time_step) -
 *                 scratch.old_old_velocity_values[q] * time_step / old_time_step) :
 *                scratch.old_velocity_values[q]);
 *   
 *           const SymmetricTensor<2, dim> extrapolated_strain_rate =
 *             (use_bdf2_scheme ?
 *                (scratch.old_strain_rates[q] * (1 + time_step / old_time_step) -
 *                 scratch.old_old_strain_rates[q] * time_step / old_time_step) :
 *                scratch.old_strain_rates[q]);
 *   
 *           const double gamma =
 *             ((EquationData::radiogenic_heating * EquationData::density(ext_T) +
 *               2 * EquationData::eta * extrapolated_strain_rate *
 *                 extrapolated_strain_rate) /
 *              (EquationData::density(ext_T) * EquationData::specific_heat));
 *   
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               data.local_rhs(i) +=
 *                 (T_term_for_rhs * scratch.phi_T[i] -
 *                  time_step * extrapolated_u * ext_grad_T * scratch.phi_T[i] -
 *                  time_step * nu * ext_grad_T * scratch.grad_phi_T[i] +
 *                  time_step * gamma * scratch.phi_T[i]) *
 *                 scratch.temperature_fe_values.JxW(q);
 *   
 *               if (temperature_constraints.is_inhomogeneously_constrained(
 *                     data.local_dof_indices[i]))
 *                 {
 *                   for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                     data.matrix_for_bc(j, i) +=
 *                       (scratch.phi_T[i] * scratch.phi_T[j] *
 *                          (use_bdf2_scheme ? ((2 * time_step + old_time_step) /
 *                                              (time_step + old_time_step)) :
 *                                             1.) +
 *                        scratch.grad_phi_T[i] * scratch.grad_phi_T[j] *
 *                          EquationData::kappa * time_step) *
 *                       scratch.temperature_fe_values.JxW(q);
 *                 }
 *             }
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::copy_local_to_global_temperature_rhs(
 *       const Assembly::CopyData::TemperatureRHS<dim> &data)
 *     {
 *       temperature_constraints.distribute_local_to_global(data.local_rhs,
 *                                                          data.local_dof_indices,
 *                                                          temperature_rhs,
 *                                                          data.matrix_for_bc);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * In the function that runs the WorkStream for actually calculating the
 * right hand side, we also generate the final matrix. As mentioned above,
 * it is a sum of the mass matrix and the Laplace matrix, times some time
 * step-dependent weight. This weight is specified by the BDF-2 time
 * integration scheme, see the introduction in @ref step_31 "step-31". What is new in this
 * tutorial program (in addition to the use of MPI parallelization and the
 * WorkStream class), is that we now precompute the temperature
 * preconditioner as well. The reason is that the setup of the Jacobi
 * preconditioner takes a noticeable time compared to the solver because we
 * usually only need between 10 and 20 iterations for solving the
 * temperature system (this might sound strange, as Jacobi really only
 * consists of a diagonal, but in Trilinos it is derived from more general
 * framework for point relaxation preconditioners which is a bit
 * inefficient). Hence, it is more efficient to precompute the
 * preconditioner, even though the matrix entries may slightly change
 * because the time step might change. This is not too big a problem because
 * we remesh every few time steps (and regenerate the preconditioner then).
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::assemble_temperature_system(
 *       const double maximal_velocity)
 *     {
 *       const bool use_bdf2_scheme = (timestep_number != 0);
 *   
 *       if (use_bdf2_scheme == true)
 *         {
 *           temperature_matrix.copy_from(temperature_mass_matrix);
 *           temperature_matrix *=
 *             (2 * time_step + old_time_step) / (time_step + old_time_step);
 *           temperature_matrix.add(time_step, temperature_stiffness_matrix);
 *         }
 *       else
 *         {
 *           temperature_matrix.copy_from(temperature_mass_matrix);
 *           temperature_matrix.add(time_step, temperature_stiffness_matrix);
 *         }
 *   
 *       if (rebuild_temperature_preconditioner == true)
 *         {
 *           T_preconditioner =
 *             std::make_shared<TrilinosWrappers::PreconditionJacobi>();
 *           T_preconditioner->initialize(temperature_matrix);
 *           rebuild_temperature_preconditioner = false;
 *         }
 *   
 * @endcode
 * 
 * The next part is computing the right hand side vectors.  To do so, we
 * first compute the average temperature @f$T_m@f$ that we use for evaluating
 * the artificial viscosity stabilization through the residual @f$E(T) =
 * (T-T_m)^2@f$. We do this by defining the midpoint between maximum and
 * minimum temperature as average temperature in the definition of the
 * entropy viscosity. An alternative would be to use the integral average,
 * but the results are not very sensitive to this choice. The rest then
 * only requires calling WorkStream::run again, binding the arguments to
 * the <code>local_assemble_temperature_rhs</code> function that are the
 * same in every call to the correct values:
 * 
 * @code
 *       temperature_rhs = 0;
 *   
 *       const QGauss<dim> quadrature_formula(parameters.temperature_degree + 2);
 *       const std::pair<double, double> global_T_range =
 *         get_extrapolated_temperature_range();
 *   
 *       const double average_temperature =
 *         0.5 * (global_T_range.first + global_T_range.second);
 *       const double global_entropy_variation =
 *         get_entropy_variation(average_temperature);
 *   
 *       using CellFilter =
 *         FilteredIterator<typename DoFHandler<2>::active_cell_iterator>;
 *   
 *       auto worker =
 *         [this, global_T_range, maximal_velocity, global_entropy_variation](
 *           const typename DoFHandler<dim>::active_cell_iterator &cell,
 *           Assembly::Scratch::TemperatureRHS<dim>               &scratch,
 *           Assembly::CopyData::TemperatureRHS<dim>              &data) {
 *           this->local_assemble_temperature_rhs(global_T_range,
 *                                                maximal_velocity,
 *                                                global_entropy_variation,
 *                                                cell,
 *                                                scratch,
 *                                                data);
 *         };
 *   
 *       auto copier = [this](const Assembly::CopyData::TemperatureRHS<dim> &data) {
 *         this->copy_local_to_global_temperature_rhs(data);
 *       };
 *   
 *       WorkStream::run(CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                                  temperature_dof_handler.begin_active()),
 *                       CellFilter(IteratorFilters::LocallyOwnedCell(),
 *                                  temperature_dof_handler.end()),
 *                       worker,
 *                       copier,
 *                       Assembly::Scratch::TemperatureRHS<dim>(
 *                         temperature_fe, stokes_fe, mapping, quadrature_formula),
 *                       Assembly::CopyData::TemperatureRHS<dim>(temperature_fe));
 *   
 *       temperature_rhs.compress(VectorOperation::add);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemsolve"></a> 
 * <h4>BoussinesqFlowProblem::solve</h4>
 * 
 
 * 
 * This function solves the linear systems in each time step of the
 * Boussinesq problem. First, we work on the Stokes system and then on the
 * temperature system. In essence, it does the same things as the respective
 * function in @ref step_31 "step-31". However, there are a few changes here.
 *   
 
 * 
 * The first change is related to the way we store our solution: we keep the
 * vectors with locally owned degrees of freedom plus ghost nodes on each
 * MPI node. When we enter a solver which is supposed to perform
 * matrix-vector products with a distributed matrix, this is not the
 * appropriate form, though. There, we will want to have the solution vector
 * to be distributed in the same way as the matrix, i.e. without any
 * ghosts. So what we do first is to generate a distributed vector called
 * <code>distributed_stokes_solution</code> and put only the locally owned
 * dofs into that, which is neatly done by the <code>operator=</code> of the
 * Trilinos vector.
 *   
 
 * 
 * Next, we scale the pressure solution (or rather, the initial guess) for
 * the solver so that it matches with the length scales in the matrices, as
 * discussed in the introduction. We also immediately scale the pressure
 * solution back to the correct units after the solution is completed.  We
 * also need to set the pressure values at hanging nodes to zero. This we
 * also did in @ref step_31 "step-31" in order not to disturb the Schur complement by some
 * vector entries that actually are irrelevant during the solve stage. As a
 * difference to @ref step_31 "step-31", here we do it only for the locally owned pressure
 * dofs. After solving for the Stokes solution, each processor copies the
 * distributed solution back into the solution vector that also includes
 * ghost elements.
 *   
 
 * 
 * The third and most obvious change is that we have two variants for the
 * Stokes solver: A fast solver that sometimes breaks down, and a robust
 * solver that is slower. This is what we already discussed in the
 * introduction. Here is how we realize it: First, we perform 30 iterations
 * with the fast solver based on the simple preconditioner based on the AMG
 * V-cycle instead of an approximate solve (this is indicated by the
 * <code>false</code> argument to the
 * <code>LinearSolvers::BlockSchurPreconditioner</code> object). If we
 * converge, everything is fine. If we do not converge, the solver control
 * object will throw an exception SolverControl::NoConvergence. Usually,
 * this would abort the program because we don't catch them in our usual
 * <code>solve()</code> functions. This is certainly not what we want to
 * happen here. Rather, we want to switch to the strong solver and continue
 * the solution process with whatever vector we got so far. Hence, we catch
 * the exception with the C++ try/catch mechanism. We then simply go through
 * the same solver sequence again in the <code>catch</code> clause, this
 * time passing the @p true flag to the preconditioner for the strong
 * solver, signaling an approximate CG solve.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::solve()
 *     {
 *       {
 *         TimerOutput::Scope timer_section(computing_timer,
 *                                          "   Solve Stokes system");
 *   
 *         pcout << "   Solving Stokes system... " << std::flush;
 *   
 *         TrilinosWrappers::MPI::BlockVector distributed_stokes_solution(
 *           stokes_rhs);
 *         distributed_stokes_solution = stokes_solution;
 *   
 *         distributed_stokes_solution.block(1) /= EquationData::pressure_scaling;
 *   
 *         const unsigned int
 *           start = (distributed_stokes_solution.block(0).size() +
 *                    distributed_stokes_solution.block(1).local_range().first),
 *           end   = (distributed_stokes_solution.block(0).size() +
 *                  distributed_stokes_solution.block(1).local_range().second);
 *         for (unsigned int i = start; i < end; ++i)
 *           if (stokes_constraints.is_constrained(i))
 *             distributed_stokes_solution(i) = 0;
 *   
 *   
 *         PrimitiveVectorMemory<TrilinosWrappers::MPI::BlockVector> mem;
 *   
 *         unsigned int  n_iterations     = 0;
 *         const double  solver_tolerance = 1e-8 * stokes_rhs.l2_norm();
 *         SolverControl solver_control(30, solver_tolerance);
 *   
 *         try
 *           {
 *             const LinearSolvers::BlockSchurPreconditioner<
 *               TrilinosWrappers::PreconditionAMG,
 *               TrilinosWrappers::PreconditionJacobi>
 *               preconditioner(stokes_matrix,
 *                              stokes_preconditioner_matrix,
 *                              *Mp_preconditioner,
 *                              *Amg_preconditioner,
 *                              false);
 *   
 *             SolverFGMRES<TrilinosWrappers::MPI::BlockVector> solver(
 *               solver_control,
 *               mem,
 *               SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(
 *                 30));
 *             solver.solve(stokes_matrix,
 *                          distributed_stokes_solution,
 *                          stokes_rhs,
 *                          preconditioner);
 *   
 *             n_iterations = solver_control.last_step();
 *           }
 *   
 *         catch (SolverControl::NoConvergence &)
 *           {
 *             const LinearSolvers::BlockSchurPreconditioner<
 *               TrilinosWrappers::PreconditionAMG,
 *               TrilinosWrappers::PreconditionJacobi>
 *               preconditioner(stokes_matrix,
 *                              stokes_preconditioner_matrix,
 *                              *Mp_preconditioner,
 *                              *Amg_preconditioner,
 *                              true);
 *   
 *             SolverControl solver_control_refined(stokes_matrix.m(),
 *                                                  solver_tolerance);
 *             SolverFGMRES<TrilinosWrappers::MPI::BlockVector> solver(
 *               solver_control_refined,
 *               mem,
 *               SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(
 *                 50));
 *             solver.solve(stokes_matrix,
 *                          distributed_stokes_solution,
 *                          stokes_rhs,
 *                          preconditioner);
 *   
 *             n_iterations =
 *               (solver_control.last_step() + solver_control_refined.last_step());
 *           }
 *   
 *   
 *         stokes_constraints.distribute(distributed_stokes_solution);
 *   
 *         distributed_stokes_solution.block(1) *= EquationData::pressure_scaling;
 *   
 *         stokes_solution = distributed_stokes_solution;
 *         pcout << n_iterations << " iterations." << std::endl;
 *       }
 *   
 *   
 * @endcode
 * 
 * Now let's turn to the temperature part: First, we compute the time step
 * size. We found that we need smaller time steps for 3d than for 2d for
 * the shell geometry. This is because the cells are more distorted in
 * that case (it is the smallest edge length that determines the CFL
 * number). Instead of computing the time step from maximum velocity and
 * minimal mesh size as in @ref step_31 "step-31", we compute local CFL numbers, i.e., on
 * each cell we compute the maximum velocity times the mesh size, and
 * compute the maximum of them. Hence, we need to choose the factor in
 * front of the time step slightly smaller. (We later re-considered this
 * approach towards time stepping. If you're curious about this, you may
 * want to read the time stepping section in @cite HDGB17 .)
 *     
 
 * 
 * After temperature right hand side assembly, we solve the linear
 * system for temperature (with fully distributed vectors without
 * ghost elements and using the solution from the last timestep as
 * our initial guess for the iterative solver), apply constraints,
 * and copy the vector back to one with ghosts.
 *     
 
 * 
 * In the end, we extract the temperature range similarly to @ref step_31 "step-31" to
 * produce some output (for example in order to help us choose the
 * stabilization constants, as discussed in the introduction). The only
 * difference is that we need to exchange maxima over all processors.
 * 
 * @code
 *       {
 *         TimerOutput::Scope timer_section(computing_timer,
 *                                          "   Assemble temperature rhs");
 *   
 *         old_time_step = time_step;
 *   
 *         const double scaling = (dim == 3 ? 0.25 : 1.0);
 *         time_step            = (scaling / (2.1 * dim * std::sqrt(1. * dim)) /
 *                      (parameters.temperature_degree * get_cfl_number()));
 *   
 *         const double maximal_velocity = get_maximal_velocity();
 *         pcout << "   Maximal velocity: "
 *               << maximal_velocity * EquationData::year_in_seconds * 100
 *               << " cm/year" << std::endl;
 *         pcout << "   "
 *               << "Time step: " << time_step / EquationData::year_in_seconds
 *               << " years" << std::endl;
 *   
 *         assemble_temperature_system(maximal_velocity);
 *       }
 *   
 *       {
 *         TimerOutput::Scope timer_section(computing_timer,
 *                                          "   Solve temperature system");
 *   
 *         SolverControl solver_control(temperature_matrix.m(),
 *                                      1e-12 * temperature_rhs.l2_norm());
 *         SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
 *   
 *         TrilinosWrappers::MPI::Vector distributed_temperature_solution(
 *           temperature_rhs);
 *         distributed_temperature_solution = old_temperature_solution;
 *   
 *         cg.solve(temperature_matrix,
 *                  distributed_temperature_solution,
 *                  temperature_rhs,
 *                  *T_preconditioner);
 *   
 *         temperature_constraints.distribute(distributed_temperature_solution);
 *         temperature_solution = distributed_temperature_solution;
 *   
 *         pcout << "   " << solver_control.last_step()
 *               << " CG iterations for temperature" << std::endl;
 *   
 *         double temperature[2] = {std::numeric_limits<double>::max(),
 *                                  -std::numeric_limits<double>::max()};
 *         double global_temperature[2];
 *   
 *         for (unsigned int i =
 *                distributed_temperature_solution.local_range().first;
 *              i < distributed_temperature_solution.local_range().second;
 *              ++i)
 *           {
 *             temperature[0] =
 *               std::min<double>(temperature[0],
 *                                distributed_temperature_solution(i));
 *             temperature[1] =
 *               std::max<double>(temperature[1],
 *                                distributed_temperature_solution(i));
 *           }
 *   
 *         temperature[0] *= -1.0;
 *         Utilities::MPI::max(temperature, MPI_COMM_WORLD, global_temperature);
 *         global_temperature[0] *= -1.0;
 *   
 *         pcout << "   Temperature range: " << global_temperature[0] << ' '
 *               << global_temperature[1] << std::endl;
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemoutput_results"></a> 
 * <h4>BoussinesqFlowProblem::output_results</h4>
 * 
 
 * 
 * Next comes the function that generates the output. The quantities to
 * output could be introduced manually like we did in @ref step_31 "step-31". An
 * alternative is to hand this task over to a class PostProcessor that
 * inherits from the class DataPostprocessor, which can be attached to
 * DataOut. This allows us to output derived quantities from the solution,
 * like the friction heating included in this example. It overloads the
 * virtual function DataPostprocessor::evaluate_vector_field(),
 * which is then internally called from DataOut::build_patches(). We have to
 * give it values of the numerical solution, its derivatives, normals to the
 * cell, the actual evaluation points and any additional quantities. This
 * follows the same procedure as discussed in @ref step_29 "step-29" and other programs.
 * 
 * @code
 *     template <int dim>
 *     class BoussinesqFlowProblem<dim>::Postprocessor
 *       : public DataPostprocessor<dim>
 *     {
 *     public:
 *       Postprocessor(const unsigned int partition, const double minimal_pressure);
 *   
 *       virtual void evaluate_vector_field(
 *         const DataPostprocessorInputs::Vector<dim> &inputs,
 *         std::vector<Vector<double>> &computed_quantities) const override;
 *   
 *       virtual std::vector<std::string> get_names() const override;
 *   
 *       virtual std::vector<
 *         DataComponentInterpretation::DataComponentInterpretation>
 *       get_data_component_interpretation() const override;
 *   
 *       virtual UpdateFlags get_needed_update_flags() const override;
 *   
 *     private:
 *       const unsigned int partition;
 *       const double       minimal_pressure;
 *     };
 *   
 *   
 *     template <int dim>
 *     BoussinesqFlowProblem<dim>::Postprocessor::Postprocessor(
 *       const unsigned int partition,
 *       const double       minimal_pressure)
 *       : partition(partition)
 *       , minimal_pressure(minimal_pressure)
 *     {}
 *   
 *   
 * @endcode
 * 
 * Here we define the names for the variables we want to output. These are
 * the actual solution values for velocity, pressure, and temperature, as
 * well as the friction heating and to each cell the number of the processor
 * that owns it. This allows us to visualize the partitioning of the domain
 * among the processors. Except for the velocity, which is vector-valued,
 * all other quantities are scalar.
 * 
 * @code
 *     template <int dim>
 *     std::vector<std::string>
 *     BoussinesqFlowProblem<dim>::Postprocessor::get_names() const
 *     {
 *       std::vector<std::string> solution_names(dim, "velocity");
 *       solution_names.emplace_back("p");
 *       solution_names.emplace_back("T");
 *       solution_names.emplace_back("friction_heating");
 *       solution_names.emplace_back("partition");
 *   
 *       return solution_names;
 *     }
 *   
 *   
 *     template <int dim>
 *     std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *     BoussinesqFlowProblem<dim>::Postprocessor::get_data_component_interpretation()
 *       const
 *     {
 *       std::vector<DataComponentInterpretation::DataComponentInterpretation>
 *         interpretation(dim,
 *                        DataComponentInterpretation::component_is_part_of_vector);
 *   
 *       interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 *       interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 *       interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 *       interpretation.push_back(DataComponentInterpretation::component_is_scalar);
 *   
 *       return interpretation;
 *     }
 *   
 *   
 *     template <int dim>
 *     UpdateFlags
 *     BoussinesqFlowProblem<dim>::Postprocessor::get_needed_update_flags() const
 *     {
 *       return update_values | update_gradients | update_quadrature_points;
 *     }
 *   
 *   
 * @endcode
 * 
 * Now we implement the function that computes the derived quantities. As we
 * also did for the output, we rescale the velocity from its SI units to
 * something more readable, namely cm/year. Next, the pressure is scaled to
 * be between 0 and the maximum pressure. This makes it more easily
 * comparable -- in essence making all pressure variables positive or
 * zero. Temperature is taken as is, and the friction heating is computed as
 * @f$2 \eta \varepsilon(\mathbf{u}) \cdot \varepsilon(\mathbf{u})@f$.
 *   
 
 * 
 * The quantities we output here are more for illustration, rather than for
 * actual scientific value. We come back to this briefly in the results
 * section of this program and explain what one may in fact be interested in.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::Postprocessor::evaluate_vector_field(
 *       const DataPostprocessorInputs::Vector<dim> &inputs,
 *       std::vector<Vector<double>>                &computed_quantities) const
 *     {
 *       const unsigned int n_evaluation_points = inputs.solution_values.size();
 *       Assert(inputs.solution_gradients.size() == n_evaluation_points,
 *              ExcInternalError());
 *       Assert(computed_quantities.size() == n_evaluation_points,
 *              ExcInternalError());
 *       Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError());
 *   
 *       for (unsigned int p = 0; p < n_evaluation_points; ++p)
 *         {
 *           for (unsigned int d = 0; d < dim; ++d)
 *             computed_quantities[p](d) = (inputs.solution_values[p](d) *
 *                                          EquationData::year_in_seconds * 100);
 *   
 *           const double pressure =
 *             (inputs.solution_values[p](dim) - minimal_pressure);
 *           computed_quantities[p](dim) = pressure;
 *   
 *           const double temperature        = inputs.solution_values[p](dim + 1);
 *           computed_quantities[p](dim + 1) = temperature;
 *   
 *           Tensor<2, dim> grad_u;
 *           for (unsigned int d = 0; d < dim; ++d)
 *             grad_u[d] = inputs.solution_gradients[p][d];
 *           const SymmetricTensor<2, dim> strain_rate = symmetrize(grad_u);
 *           computed_quantities[p](dim + 2) =
 *             2 * EquationData::eta * strain_rate * strain_rate;
 *   
 *           computed_quantities[p](dim + 3) = partition;
 *         }
 *     }
 *   
 *   
 * @endcode
 * 
 * The <code>output_results()</code> function has a similar task to the one
 * in @ref step_31 "step-31". However, here we are going to demonstrate a different
 * technique on how to merge output from different DoFHandler objects. The
 * way we're going to achieve this recombination is to create a joint
 * DoFHandler that collects both components, the Stokes solution and the
 * temperature solution. This can be nicely done by combining the finite
 * elements from the two systems to form one FESystem, and let this
 * collective system define a new DoFHandler object. To be sure that
 * everything was done correctly, we perform a sanity check that ensures
 * that we got all the dofs from both Stokes and temperature even in the
 * combined system. We then combine the data vectors. Unfortunately, there
 * is no straight-forward relation that tells us how to sort Stokes and
 * temperature vector into the joint vector. The way we can get around this
 * trouble is to rely on the information collected in the FESystem. For each
 * dof on a cell, the joint finite element knows to which equation component
 * (velocity component, pressure, or temperature) it belongs – that's the
 * information we need! So we step through all cells (with iterators into
 * all three DoFHandlers moving in sync), and for each joint cell dof, we
 * read out that component using the FiniteElement::system_to_base_index
 * function (see there for a description of what the various parts of its
 * return value contain). We also need to keep track whether we're on a
 * Stokes dof or a temperature dof, which is contained in
 * joint_fe.system_to_base_index(i).first.first. Eventually, the dof_indices
 * data structures on either of the three systems tell us how the relation
 * between global vector and local dofs looks like on the present cell,
 * which concludes this tedious work. We make sure that each processor only
 * works on the subdomain it owns locally (and not on ghost or artificial
 * cells) when building the joint solution vector. The same will then have
 * to be done in DataOut::build_patches(), but that function does so
 * automatically.
 *   
 
 * 
 * What we end up with is a set of patches that we can write using the
 * functions in DataOutBase in a variety of output formats. Here, we then
 * have to pay attention that what each processor writes is really only its
 * own part of the domain, i.e. we will want to write each processor's
 * contribution into a separate file. This we do by adding an additional
 * number to the filename when we write the solution. This is not really
 * new, we did it similarly in @ref step_40 "step-40". Note that we write in the compressed
 * format @p .vtu instead of plain vtk files, which saves quite some
 * storage.
 *   
 
 * 
 * All the rest of the work is done in the PostProcessor class.
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::output_results()
 *     {
 *       TimerOutput::Scope timer_section(computing_timer, "Postprocessing");
 *   
 *       const FESystem<dim> joint_fe(stokes_fe, 1, temperature_fe, 1);
 *   
 *       DoFHandler<dim> joint_dof_handler(triangulation);
 *       joint_dof_handler.distribute_dofs(joint_fe);
 *       Assert(joint_dof_handler.n_dofs() ==
 *                stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(),
 *              ExcInternalError());
 *   
 *       TrilinosWrappers::MPI::Vector joint_solution;
 *       joint_solution.reinit(joint_dof_handler.locally_owned_dofs(),
 *                             MPI_COMM_WORLD);
 *   
 *       {
 *         std::vector<types::global_dof_index> local_joint_dof_indices(
 *           joint_fe.n_dofs_per_cell());
 *         std::vector<types::global_dof_index> local_stokes_dof_indices(
 *           stokes_fe.n_dofs_per_cell());
 *         std::vector<types::global_dof_index> local_temperature_dof_indices(
 *           temperature_fe.n_dofs_per_cell());
 *   
 *         typename DoFHandler<dim>::active_cell_iterator
 *           joint_cell       = joint_dof_handler.begin_active(),
 *           joint_endc       = joint_dof_handler.end(),
 *           stokes_cell      = stokes_dof_handler.begin_active(),
 *           temperature_cell = temperature_dof_handler.begin_active();
 *         for (; joint_cell != joint_endc;
 *              ++joint_cell, ++stokes_cell, ++temperature_cell)
 *           if (joint_cell->is_locally_owned())
 *             {
 *               joint_cell->get_dof_indices(local_joint_dof_indices);
 *               stokes_cell->get_dof_indices(local_stokes_dof_indices);
 *               temperature_cell->get_dof_indices(local_temperature_dof_indices);
 *   
 *               for (unsigned int i = 0; i < joint_fe.n_dofs_per_cell(); ++i)
 *                 if (joint_fe.system_to_base_index(i).first.first == 0)
 *                   {
 *                     Assert(joint_fe.system_to_base_index(i).second <
 *                              local_stokes_dof_indices.size(),
 *                            ExcInternalError());
 *   
 *                     joint_solution(local_joint_dof_indices[i]) = stokes_solution(
 *                       local_stokes_dof_indices[joint_fe.system_to_base_index(i)
 *                                                  .second]);
 *                   }
 *                 else
 *                   {
 *                     Assert(joint_fe.system_to_base_index(i).first.first == 1,
 *                            ExcInternalError());
 *                     Assert(joint_fe.system_to_base_index(i).second <
 *                              local_temperature_dof_indices.size(),
 *                            ExcInternalError());
 *                     joint_solution(local_joint_dof_indices[i]) =
 *                       temperature_solution(
 *                         local_temperature_dof_indices
 *                           [joint_fe.system_to_base_index(i).second]);
 *                   }
 *             }
 *       }
 *   
 *       joint_solution.compress(VectorOperation::insert);
 *   
 *       const IndexSet locally_relevant_joint_dofs =
 *         DoFTools::extract_locally_relevant_dofs(joint_dof_handler);
 *       TrilinosWrappers::MPI::Vector locally_relevant_joint_solution;
 *       locally_relevant_joint_solution.reinit(locally_relevant_joint_dofs,
 *                                              MPI_COMM_WORLD);
 *       locally_relevant_joint_solution = joint_solution;
 *   
 *       Postprocessor postprocessor(Utilities::MPI::this_mpi_process(
 *                                     MPI_COMM_WORLD),
 *                                   stokes_solution.block(1).min());
 *   
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(joint_dof_handler);
 *       data_out.add_data_vector(locally_relevant_joint_solution, postprocessor);
 *       data_out.build_patches();
 *   
 *       static int out_index = 0;
 *       data_out.write_vtu_with_pvtu_record(
 *         "./", "solution", out_index, MPI_COMM_WORLD, 5);
 *   
 *       ++out_index;
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemrefine_mesh"></a> 
 * <h4>BoussinesqFlowProblem::refine_mesh</h4>
 * 
 
 * 
 * This function isn't really new either. Since the <code>setup_dofs</code>
 * function that we call in the middle has its own timer section, we split
 * timing this function into two sections. It will also allow us to easily
 * identify which of the two is more expensive.
 *   
 
 * 
 * One thing of note, however, is that we only want to compute error
 * indicators on the locally owned subdomain. In order to achieve this, we
 * pass one additional argument to the KellyErrorEstimator::estimate
 * function. Note that the vector for error estimates is resized to the
 * number of active cells present on the current process, which is less than
 * the total number of active cells on all processors (but more than the
 * number of locally owned active cells); each processor only has a few
 * coarse cells around the locally owned ones, as also explained in @ref step_40 "step-40".
 *   
 
 * 
 * The local error estimates are then handed to a %parallel version of
 * GridRefinement (in namespace parallel::distributed::GridRefinement, see
 * also @ref step_40 "step-40") which looks at the errors and finds the cells that need
 * refinement by comparing the error values across processors. As in
 * @ref step_31 "step-31", we want to limit the maximum grid level. So in case some cells
 * have been marked that are already at the finest level, we simply clear
 * the refine flags.
 * 
 * @code
 *     template <int dim>
 *     void
 *     BoussinesqFlowProblem<dim>::refine_mesh(const unsigned int max_grid_level)
 *     {
 *       parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector>
 *         temperature_trans(temperature_dof_handler);
 *       parallel::distributed::SolutionTransfer<dim,
 *                                               TrilinosWrappers::MPI::BlockVector>
 *         stokes_trans(stokes_dof_handler);
 *   
 *       {
 *         TimerOutput::Scope timer_section(computing_timer,
 *                                          "Refine mesh structure, part 1");
 *   
 *         Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *         KellyErrorEstimator<dim>::estimate(
 *           temperature_dof_handler,
 *           QGauss<dim - 1>(parameters.temperature_degree + 1),
 *           std::map<types::boundary_id, const Function<dim> *>(),
 *           temperature_solution,
 *           estimated_error_per_cell,
 *           ComponentMask(),
 *           nullptr,
 *           0,
 *           triangulation.locally_owned_subdomain());
 *   
 *         parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction(
 *           triangulation, estimated_error_per_cell, 0.3, 0.1);
 *   
 *         if (triangulation.n_levels() > max_grid_level)
 *           for (typename Triangulation<dim>::active_cell_iterator cell =
 *                  triangulation.begin_active(max_grid_level);
 *                cell != triangulation.end();
 *                ++cell)
 *             cell->clear_refine_flag();
 *   
 * @endcode
 * 
 * With all flags marked as necessary, we can then tell the
 * parallel::distributed::SolutionTransfer objects to get ready to
 * transfer data from one mesh to the next, which they will do when
 * notified by
 * Triangulation as part of the @p execute_coarsening_and_refinement() call.
 * The syntax is similar to the non-%parallel solution transfer (with the
 * exception that here a pointer to the vector entries is enough). The
 * remainder of the function further down below is then concerned with
 * setting up the data structures again after mesh refinement and
 * restoring the solution vectors on the new mesh.
 * 
 * @code
 *         const std::vector<const TrilinosWrappers::MPI::Vector *> x_temperature = {
 *           &temperature_solution, &old_temperature_solution};
 *         const std::vector<const TrilinosWrappers::MPI::BlockVector *> x_stokes = {
 *           &stokes_solution, &old_stokes_solution};
 *   
 *         triangulation.prepare_coarsening_and_refinement();
 *   
 *         temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
 *         stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
 *   
 *         triangulation.execute_coarsening_and_refinement();
 *       }
 *   
 *       setup_dofs();
 *   
 *       {
 *         TimerOutput::Scope timer_section(computing_timer,
 *                                          "Refine mesh structure, part 2");
 *   
 *         {
 *           TrilinosWrappers::MPI::Vector distributed_temp1(temperature_rhs);
 *           TrilinosWrappers::MPI::Vector distributed_temp2(temperature_rhs);
 *   
 *           std::vector<TrilinosWrappers::MPI::Vector *> tmp = {&distributed_temp1,
 *                                                               &distributed_temp2};
 *           temperature_trans.interpolate(tmp);
 *   
 * @endcode
 * 
 * enforce constraints to make the interpolated solution conforming on
 * the new mesh:
 * 
 * @code
 *           temperature_constraints.distribute(distributed_temp1);
 *           temperature_constraints.distribute(distributed_temp2);
 *   
 *           temperature_solution     = distributed_temp1;
 *           old_temperature_solution = distributed_temp2;
 *         }
 *   
 *         {
 *           TrilinosWrappers::MPI::BlockVector distributed_stokes(stokes_rhs);
 *           TrilinosWrappers::MPI::BlockVector old_distributed_stokes(stokes_rhs);
 *   
 *           std::vector<TrilinosWrappers::MPI::BlockVector *> stokes_tmp = {
 *             &distributed_stokes, &old_distributed_stokes};
 *   
 *           stokes_trans.interpolate(stokes_tmp);
 *   
 * @endcode
 * 
 * enforce constraints to make the interpolated solution conforming on
 * the new mesh:
 * 
 * @code
 *           stokes_constraints.distribute(distributed_stokes);
 *           stokes_constraints.distribute(old_distributed_stokes);
 *   
 *           stokes_solution     = distributed_stokes;
 *           old_stokes_solution = old_distributed_stokes;
 *         }
 *       }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-BoussinesqFlowProblemrun"></a> 
 * <h4>BoussinesqFlowProblem::run</h4>
 * 
 
 * 
 * This is the final and controlling function in this class. It, in fact,
 * runs the entire rest of the program and is, once more, very similar to
 * @ref step_31 "step-31". The only substantial difference is that we use a different mesh
 * now (a GridGenerator::hyper_shell instead of a simple cube geometry).
 * 
 * @code
 *     template <int dim>
 *     void BoussinesqFlowProblem<dim>::run()
 *     {
 *       GridGenerator::hyper_shell(triangulation,
 *                                  Point<dim>(),
 *                                  EquationData::R0,
 *                                  EquationData::R1,
 *                                  (dim == 3) ? 96 : 12,
 *                                  true);
 *   
 *       global_Omega_diameter = GridTools::diameter(triangulation);
 *   
 *       triangulation.refine_global(parameters.initial_global_refinement);
 *   
 *       setup_dofs();
 *   
 *       unsigned int pre_refinement_step = 0;
 *   
 *     start_time_iteration:
 *   
 *       {
 *         TrilinosWrappers::MPI::Vector solution(
 *           temperature_dof_handler.locally_owned_dofs());
 * @endcode
 * 
 * VectorTools::project supports parallel vector classes with most
 * standard finite elements via deal.II's own native MatrixFree framework:
 * since we use standard Lagrange elements of moderate order this function
 * works well here.
 * 
 * @code
 *         VectorTools::project(temperature_dof_handler,
 *                              temperature_constraints,
 *                              QGauss<dim>(parameters.temperature_degree + 2),
 *                              EquationData::TemperatureInitialValues<dim>(),
 *                              solution);
 * @endcode
 * 
 * Having so computed the current temperature field, let us set the member
 * variable that holds the temperature nodes. Strictly speaking, we really
 * only need to set <code>old_temperature_solution</code> since the first
 * thing we will do is to compute the Stokes solution that only requires
 * the previous time step's temperature field. That said, nothing good can
 * come from not initializing the other vectors as well (especially since
 * it's a relatively cheap operation and we only have to do it once at the
 * beginning of the program) if we ever want to extend our numerical
 * method or physical model, and so we initialize
 * <code>old_temperature_solution</code> and
 * <code>old_old_temperature_solution</code> as well. The assignment makes
 * sure that the vectors on the left hand side (which where initialized to
 * contain ghost elements as well) also get the correct ghost elements. In
 * other words, the assignment here requires communication between
 * processors:
 * 
 * @code
 *         temperature_solution         = solution;
 *         old_temperature_solution     = solution;
 *         old_old_temperature_solution = solution;
 *       }
 *   
 *       timestep_number = 0;
 *       time_step = old_time_step = 0;
 *   
 *       double time = 0;
 *   
 *       do
 *         {
 *           pcout << "Timestep " << timestep_number
 *                 << ":  t=" << time / EquationData::year_in_seconds << " years"
 *                 << std::endl;
 *   
 *           assemble_stokes_system();
 *           build_stokes_preconditioner();
 *           assemble_temperature_matrix();
 *   
 *           solve();
 *   
 *           pcout << std::endl;
 *   
 *           if ((timestep_number == 0) &&
 *               (pre_refinement_step < parameters.initial_adaptive_refinement))
 *             {
 *               refine_mesh(parameters.initial_global_refinement +
 *                           parameters.initial_adaptive_refinement);
 *               ++pre_refinement_step;
 *               goto start_time_iteration;
 *             }
 *           else if ((timestep_number > 0) &&
 *                    (timestep_number % parameters.adaptive_refinement_interval ==
 *                     0))
 *             refine_mesh(parameters.initial_global_refinement +
 *                         parameters.initial_adaptive_refinement);
 *   
 *           if ((parameters.generate_graphical_output == true) &&
 *               (timestep_number % parameters.graphical_output_interval == 0))
 *             output_results();
 *   
 * @endcode
 * 
 * In order to speed up linear solvers, we extrapolate the solutions
 * from the old time levels to the new one. This gives a very good
 * initial guess, cutting the number of iterations needed in solvers
 * by more than one half. We do not need to extrapolate in the last
 * iteration, so if we reached the final time, we stop here.
 *         
 
 * 
 * As the last thing during a time step (before actually bumping up
 * the number of the time step), we check whether the current time
 * step number is divisible by 100, and if so we let the computing
 * timer print a summary of CPU times spent so far.
 * 
 * @code
 *           if (time > parameters.end_time * EquationData::year_in_seconds)
 *             break;
 *   
 *           TrilinosWrappers::MPI::BlockVector old_old_stokes_solution;
 *           old_old_stokes_solution      = old_stokes_solution;
 *           old_stokes_solution          = stokes_solution;
 *           old_old_temperature_solution = old_temperature_solution;
 *           old_temperature_solution     = temperature_solution;
 *           if (old_time_step > 0)
 *             {
 * @endcode
 * 
 * Trilinos sadd does not like ghost vectors even as input. Copy
 * into distributed vectors for now:
 * 
 * @code
 *               {
 *                 TrilinosWrappers::MPI::BlockVector distr_solution(stokes_rhs);
 *                 distr_solution = stokes_solution;
 *                 TrilinosWrappers::MPI::BlockVector distr_old_solution(stokes_rhs);
 *                 distr_old_solution = old_old_stokes_solution;
 *                 distr_solution.sadd(1. + time_step / old_time_step,
 *                                     -time_step / old_time_step,
 *                                     distr_old_solution);
 *                 stokes_solution = distr_solution;
 *               }
 *               {
 *                 TrilinosWrappers::MPI::Vector distr_solution(temperature_rhs);
 *                 distr_solution = temperature_solution;
 *                 TrilinosWrappers::MPI::Vector distr_old_solution(temperature_rhs);
 *                 distr_old_solution = old_old_temperature_solution;
 *                 distr_solution.sadd(1. + time_step / old_time_step,
 *                                     -time_step / old_time_step,
 *                                     distr_old_solution);
 *                 temperature_solution = distr_solution;
 *               }
 *             }
 *   
 *           if ((timestep_number > 0) && (timestep_number % 100 == 0))
 *             computing_timer.print_summary();
 *   
 *           time += time_step;
 *           ++timestep_number;
 *         }
 *       while (true);
 *   
 * @endcode
 * 
 * If we are generating graphical output, do so also for the last time
 * step unless we had just done so before we left the do-while loop
 * 
 * @code
 *       if ((parameters.generate_graphical_output == true) &&
 *           !((timestep_number - 1) % parameters.graphical_output_interval == 0))
 *         output_results();
 *     }
 *   } // namespace Step32
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_32-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * The main function is short as usual and very similar to the one in
 * @ref step_31 "step-31". Since we use a parameter file which is specified as an argument in
 * the command line, we have to read it in here and pass it on to the
 * Parameters class for parsing. If no filename is given in the command line,
 * we simply use the <code>step-32.prm</code> file which is distributed
 * together with the program.
 * 
 
 * 
 * Because 3d computations are simply very slow unless you throw a lot of
 * processors at them, the program defaults to 2d. You can get the 3d version
 * by changing the constant dimension below to 3.
 * 
 * @code
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         using namespace Step32;
 *         using namespace dealii;
 *   
 *         Utilities::MPI::MPI_InitFinalize mpi_initialization(
 *           argc, argv, numbers::invalid_unsigned_int);
 *   
 *         std::string parameter_filename;
 *         if (argc >= 2)
 *           parameter_filename = argv[1];
 *         else
 *           parameter_filename = "step-32.prm";
 *   
 *         const int                              dim = 2;
 *         BoussinesqFlowProblem<dim>::Parameters parameters(parameter_filename);
 *         BoussinesqFlowProblem<dim>             flow_problem(parameters);
 *         flow_problem.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="step_32-Results"></a><h1>Results</h1>
 
 
When run, the program simulates convection in 3d in much the same way
as @ref step_31 "step-31" did, though with an entirely different testcase.
 
 
<a name="step_32-Comparisonofresultswithstep31"></a><h3>Comparison of results with step-31</h3>
 
 
Before we go to this testcase, however, let us show a few results from a
slightly earlier version of this program that was solving exactly the
testcase we used in @ref step_31 "step-31", just that we now solve it in parallel and with
much higher resolution. We show these results mainly for comparison.
 
Here are two images that show this higher resolution if we choose a 3d
computation in <code>main()</code> and if we set
<code>initial_refinement=3</code> and
<code>n_pre_refinement_steps=4</code>. At the time steps shown, the
meshes had around 72,000 and 236,000 cells, for a total of 2,680,000
and 8,250,000 degrees of freedom, respectively, more than an order of
magnitude more than we had available in @ref step_31 "step-31":
 
<table align="center" class="doxtable">
  <tr>
    <td>
        <img src="https://www.dealii.org/images/steps/developer/step-32.3d.cube.0.png" alt="">
    </td>
  </tr>
  <tr>
    <td>
        <img src="https://www.dealii.org/images/steps/developer/step-32.3d.cube.1.png" alt="">
    </td>
  </tr>
</table>
 
The computation was done on a subset of 50 processors of the Brazos
cluster at Texas A&amp;M University.
 
 
<a name="step_32-Resultsfora2dcircularshelltestcase"></a><h3>Results for a 2d circular shell testcase</h3>
 
 
Next, we will run @ref step_32 "step-32" with the parameter file in the directory with one
change: we increase the final time to 1e9. Here we are using 16 processors. The
command to launch is (note that @ref step_32 "step-32".prm is the default):
 
<code>
<pre>
\$ mpirun -np 16 ./step-32
</pre>
</code>
 
Note that running a job on a cluster typically requires going through a job
scheduler, which we won't discuss here. The output will look roughly like
this:
 
<code>
<pre>
\$ mpirun -np 16 ./step-32
Number of active cells: 12,288 (on 6 levels)
Number of degrees of freedom: 186,624 (99,840+36,864+49,920)
 
Timestep 0:  t=0 years
 
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 41 iterations.
   Maximal velocity: 60.4935 cm/year
   Time step: 18166.9 years
   17 CG iterations for temperature
   Temperature range: 973 4273.16
 
Number of active cells: 15,921 (on 7 levels)
Number of degrees of freedom: 252,723 (136,640+47,763+68,320)
 
Timestep 0:  t=0 years
 
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 50 iterations.
   Maximal velocity: 60.3223 cm/year
   Time step: 10557.6 years
   19 CG iterations for temperature
   Temperature range: 973 4273.16
 
Number of active cells: 19,926 (on 8 levels)
Number of degrees of freedom: 321,246 (174,312+59,778+87,156)
 
Timestep 0:  t=0 years
 
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 50 iterations.
   Maximal velocity: 57.8396 cm/year
   Time step: 5453.78 years
   18 CG iterations for temperature
   Temperature range: 973 4273.16
 
Timestep 1:  t=5453.78 years
 
   Solving Stokes system... 49 iterations.
   Maximal velocity: 59.0231 cm/year
   Time step: 5345.86 years
   18 CG iterations for temperature
   Temperature range: 973 4273.16
 
Timestep 2:  t=10799.6 years
 
   Solving Stokes system... 24 iterations.
   Maximal velocity: 60.2139 cm/year
   Time step: 5241.51 years
   17 CG iterations for temperature
   Temperature range: 973 4273.16
 
[...]
 
Timestep 100:  t=272151 years
 
   Solving Stokes system... 21 iterations.
   Maximal velocity: 161.546 cm/year
   Time step: 1672.96 years
   17 CG iterations for temperature
   Temperature range: 973 4282.57
 
Number of active cells: 56,085 (on 8 levels)
Number of degrees of freedom: 903,408 (490,102+168,255+245,051)
 
 
 
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |       115s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |       103 |      2.82s |       2.5% |
| Assemble temperature matrices   |        12 |     0.452s |      0.39% |
| Assemble temperature rhs        |       103 |      11.5s |        10% |
| Build Stokes preconditioner     |        12 |      2.09s |       1.8% |
| Solve Stokes system             |       103 |      90.4s |        79% |
| Solve temperature system        |       103 |      1.53s |       1.3% |
| Postprocessing                  |         3 |     0.532s |      0.46% |
| Refine mesh structure, part 1   |        12 |      0.93s |      0.81% |
| Refine mesh structure, part 2   |        12 |     0.384s |      0.33% |
| Setup dof systems               |        13 |      2.96s |       2.6% |
+---------------------------------+-----------+------------+------------+
 
[...]
 
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start    |  9.14e+04s |            |
|                                             |            |            |
| Section                         | no. calls |  wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble Stokes system          |     47045 |  2.05e+03s |       2.2% |
| Assemble temperature matrices   |      4707 |       310s |      0.34% |
| Assemble temperature rhs        |     47045 |   8.7e+03s |       9.5% |
| Build Stokes preconditioner     |      4707 |  1.48e+03s |       1.6% |
| Solve Stokes system             |     47045 |  7.34e+04s |        80% |
| Solve temperature system        |     47045 |  1.46e+03s |       1.6% |
| Postprocessing                  |      1883 |       222s |      0.24% |
| Refine mesh structure, part 1   |      4706 |       641s |       0.7% |
| Refine mesh structure, part 2   |      4706 |       259s |      0.28% |
| Setup dof systems               |      4707 |  1.86e+03s |         2% |
+---------------------------------+-----------+------------+------------+
</pre>
</code>
 
The simulation terminates when the time reaches the 1 billion years
selected in the input file.  You can extrapolate from this how long a
simulation would take for a different final time (the time step size
ultimately settles on somewhere around 20,000 years, so computing for
two billion years will take 100,000 time steps, give or take 20%).  As
can be seen here, we spend most of the compute time in assembling
linear systems and &mdash; above all &mdash; in solving Stokes
systems.
 
 
To demonstrate the output we show the output from every 1250th time step here:
<table>
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-000.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-050.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-100.png" alt="">
    </td>
  </tr>
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-150.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-200.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-250.png" alt="">
    </td>
  </tr>
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-300.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-350.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-400.png" alt="">
    </td>
  </tr>
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-450.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-500.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-550.png" alt="">
    </td>
  </tr>
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-time-600.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-cells.png" alt="">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-32-2d-partition.png" alt="">
    </td>
  </tr>
</table>
 
The last two images show the grid as well as the partitioning of the mesh for
the same computation with 16 subdomains and 16 processors. The full dynamics of
this simulation are really only visible by looking at an animation, for example
the one <a
href="https://www.dealii.org/images/steps/developer/step-32-2d-temperature.webm">shown
on this site</a>. This image is well worth watching due to its artistic quality
and entrancing depiction of the evolution of the magma plumes.
 
If you watch the movie, you'll see that the convection pattern goes
through several stages: First, it gets rid of the instable temperature
layering with the hot material overlain by the dense cold
material. After this great driver is removed and we have a sort of
stable situation, a few blobs start to separate from the hot boundary
layer at the inner ring and rise up, with a few cold fingers also
dropping down from the outer boundary layer. During this phase, the solution
remains mostly symmetric, reflecting the 12-fold symmetry of the
original mesh. In a final phase, the fluid enters vigorous chaotic
stirring in which all symmetries are lost. This is a pattern that then
continues to dominate flow.
 
These different phases can also be identified if we look at the
maximal velocity as a function of time in the simulation:
 
<img src="https://www.dealii.org/images/steps/developer/step-32.2d.t_vs_vmax.png" alt="">
 
Here, the velocity (shown in centimeters per year) becomes very large,
to the order of several meters per year) at the beginning when the
temperature layering is instable. It then calms down to relatively
small values before picking up again in the chaotic stirring
regime. There, it remains in the range of 10-40 centimeters per year,
quite within the physically expected region.
 
 
<a name="step_32-Resultsfora3dsphericalshelltestcase"></a><h3>Results for a 3d spherical shell testcase</h3>
 
 
3d computations are very expensive computationally. Furthermore, as
seen above, interesting behavior only starts after quite a long time
requiring more CPU hours than is available on a typical
cluster. Consequently, rather than showing a complete simulation here,
let us simply show a couple of pictures we have obtained using the
successor to this program, called <i>ASPECT</i> (short for <i>Advanced
%Solver for Problems in Earth's ConvecTion</i>), that is being
developed independently of deal.II and that already incorporates some
of the extensions discussed below. The following two pictures show
isocontours of the temperature and the partition of the domain (along
with the mesh) onto 512 processors:
 
<p align="center">
<img src="https://www.dealii.org/images/steps/developer/step-32.3d-sphere.solution.png" alt="">
 
<img src="https://www.dealii.org/images/steps/developer/step-32.3d-sphere.partition.png" alt="">
</p>
 
 
<a name="step-32-extensions"></a>
<a name="step_32-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
 
 
There are many directions in which this program could be extended. As
mentioned at the end of the introduction, most of these are under active
development in the <i>ASPECT</i> (short for <i>Advanced %Solver for Problems
in Earth's ConvecTion</i>) code at the time this tutorial program is being
finished. Specifically, the following are certainly topics that one should
address to make the program more useful:
 
<ul>
  <li> <b>Adiabatic heating/cooling:</b>
  The temperature field we get in our simulations after a while
  is mostly constant with boundary layers at the inner and outer
  boundary, and streamers of cold and hot material mixing
  everything. Yet, this doesn't match our expectation that things
  closer to the earth core should be hotter than closer to the
  surface. The reason is that the energy equation we have used does
  not include a term that describes adiabatic cooling and heating:
  rock, like gas, heats up as you compress it. Consequently, material
  that rises up cools adiabatically, and cold material that sinks down
  heats adiabatically. The correct temperature equation would
  therefore look somewhat like this:
  @f{eqnarray*}{
    \frac{D T}{Dt}
    -
    \nabla \cdot \kappa \nabla T &=& \gamma + \tau\frac{Dp}{Dt},
  @f}
  or, expanding the advected derivative @f$\frac{D}{Dt} =
  \frac{\partial}{\partial t} + \mathbf u \cdot \nabla@f$:
  @f{eqnarray*}{
    \frac{\partial T}{\partial t}
    +
    {\mathbf u} \cdot \nabla T
    -
    \nabla \cdot \kappa \nabla T &=& \gamma +
    \tau\left\{\frac{\partial
    p}{\partial t} + \mathbf u \cdot \nabla p \right\}.
  @f}
  In other words, as pressure increases in a rock volume
  (@f$\frac{Dp}{Dt}>0@f$) we get an additional heat source, and vice
  versa.
 
  The time derivative of the pressure is a bit awkward to
  implement. If necessary, one could approximate using the fact
  outlined in the introduction that the pressure can be decomposed
  into a dynamic component due to temperature differences and the
  resulting flow, and a static component that results solely from the
  static pressure of the overlying rock. Since the latter is much
  bigger, one may approximate @f$p\approx p_{\text{static}}=-\rho_{\text{ref}}
  [1+\beta T_{\text{ref}}] \varphi@f$, and consequently
  @f$\frac{Dp}{Dt} \approx \left\{- \mathbf u \cdot \nabla \rho_{\text{ref}}
  [1+\beta T_{\text{ref}}]\varphi\right\} = \rho_{\text{ref}}
  [1+\beta T_{\text{ref}}] \mathbf u \cdot \mathbf g@f$.
  In other words, if the fluid is moving in the direction of gravity
  (downward) it will be compressed and because in that case @f$\mathbf u
  \cdot \mathbf g > 0@f$ we get a positive heat source. Conversely, the
  fluid will cool down if it moves against the direction of gravity.
 
<li> <b>Compressibility:</b>
  As already hinted at in the temperature model above,
  mantle rocks are not incompressible. Rather, given the enormous pressures in
  the earth mantle (at the core-mantle boundary, the pressure is approximately
  140 GPa, equivalent to 1,400,000 times atmospheric pressure), rock actually
  does compress to something around 1.5 times the density it would have
  at surface pressure. Modeling this presents any number of
  difficulties. Primarily, the mass conservation equation is no longer
  @f$\textrm{div}\;\mathbf u=0@f$ but should read
  @f$\textrm{div}(\rho\mathbf u)=0@f$ where the density @f$\rho@f$ is now no longer
  spatially constant but depends on temperature and pressure. A consequence is
  that the model is now no longer linear; a linearized version of the Stokes
  equation is also no longer symmetric requiring us to rethink preconditioners
  and, possibly, even the discretization. We won't go into detail here as to
  how this can be resolved.
 
<li> <b>Nonlinear material models:</b> As already hinted at in various places,
  material parameters such as the density, the viscosity, and the various
  thermal parameters are not constant throughout the earth mantle. Rather,
  they nonlinearly depend on the pressure and temperature, and in the case of
  the viscosity on the strain rate @f$\varepsilon(\mathbf u)@f$. For complicated
  models, the only way to solve such models accurately may be to actually
  iterate this dependence out in each time step, rather than simply freezing
  coefficients at values extrapolated from the previous time step(s).
 
<li> <b>Checkpoint/restart:</b> Running this program in 2d on a number of
  processors allows solving realistic models in a day or two. However, in 3d,
  compute times are so large that one runs into two typical problems: (i) On
  most compute clusters, the queuing system limits run times for individual
  jobs are to 2 or 3 days; (ii) losing the results of a computation due to
  hardware failures, misconfigurations, or power outages is a shame when
  running on hundreds of processors for a couple of days. Both of these
  problems can be addressed by periodically saving the state of the program
  and, if necessary, restarting the program at this point. This technique is
  commonly called <i>checkpoint/restart</i> and it requires that the entire
  state of the program is written to a permanent storage location (e.g. a hard
  drive). Given the complexity of the data structures of this program, this is
  not entirely trivial (it may also involve writing gigabytes or more of
  data), but it can be made easier by realizing that one can save the state
  between two time steps where it essentially only consists of the mesh and
  solution vectors; during restart one would then first re-enumerate degrees
  of freedom in the same way as done before and then re-assemble
  matrices. Nevertheless, given the distributed nature of the data structures
  involved here, saving and restoring the state of a program is not
  trivial. An additional complexity is introduced by the fact that one may
  want to change the number of processors between runs, for example because
  one may wish to continue computing on a mesh that is finer than the one used
  to precompute a starting temperature field at an intermediate time.
 
<li> <b>Predictive postprocessing:</b> The point of computations like this is
  not simply to solve the equations. Rather, it is typically the exploration
  of different physical models and their comparison with things that we can
  measure at the earth surface, in order to find which models are realistic
  and which are contradicted by reality. To this end, we need to compute
  quantities from our solution vectors that are related to what we can
  observe. Among these are, for example, heatfluxes at the surface of the
  earth, as well as seismic velocities throughout the mantle as these affect
  earthquake waves that are recorded by seismographs.
 
<li> <b>Better refinement criteria:</b> As can be seen above for the
3d case, the mesh in 3d is primarily refined along the inner
boundary. This is because the boundary layer there is stronger than
any other transition in the domain, leading us to refine there almost
exclusively and basically not at all following the plumes. One
certainly needs better refinement criteria to track the parts of the
solution we are really interested in better than the criterion used
here, namely the KellyErrorEstimator applied to the temperature, is
able to.
</ul>
 
 
There are many other ways to extend the current program. However, rather than
discussing them here, let us point to the much larger open
source code ASPECT (see https://aspect.geodynamics.org/ ) that constitutes the
further development of @ref step_32 "step-32" and that already includes many such possible
extensions.
 *
 *
<a name="step_32-PlainProg"></a>
<h1> The plain program</h1>
@include "step-32.cc"
*/