<commit_before>/* $Id$ */
/* Author: <NAME>, Uppsala University,
           <NAME>, Texas A&M University 2007, 2008, 2009 */
/*                                                                */
/*    Copyright (C) 2008, 2009 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

				 // @sect3{Include files}

				 // We include the functionality
				 // of these well-known deal.II
				 // library files and some C++
				 // header files.
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/function.h>
#include <base/utilities.h>
#include <base/conditional_ostream.h>
#include <base/work_stream.h>
#include <base/timer.h>

#include <lac/full_matrix.h>
#include <lac/solver_bicgstab.h>
#include <lac/solver_cg.h>
#include <lac/constraint_matrix.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_precondition.h>

#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/filtered_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>

#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>

#include <fe/fe_q.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <fe/mapping_q1.h>

#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numerics/solution_transfer.h>

#include <fstream>
#include <iostream>
#include <sstream>

				 // This is the only include file that is
				 // new: We use Trilinos for defining the
				 // %parallel partitioning of the matrices
				 // and vectors, and as explained in the
				 // introduction, an <code>Epetra_Map</code>
				 // is the Trilinos data structure for the
				 // definition of which part of a
				 // distributed vector is stored locally:
#include <Epetra_Map.h>


				 // Next, we import all deal.II
				 // names into global namespace:
using namespace dealii;

				 // @sect3{Equation data}

				 // This program is mainly an extension of
				 // step-31 to operate in %parallel, so the
				 // equation data remains the same.
namespace EquationData
{
  const double eta = 1;
  const double kappa = 1e-6;
  const double Rayleigh_number = 10;



  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;

      virtual void vector_value (const Point<dim> &p,
				 Vector<double>   &value) const;
  };



  template <int dim>
  double
  TemperatureInitialValues<dim>::value (const Point<dim>  &,
					const unsigned int) const
  {
				     /* Data for shell problem */
				     /*return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);*/

				     /* Data for cube problem */
    return 0.;
  }


  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);
  }



  template <int dim>
  class TemperatureRightHandSide : public Function<dim>
  {
    public:
      TemperatureRightHandSide () : Function<dim>(1) {}

      virtual double value (const Point<dim>   &p,
			    const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
				 Vector<double>   &value) const;
  };



  template <int dim>
  double
  TemperatureRightHandSide<dim>::value (const Point<dim>  &p,
					const unsigned int component) const
  {
				     /* Data for shell problem. */
				     /*    return 0; */

				     /* Data for cube problem. */
    Assert (component == 0,
	    ExcMessage ("Invalid operation for a scalar function."));
    
    Assert ((dim==2) || (dim==3), ExcNotImplemented());
    
    static const Point<dim> source_centers[3]
      = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
	  (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
	  (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
    static const double source_radius
      = (dim == 2 ? 1./32 : 1./8);
      
    return ((source_centers[0].distance (p) < source_radius)
	    ||
	    (source_centers[1].distance (p) < source_radius)
	    ||
	    (source_centers[2].distance (p) < source_radius)
	    ?
	    1
	    :
	    0);
  }


  template <int dim>
  void
  TemperatureRightHandSide<dim>::vector_value (const Point<dim> &p,
					       Vector<double>   &values) const
  {
    for (unsigned int c=0; c<this->n_components; ++c)
      values(c) = TemperatureRightHandSide<dim>::value (p, c);
  }
}



				 // @sect3{Linear solvers and preconditioners}

				 // In comparison to step-31, we did one
				 // change in the linear algebra of the
				 // problem: We exchange the
				 // <code>InverseMatrix</code> that
				 // previously held the approximation of the
				 // Schur complement by a preconditioner
				 // only (we will choose ILU in the
				 // application code below). This is the
				 // same trick we already did for the
				 // velocity block - the idea of this is
				 // that the outer iterations will
				 // eventually also make the inner
				 // approximation for the Schur complement
				 // good. If the preconditioner we're using
				 // is good enough, there will be no
				 // increase in the (outer) iteration
				 // count. All we need to do for
				 // implementing this change here is to give
				 // the respective variable in the
				 // BlockSchurPreconditioner class another
				 // name.
namespace LinearSolvers
{
  template <class PreconditionerA, class PreconditionerMp>
  class BlockSchurPreconditioner : public Subscriptor
  {
    public:
      BlockSchurPreconditioner (
	const TrilinosWrappers::BlockSparseMatrix  &S,
	const PreconditionerMp                     &Mppreconditioner,
	const PreconditionerA                      &Apreconditioner);

      void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
		  const TrilinosWrappers::MPI::BlockVector &src) const;

    private:
      const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
      const PreconditionerMp &mp_preconditioner;
      const PreconditionerA  &a_preconditioner;
      mutable TrilinosWrappers::MPI::Vector tmp;
  };



  template <class PreconditionerA, class PreconditionerMp>
  BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
  BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
			   const PreconditionerMp                    &Mppreconditioner,
			   const PreconditionerA                     &Apreconditioner)
		  :
		  stokes_matrix     (&S),
		  mp_preconditioner (Mppreconditioner),
		  a_preconditioner  (Apreconditioner),
		  tmp               (stokes_matrix->block(1,1).row_partitioner())
  {}



  template <class PreconditionerA, class PreconditionerMp>
  void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
    TrilinosWrappers::MPI::BlockVector       &dst,
    const TrilinosWrappers::MPI::BlockVector &src) const
  {
    a_preconditioner.vmult (dst.block(0), src.block(0));
    stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
    tmp *= -1;
    mp_preconditioner.vmult (dst.block(1), tmp);
  }
}



				 // @sect3{Definition of assembly data structures}
				 //
				 // As described in the introduction, we will
				 // use the WorkStream mechanism discussed in
				 // the @ref threads module 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.
				 //
				 // Using the program in step-31, we have
				 // 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 organize the scratch
				 // arrays and a CopyData arrays for each
				 // of those four assembly components
				 // using a <code>struct</code>
				 // environment.
				 //
				 // Regarding the Scratch array, each
				 // struct is equipped with a constructor
				 // that create an FEValues object for a
				 // @ref FiniteElement "finite element", a
				 // @ref Quadrature "quadrature formula"
				 // and some
				 // @ref UpdateFlags "update flags".
				 // 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 improve
				 // performance of assembly.
namespace Assembly
{
  namespace Scratch
  {
    template <int dim>
    struct StokesPreconditioner 
    {
	StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
			      const Quadrature<dim>    &stokes_quadrature,
			      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 UpdateFlags         update_flags)
		    :
		    stokes_fe_values (stokes_fe, stokes_quadrature,
				      update_flags),
		    grad_phi_u (stokes_fe.dofs_per_cell),
		    phi_p (stokes_fe.dofs_per_cell)
    {}



    template <int dim>
    StokesPreconditioner<dim>::
    StokesPreconditioner (const StokesPreconditioner &scratch)
		    :
		    stokes_fe_values (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)
    {}



				     // Observe that we derive the
				     // StokesSystem scratch array from the
				     // StokesPreconditioner array. 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.
    template <int dim>
    struct StokesSystem : public StokesPreconditioner<dim>
    {
	StokesSystem (const FiniteElement<dim> &stokes_fe,
		      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 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,
					       stokes_update_flags),
		    temperature_fe_values (temperature_fe, stokes_quadrature,
					   temperature_update_flags),
		    phi_u (stokes_fe.dofs_per_cell),
		    grads_phi_u (stokes_fe.dofs_per_cell),
		    div_phi_u (stokes_fe.dofs_per_cell),
		    old_temperature_values (stokes_quadrature.n_quadrature_points)
    {}


    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const StokesSystem<dim> &scratch)
		    :
		    StokesPreconditioner<dim> (scratch),
		    temperature_fe_values (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)
    {}



    template <int dim>
    struct TemperatureMatrix 
    {
	TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
			   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 Quadrature<dim>    &temperature_quadrature)
		    :
		    temperature_fe_values (temperature_fe, temperature_quadrature,
					   update_values    | update_gradients |
					   update_JxW_values),
		    phi_T (temperature_fe.dofs_per_cell),
		    grad_phi_T (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const TemperatureMatrix &scratch)
		    :
		    temperature_fe_values (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)
    {}


    template <int dim>
    struct TemperatureRHS 
    {
	TemperatureRHS (const FiniteElement<dim> &temperature_fe,
			const FiniteElement<dim> &stokes_fe,
			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<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;

	std::vector<double>         gamma_values;
    };

    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const FiniteElement<dim> &temperature_fe,
		    const FiniteElement<dim> &stokes_fe,
		    const Quadrature<dim>    &quadrature)
		    :
		    temperature_fe_values (temperature_fe, quadrature,
					   update_values    |
					   update_gradients |
					   update_hessians  |
					   update_quadrature_points |
					   update_JxW_values),
		    stokes_fe_values (stokes_fe, quadrature,
				      update_values),
		    phi_T (temperature_fe.dofs_per_cell),
		    grad_phi_T (temperature_fe.dofs_per_cell),
		
		    old_velocity_values (quadrature.n_quadrature_points),
		    old_old_velocity_values (quadrature.n_quadrature_points),
		
		    old_temperature_values (quadrature.n_quadrature_points),
		    old_old_temperature_values(quadrature.n_quadrature_points),
		    old_temperature_grads(quadrature.n_quadrature_points),
		    old_old_temperature_grads(quadrature.n_quadrature_points),
		    old_temperature_laplacians(quadrature.n_quadrature_points),
		    old_old_temperature_laplacians(quadrature.n_quadrature_points),
		
		    gamma_values (quadrature.n_quadrature_points)		
    {}


    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const TemperatureRHS &scratch)
		    :
		    temperature_fe_values (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_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_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),
		
		    gamma_values (scratch.gamma_values)
    {}
  }


				   // The CopyData arrays are similar to the
				   // Scratch arrays. They provide a
				   // constructor, a copy operation, and
				   // some arrays for local matrix, local
				   // vectors and the relation between local
				   // and global degrees of freedom (aka
				   // <code>local_dof_indices</code>).
  namespace CopyData
  {
    template <int dim>
    struct StokesPreconditioner 
    {
	StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
	StokesPreconditioner (const StokesPreconditioner &data);

	FullMatrix<double>          local_matrix;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    StokesPreconditioner<dim>::
    StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
		    :
		    local_matrix (stokes_fe.dofs_per_cell, 
				  stokes_fe.dofs_per_cell),
		    local_dof_indices (stokes_fe.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);
	StokesSystem (const StokesSystem<dim> &data);
      
	Vector<double> local_rhs;
    };
  

    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const FiniteElement<dim> &stokes_fe)
		    :
		    StokesPreconditioner<dim> (stokes_fe),
		    local_rhs (stokes_fe.dofs_per_cell)
    {}


    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const StokesSystem<dim> &data)
		    :
		    StokesPreconditioner<dim> (data),
		    local_rhs (data.local_rhs)
    {}



    template <int dim>
    struct TemperatureMatrix 
    {
	TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
	TemperatureMatrix (const TemperatureMatrix &data);
	
	FullMatrix<double>          local_mass_matrix;
	FullMatrix<double>          local_stiffness_matrix;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
		    :
		    local_mass_matrix (temperature_fe.dofs_per_cell,
				       temperature_fe.dofs_per_cell),
		    local_stiffness_matrix (temperature_fe.dofs_per_cell,
					    temperature_fe.dofs_per_cell),
		    local_dof_indices (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const TemperatureMatrix &data)
		    :
		    local_mass_matrix (data.local_mass_matrix),
		    local_stiffness_matrix (data.local_stiffness_matrix),
		    local_dof_indices (data.local_dof_indices)
    {}


    template <int dim>
    struct TemperatureRHS 
    {
	TemperatureRHS (const FiniteElement<dim> &temperature_fe);
	TemperatureRHS (const TemperatureRHS &data);
	
	Vector<double>              local_rhs;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const FiniteElement<dim> &temperature_fe)
		    :
		    local_rhs (temperature_fe.dofs_per_cell),
		    local_dof_indices (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const TemperatureRHS &data)
		    :
		    local_rhs (data.local_rhs),
		    local_dof_indices (data.local_dof_indices)
    {}
  }  
}



				 // @sect3{The <code>BoussinesqFlowProblem</code> class template}
				 // 
				 // This is the declaration of the main
				 // class. It is very similar to
				 // step-31. 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.
				 // 
				 // Moreover, we include an MPI communicator
				 // and an Epetra_Map (see the introduction)
				 // that are needed for communication and
				 // data exchange if the Trilinos matrices
				 // and vectors are distributed over several
				 // processors. Finally, 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 produce the same output
				 // it will just be replicated many times
				 // over; with the ConditionalOStream class,
				 // only the output generated by one task
				 // will actually be printed to screen,
				 // whereas the output by all the other
				 // threads will simply be forgotten.
				 //
				 // 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). For
				 // the matrices, we happen to use the same
				 // class names for %parallel and sequential
				 // data structures, i.e. all matrices will
				 // actually be considered parallel
				 // below. On the other hand, for vectors,
				 // only those from namespace
				 // TrilinosWrappers::MPI are actually
				 // distributed. In particular, we will
				 // frequently have to query velocities and
				 // temperatures at arbitrary quadrature
				 // points; consequently, rather than
				 // "localizing" a vector whenever we need a
				 // localized vector, we solve linear
				 // systems in %parallel but then immediately
				 // localize the solution for further
				 // processing. The various
				 // <code>*_solution</code> vectors are
				 // therefore filled immediately after
				 // solving their respective linear system
				 // in %parallel.
				 //
				 // The only other new data member is
				 // <code>computing_timer</code>. Its class
				 // type, TimerOutput, can be 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. At the end of
				 // the program, the destructor of the
				 // TimerOutput class will automatically
				 // produce a nice summary of the times spent
				 // in all the sections. For this output, one
				 // can choose whether wall clock or CPU times
				 // are to be printed, as well as whether we
				 // want a percentage breakdown of where the
				 // time was spent (this choice is made in the
				 // constructor of TimerOutput, which is
				 // called in the constructor of
				 // <code>BoussinesqFlowProblem</code>). You
				 // can take a look at the output generated
				 // from this variable in the results section
				 // of this tutorial program.
template <int dim>
class BoussinesqFlowProblem
{
  public:
    BoussinesqFlowProblem ();
    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);
    void project_temperature_field ();
    double get_maximal_velocity () const;
    std::pair<double,double> get_extrapolated_temperature_range () const;
    void solve ();
    void output_results () const;
    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<double>          &gamma_values,
		      const double                        global_u_infty,
		      const double                        global_T_variation,
		      const double                        cell_diameter) const;


    const Epetra_Comm                  &trilinos_communicator;
    
    ConditionalOStream                  pcout;

    Triangulation<dim>                  triangulation;
    double                              global_Omega_diameter;

    const unsigned int                  stokes_degree;
    FESystem<dim>                       stokes_fe;
    DoFHandler<dim>                     stokes_dof_handler;
    ConstraintMatrix                    stokes_constraints;

    std::vector<Epetra_Map>             stokes_partitioner;
    TrilinosWrappers::BlockSparseMatrix stokes_matrix;
    TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;

    TrilinosWrappers::BlockVector       stokes_solution;
    TrilinosWrappers::BlockVector       old_stokes_solution;
    TrilinosWrappers::MPI::BlockVector  stokes_rhs;


    const unsigned int                  temperature_degree;    
    FE_Q<dim>                           temperature_fe;
    DoFHandler<dim>                     temperature_dof_handler;
    ConstraintMatrix                    temperature_constraints;

    Epetra_Map                          temperature_partitioner;
    TrilinosWrappers::SparseMatrix      temperature_mass_matrix;
    TrilinosWrappers::SparseMatrix      temperature_stiffness_matrix;
    TrilinosWrappers::SparseMatrix      temperature_matrix;

    TrilinosWrappers::Vector            temperature_solution;
    TrilinosWrappers::Vector            old_temperature_solution;
    TrilinosWrappers::Vector            old_old_temperature_solution;
    TrilinosWrappers::MPI::Vector       temperature_rhs;


    double time_step;
    double old_time_step;
    unsigned int timestep_number;

    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;

    bool rebuild_stokes_matrix;
    bool rebuild_stokes_preconditioner;
    bool rebuild_temperature_matrices;
    bool rebuild_temperature_preconditioner;

    TimerOutput computing_timer;

    void setup_stokes_matrix ();
    void setup_stokes_preconditioner ();
    void setup_temperature_matrices ();
    
    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 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);
};


				 // @sect3{BoussinesqFlowProblem class implementation}

				 // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem}
				 // 
				 // The constructor of the problem is very
				 // similar to the constructor in
				 // step-31. What is different is the
				 // %parallel communication: Trilins uses a
				 // message passing interface (MPI) for
				 // data distribution. When entering the
				 // BoussinesqFlowProblem class, we have
				 // to decide how the parallization is to
				 // be done. We choose a rather simple
				 // strategy and let all processors
				 // running the program work together,
				 // specified by the communicator
				 // <code>comm_world()</code>. Next, we
				 // create some modified output stream as
				 // we already did in step-18. In MPI, all
				 // the processors run the same program
				 // individually (they simply operate on
				 // different chunks of data and exchange
				 // some data from time to time). Since we
				 // do not want each processor to write
				 // the same information to screen (like
				 // the number of degrees of freedom), we
				 // only use one processor for writing
				 // that output to terminal windows. The
				 // implementation of this idea is to<commit_msg>One typo fix. Enough for today.<commit_after>/* $Id$ */
/* Author: <NAME>, Uppsala University,
           <NAME>, Texas A&M University 2007, 2008, 2009 */
/*                                                                */
/*    Copyright (C) 2008, 2009 by the deal.II authors */
/*                                                                */
/*    This file is subject to QPL and may not be  distributed     */
/*    without copyright and license information. Please refer     */
/*    to the file deal.II/doc/license.html for the  text  and     */
/*    further information on this license.                        */

				 // @sect3{Include files}

				 // We include the functionality
				 // of these well-known deal.II
				 // library files and some C++
				 // header files.
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/function.h>
#include <base/utilities.h>
#include <base/conditional_ostream.h>
#include <base/work_stream.h>
#include <base/timer.h>

#include <lac/full_matrix.h>
#include <lac/solver_bicgstab.h>
#include <lac/solver_cg.h>
#include <lac/constraint_matrix.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_precondition.h>

#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/filtered_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>

#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>

#include <fe/fe_q.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <fe/mapping_q1.h>

#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numerics/solution_transfer.h>

#include <fstream>
#include <iostream>
#include <sstream>

				 // This is the only include file that is
				 // new: We use Trilinos for defining the
				 // %parallel partitioning of the matrices
				 // and vectors, and as explained in the
				 // introduction, an <code>Epetra_Map</code>
				 // is the Trilinos data structure for the
				 // definition of which part of a
				 // distributed vector is stored locally:
#include <Epetra_Map.h>


				 // Next, we import all deal.II
				 // names into global namespace:
using namespace dealii;

				 // @sect3{Equation data}

				 // This program is mainly an extension of
				 // step-31 to operate in %parallel, so the
				 // equation data remains the same.
namespace EquationData
{
  const double eta = 1;
  const double kappa = 1e-6;
  const double Rayleigh_number = 10;



  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;

      virtual void vector_value (const Point<dim> &p,
				 Vector<double>   &value) const;
  };



  template <int dim>
  double
  TemperatureInitialValues<dim>::value (const Point<dim>  &,
					const unsigned int) const
  {
				     /* Data for shell problem */
				     /*return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);*/

				     /* Data for cube problem */
    return 0.;
  }


  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);
  }



  template <int dim>
  class TemperatureRightHandSide : public Function<dim>
  {
    public:
      TemperatureRightHandSide () : Function<dim>(1) {}

      virtual double value (const Point<dim>   &p,
			    const unsigned int  component = 0) const;

      virtual void vector_value (const Point<dim> &p,
				 Vector<double>   &value) const;
  };



  template <int dim>
  double
  TemperatureRightHandSide<dim>::value (const Point<dim>  &p,
					const unsigned int component) const
  {
				     /* Data for shell problem. */
				     /*    return 0; */

				     /* Data for cube problem. */
    Assert (component == 0,
	    ExcMessage ("Invalid operation for a scalar function."));
    
    Assert ((dim==2) || (dim==3), ExcNotImplemented());
    
    static const Point<dim> source_centers[3]
      = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
	  (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
	  (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
    static const double source_radius
      = (dim == 2 ? 1./32 : 1./8);
      
    return ((source_centers[0].distance (p) < source_radius)
	    ||
	    (source_centers[1].distance (p) < source_radius)
	    ||
	    (source_centers[2].distance (p) < source_radius)
	    ?
	    1
	    :
	    0);
  }


  template <int dim>
  void
  TemperatureRightHandSide<dim>::vector_value (const Point<dim> &p,
					       Vector<double>   &values) const
  {
    for (unsigned int c=0; c<this->n_components; ++c)
      values(c) = TemperatureRightHandSide<dim>::value (p, c);
  }
}



				 // @sect3{Linear solvers and preconditioners}

				 // In comparison to step-31, we did one
				 // change in the linear algebra of the
				 // problem: We exchange the
				 // <code>InverseMatrix</code> that
				 // previously held the approximation of the
				 // Schur complement by a preconditioner
				 // only (we will choose ILU in the
				 // application code below). This is the
				 // same trick we already did for the
				 // velocity block - the idea of this is
				 // that the outer iterations will
				 // eventually also make the inner
				 // approximation for the Schur complement
				 // good. If the preconditioner we're using
				 // is good enough, there will be no
				 // increase in the (outer) iteration
				 // count. All we need to do for
				 // implementing this change here is to give
				 // the respective variable in the
				 // BlockSchurPreconditioner class another
				 // name.
namespace LinearSolvers
{
  template <class PreconditionerA, class PreconditionerMp>
  class BlockSchurPreconditioner : public Subscriptor
  {
    public:
      BlockSchurPreconditioner (
	const TrilinosWrappers::BlockSparseMatrix  &S,
	const PreconditionerMp                     &Mppreconditioner,
	const PreconditionerA                      &Apreconditioner);

      void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
		  const TrilinosWrappers::MPI::BlockVector &src) const;

    private:
      const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
      const PreconditionerMp &mp_preconditioner;
      const PreconditionerA  &a_preconditioner;
      mutable TrilinosWrappers::MPI::Vector tmp;
  };



  template <class PreconditionerA, class PreconditionerMp>
  BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
  BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S,
			   const PreconditionerMp                    &Mppreconditioner,
			   const PreconditionerA                     &Apreconditioner)
		  :
		  stokes_matrix     (&S),
		  mp_preconditioner (Mppreconditioner),
		  a_preconditioner  (Apreconditioner),
		  tmp               (stokes_matrix->block(1,1).row_partitioner())
  {}



  template <class PreconditionerA, class PreconditionerMp>
  void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
    TrilinosWrappers::MPI::BlockVector       &dst,
    const TrilinosWrappers::MPI::BlockVector &src) const
  {
    a_preconditioner.vmult (dst.block(0), src.block(0));
    stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
    tmp *= -1;
    mp_preconditioner.vmult (dst.block(1), tmp);
  }
}



				 // @sect3{Definition of assembly data structures}
				 //
				 // As described in the introduction, we will
				 // use the WorkStream mechanism discussed in
				 // the @ref threads module 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.
				 //
				 // Using the program in step-31, we have
				 // 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 organize the scratch
				 // arrays and a CopyData arrays for each
				 // of those four assembly components
				 // using a <code>struct</code>
				 // environment.
				 //
				 // Regarding the Scratch array, each
				 // struct is equipped with a constructor
				 // that create an FEValues object for a
				 // @ref FiniteElement "finite element", a
				 // @ref Quadrature "quadrature formula"
				 // and some
				 // @ref UpdateFlags "update flags".
				 // 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 improve
				 // performance of assembly.
namespace Assembly
{
  namespace Scratch
  {
    template <int dim>
    struct StokesPreconditioner 
    {
	StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
			      const Quadrature<dim>    &stokes_quadrature,
			      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 UpdateFlags         update_flags)
		    :
		    stokes_fe_values (stokes_fe, stokes_quadrature,
				      update_flags),
		    grad_phi_u (stokes_fe.dofs_per_cell),
		    phi_p (stokes_fe.dofs_per_cell)
    {}



    template <int dim>
    StokesPreconditioner<dim>::
    StokesPreconditioner (const StokesPreconditioner &scratch)
		    :
		    stokes_fe_values (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)
    {}



				     // Observe that we derive the
				     // StokesSystem scratch array from the
				     // StokesPreconditioner array. 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.
    template <int dim>
    struct StokesSystem : public StokesPreconditioner<dim>
    {
	StokesSystem (const FiniteElement<dim> &stokes_fe,
		      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 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,
					       stokes_update_flags),
		    temperature_fe_values (temperature_fe, stokes_quadrature,
					   temperature_update_flags),
		    phi_u (stokes_fe.dofs_per_cell),
		    grads_phi_u (stokes_fe.dofs_per_cell),
		    div_phi_u (stokes_fe.dofs_per_cell),
		    old_temperature_values (stokes_quadrature.n_quadrature_points)
    {}


    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const StokesSystem<dim> &scratch)
		    :
		    StokesPreconditioner<dim> (scratch),
		    temperature_fe_values (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)
    {}



    template <int dim>
    struct TemperatureMatrix 
    {
	TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
			   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 Quadrature<dim>    &temperature_quadrature)
		    :
		    temperature_fe_values (temperature_fe, temperature_quadrature,
					   update_values    | update_gradients |
					   update_JxW_values),
		    phi_T (temperature_fe.dofs_per_cell),
		    grad_phi_T (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const TemperatureMatrix &scratch)
		    :
		    temperature_fe_values (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)
    {}


    template <int dim>
    struct TemperatureRHS 
    {
	TemperatureRHS (const FiniteElement<dim> &temperature_fe,
			const FiniteElement<dim> &stokes_fe,
			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<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;

	std::vector<double>         gamma_values;
    };

    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const FiniteElement<dim> &temperature_fe,
		    const FiniteElement<dim> &stokes_fe,
		    const Quadrature<dim>    &quadrature)
		    :
		    temperature_fe_values (temperature_fe, quadrature,
					   update_values    |
					   update_gradients |
					   update_hessians  |
					   update_quadrature_points |
					   update_JxW_values),
		    stokes_fe_values (stokes_fe, quadrature,
				      update_values),
		    phi_T (temperature_fe.dofs_per_cell),
		    grad_phi_T (temperature_fe.dofs_per_cell),
		
		    old_velocity_values (quadrature.n_quadrature_points),
		    old_old_velocity_values (quadrature.n_quadrature_points),
		
		    old_temperature_values (quadrature.n_quadrature_points),
		    old_old_temperature_values(quadrature.n_quadrature_points),
		    old_temperature_grads(quadrature.n_quadrature_points),
		    old_old_temperature_grads(quadrature.n_quadrature_points),
		    old_temperature_laplacians(quadrature.n_quadrature_points),
		    old_old_temperature_laplacians(quadrature.n_quadrature_points),
		
		    gamma_values (quadrature.n_quadrature_points)		
    {}


    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const TemperatureRHS &scratch)
		    :
		    temperature_fe_values (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_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_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),
		
		    gamma_values (scratch.gamma_values)
    {}
  }


				   // The CopyData arrays are similar to the
				   // Scratch arrays. They provide a
				   // constructor, a copy operation, and
				   // some arrays for local matrix, local
				   // vectors and the relation between local
				   // and global degrees of freedom (aka
				   // <code>local_dof_indices</code>).
  namespace CopyData
  {
    template <int dim>
    struct StokesPreconditioner 
    {
	StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
	StokesPreconditioner (const StokesPreconditioner &data);

	FullMatrix<double>          local_matrix;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    StokesPreconditioner<dim>::
    StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
		    :
		    local_matrix (stokes_fe.dofs_per_cell, 
				  stokes_fe.dofs_per_cell),
		    local_dof_indices (stokes_fe.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);
	StokesSystem (const StokesSystem<dim> &data);
      
	Vector<double> local_rhs;
    };
  

    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const FiniteElement<dim> &stokes_fe)
		    :
		    StokesPreconditioner<dim> (stokes_fe),
		    local_rhs (stokes_fe.dofs_per_cell)
    {}


    template <int dim>
    StokesSystem<dim>::
    StokesSystem (const StokesSystem<dim> &data)
		    :
		    StokesPreconditioner<dim> (data),
		    local_rhs (data.local_rhs)
    {}



    template <int dim>
    struct TemperatureMatrix 
    {
	TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
	TemperatureMatrix (const TemperatureMatrix &data);
	
	FullMatrix<double>          local_mass_matrix;
	FullMatrix<double>          local_stiffness_matrix;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
		    :
		    local_mass_matrix (temperature_fe.dofs_per_cell,
				       temperature_fe.dofs_per_cell),
		    local_stiffness_matrix (temperature_fe.dofs_per_cell,
					    temperature_fe.dofs_per_cell),
		    local_dof_indices (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureMatrix<dim>::
    TemperatureMatrix (const TemperatureMatrix &data)
		    :
		    local_mass_matrix (data.local_mass_matrix),
		    local_stiffness_matrix (data.local_stiffness_matrix),
		    local_dof_indices (data.local_dof_indices)
    {}


    template <int dim>
    struct TemperatureRHS 
    {
	TemperatureRHS (const FiniteElement<dim> &temperature_fe);
	TemperatureRHS (const TemperatureRHS &data);
	
	Vector<double>              local_rhs;
	std::vector<unsigned int>   local_dof_indices;
    };

    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const FiniteElement<dim> &temperature_fe)
		    :
		    local_rhs (temperature_fe.dofs_per_cell),
		    local_dof_indices (temperature_fe.dofs_per_cell)
    {}


    template <int dim>
    TemperatureRHS<dim>::
    TemperatureRHS (const TemperatureRHS &data)
		    :
		    local_rhs (data.local_rhs),
		    local_dof_indices (data.local_dof_indices)
    {}
  }  
}



				 // @sect3{The <code>BoussinesqFlowProblem</code> class template}
				 // 
				 // This is the declaration of the main
				 // class. It is very similar to
				 // step-31. 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.
				 // 
				 // Moreover, we include an MPI communicator
				 // and an Epetra_Map (see the introduction)
				 // that are needed for communication and
				 // data exchange if the Trilinos matrices
				 // and vectors are distributed over several
				 // processors. Finally, 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 produce the same output
				 // it will just be replicated many times
				 // over; with the ConditionalOStream class,
				 // only the output generated by one task
				 // will actually be printed to screen,
				 // whereas the output by all the other
				 // threads will simply be forgotten.
				 //
				 // 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). For
				 // the matrices, we happen to use the same
				 // class names for %parallel and sequential
				 // data structures, i.e. all matrices will
				 // actually be considered parallel
				 // below. On the other hand, for vectors,
				 // only those from namespace
				 // TrilinosWrappers::MPI are actually
				 // distributed. In particular, we will
				 // frequently have to query velocities and
				 // temperatures at arbitrary quadrature
				 // points; consequently, rather than
				 // "localizing" a vector whenever we need a
				 // localized vector, we solve linear
				 // systems in %parallel but then immediately
				 // localize the solution for further
				 // processing. The various
				 // <code>*_solution</code> vectors are
				 // therefore filled immediately after
				 // solving their respective linear system
				 // in %parallel.
				 //
				 // The only other new data member is
				 // <code>computing_timer</code>. Its class
				 // type, TimerOutput, can be 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. At the end of
				 // the program, the destructor of the
				 // TimerOutput class will automatically
				 // produce a nice summary of the times spent
				 // in all the sections. For this output, one
				 // can choose whether wall clock or CPU times
				 // are to be printed, as well as whether we
				 // want a percentage breakdown of where the
				 // time was spent (this choice is made in the
				 // constructor of TimerOutput, which is
				 // called in the constructor of
				 // <code>BoussinesqFlowProblem</code>). You
				 // can take a look at the output generated
				 // from this variable in the results section
				 // of this tutorial program.
template <int dim>
class BoussinesqFlowProblem
{
  public:
    BoussinesqFlowProblem ();
    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);
    void project_temperature_field ();
    double get_maximal_velocity () const;
    std::pair<double,double> get_extrapolated_temperature_range () const;
    void solve ();
    void output_results () const;
    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<double>          &gamma_values,
		      const double                        global_u_infty,
		      const double                        global_T_variation,
		      const double                        cell_diameter) const;


    const Epetra_Comm                  &trilinos_communicator;
    
    ConditionalOStream                  pcout;

    Triangulation<dim>                  triangulation;
    double                              global_Omega_diameter;

    const unsigned int                  stokes_degree;
    FESystem<dim>                       stokes_fe;
    DoFHandler<dim>                     stokes_dof_handler;
    ConstraintMatrix                    stokes_constraints;

    std::vector<Epetra_Map>             stokes_partitioner;
    TrilinosWrappers::BlockSparseMatrix stokes_matrix;
    TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;

    TrilinosWrappers::BlockVector       stokes_solution;
    TrilinosWrappers::BlockVector       old_stokes_solution;
    TrilinosWrappers::MPI::BlockVector  stokes_rhs;


    const unsigned int                  temperature_degree;    
    FE_Q<dim>                           temperature_fe;
    DoFHandler<dim>                     temperature_dof_handler;
    ConstraintMatrix                    temperature_constraints;

    Epetra_Map                          temperature_partitioner;
    TrilinosWrappers::SparseMatrix      temperature_mass_matrix;
    TrilinosWrappers::SparseMatrix      temperature_stiffness_matrix;
    TrilinosWrappers::SparseMatrix      temperature_matrix;

    TrilinosWrappers::Vector            temperature_solution;
    TrilinosWrappers::Vector            old_temperature_solution;
    TrilinosWrappers::Vector            old_old_temperature_solution;
    TrilinosWrappers::MPI::Vector       temperature_rhs;


    double time_step;
    double old_time_step;
    unsigned int timestep_number;

    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;

    bool rebuild_stokes_matrix;
    bool rebuild_stokes_preconditioner;
    bool rebuild_temperature_matrices;
    bool rebuild_temperature_preconditioner;

    TimerOutput computing_timer;

    void setup_stokes_matrix ();
    void setup_stokes_preconditioner ();
    void setup_temperature_matrices ();
    
    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 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);
};


				 // @sect3{BoussinesqFlowProblem class implementation}

				 // @sect4{BoussinesqFlowProblem::BoussinesqFlowProblem}
				 // 
				 // The constructor of the problem is very
				 // similar to the constructor in
				 // 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 parallization is to
				 // be done. We choose a rather simple
				 // strategy and let all processors
				 // running the program work together,
				 // specified by the communicator
				 // <code>comm_world()</code>. Next, we
				 // create some modified output stream as
				 // we already did in step-18. In MPI, all
				 // the processors run the same program
				 // individually (they simply operate on
				 // different chunks of data and exchange
				 // some data from time to time). Since we
				 // do not want each processor to write
				 // the same information to screen (like
				 // the number of degrees of freedom), we
				 // only use one processor for writing
				 // that output to terminal windows. The
				 // implementation of this idea is to
{'fraction_non_alphanumeric': 0.08777468085434023, 'fraction_numerical': 0.003911362971403933, 'mean_word_length': 2.782178217821782, 'pattern_counts': {'":': 0, '<': 551, '<?xml version=': 0, '>': 547, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7539156', 'n_tokens_mistral': 22337, 'n_tokens_neox': 18510, 'n_words': 5052}