parallel_in_time.h

Go to the documentation of this file.

 
 *     double tstop;              /* evolve to this time*/
 *     int level;
 *     double deltaT;
 *   
 *     int index;
 *     braid_StepStatusGetLevel(status, &level);
 *     braid_StepStatusGetTstartTstop(status, &tstart, &tstop);
 *     braid_StepStatusGetTIndex(status, &index);
 *   
 *     deltaT = tstop - tstart;
 *   
 *     ::Vector<double>& solution = u->data;
 *   
 *     HeatEquation<2>& heateq = app->eq;
 *   
 *     heateq.step(solution, deltaT, tstart, index);
 *   
 *     return 0;
 *   }
 *   
 *   
 * @endcode
 * 
 * In this function we initialize a vector at an arbitrary time.
 * At this point we don't know anything about what the solution
 * looks like, and we can really initialize to anything, so in
 * this case use reinit to initialize the memory and set the
 * values to zero.
 * 
 * @code
 *   int
 *   my_Init(braid_App     app,
 *           double        t,
 *           braid_Vector *u_ptr)
 *   {
 *     my_Vector *u = new(my_Vector);
 *     int size = app->eq.size();
 *     u->data.reinit(size);
 *   
 *     app->eq.initialize(t, u->data);
 *   
 *     *u_ptr = u;
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * Here we need to copy the vector u into the vector v. We do this
 * by allocating a new vector, then reinitializing the deal.ii
 * vector to the correct size. The deal.ii reinitialization sets
 * every value to zero, so next we need to iterate over the vector
 * u and copy the values to the new vector v.
 * 
 * @code
 *   int
 *   my_Clone(braid_App     app,
 *            braid_Vector  u,
 *            braid_Vector *v_ptr)
 *   {
 *     UNUSED(app);
 *     my_Vector *v = new(my_Vector);
 *     int size = u->data.size();
 *     v->data.reinit(size);
 *     for(size_t i=0, end=v->data.size(); i != end; ++i)
 *       {
 *         v->data[i] = u->data[i];
 *       }
 *     *v_ptr = v;
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * Here we need to free the memory used by vector u. This is
 * pretty simple since the deal.ii vector is stored inside the
 * XBraid vector, so we just delete the XBraid vector u and it
 * puts the deal.ii vector out of scope and releases its memory.
 * 
 * @code
 *   int
 *   my_Free(braid_App    app,
 *           braid_Vector u)
 *   {
 *     UNUSED(app);
 *     delete u;
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This is to perform an axpy type operation. That is to say we
 * do @f$y = \alpha x + \beta y@f$. Fortunately deal.ii already has
 * this operation built in to its vector class, so we get the
 * reference to the vector y and call the sadd method.
 * 
 * @code
 *   int my_Sum(braid_App app,
 *              double alpha,
 *              braid_Vector x,
 *              double beta,
 *              braid_Vector y)
 *   {
 *     UNUSED(app);
 *     Vector<double>& vec = y->data;
 *     vec.sadd(beta, alpha, x->data);
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This calculates the spatial norm using the l2 norm. According
 * to XBraid, this could be just about any spatial norm but we'll
 * keep it simple and used deal.ii vector's built in l2_norm method.
 * 
 * @code
 *   int
 *   my_SpatialNorm(braid_App     app,
 *                  braid_Vector  u,
 *                  double       *norm_ptr)
 *   {
 *     UNUSED(app);
 *     double dot = 0.0;
 *     dot = u->data.l2_norm();
 *     *norm_ptr = dot;
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This function is called at various points depending on the access
 * level specified when configuring the XBraid struct. This function
 * is used to print out data during the run time, such as plots of the
 * data. The status struct contains a ton of information about the
 * simulation run. Here we get the current time and timestep number.
 * The output_results function is called to plot the solution data.
 * If the method of manufactured solutions is being used, then the
 * error of this time step is computed and processed.
 * 
 * @code
 *   int
 *   my_Access(braid_App          app,
 *             braid_Vector       u,
 *             braid_AccessStatus astatus)
 *   {
 *     double     t;
 *     int index;
 *   
 *     braid_AccessStatusGetT(astatus, &t);
 *     braid_AccessStatusGetTIndex(astatus, &index);
 *   
 *     app->eq.output_results(index, t, u->data);
 *   
 *   #if DO_MFG
 *     if(index == app->final_step)
 *       {
 *         app->eq.process_solution(t, index, u->data);
 *       }
 *   #endif
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This calculates the size of buffer needed to pack the solution
 * data into a linear buffer for transfer to another processor via
 * MPI. We query the size of the data from the HeatEquation class
 * and return the buffer size.
 * 
 * @code
 *   int
 *   my_BufSize(braid_App           app,
 *              int                 *size_ptr,
 *              braid_BufferStatus  bstatus)
 *   {
 *     UNUSED(bstatus);
 *     int size = app->eq.size();
 *     *size_ptr = (size+1)*sizeof(double);
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This function packs a linear buffer with data so that the buffer
 * may be sent to another processor via MPI. The buffer is cast to
 * a type we can work with. The first element of the buffer is the
 * size of the buffer. Then we iterate over soltuion vector u and
 * fill the buffer with our solution data. Finally we tell XBraid
 * how much data we wrote.
 * 
 * @code
 *   int
 *   my_BufPack(braid_App           app,
 *              braid_Vector        u,
 *              void               *buffer,
 *              braid_BufferStatus  bstatus)
 *   {
 *   
 *     UNUSED(app);
 *     double *dbuffer = (double*)buffer;
 *     int size = u->data.size();
 *     dbuffer[0] = size;
 *     for(int i=0; i != size; ++i)
 *       {
 *         dbuffer[i+1] = (u->data)[i];
 *       }
 *     braid_BufferStatusSetSize(bstatus, (size+1)*sizeof(double));
 *   
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * This function unpacks a buffer that was recieved from a different
 * processor via MPI. The size of the buffer is read from the first
 * element, then we iterate over the size of the buffer and fill
 * the values of solution vector u with the data in the buffer.
 * 
 * @code
 *   int
 *   my_BufUnpack(braid_App           app,
 *                void               *buffer,
 *                braid_Vector       *u_ptr,
 *                braid_BufferStatus  bstatus)
 *   {
 *     UNUSED(app);
 *     UNUSED(bstatus);
 *   
 *     my_Vector *u = NULL;
 *     double *dbuffer = (double*)buffer;
 *     int size = static_cast<int>(dbuffer[0]);
 *     u = new(my_Vector);
 *     u->data.reinit(size);
 *   
 *     for(int i = 0; i != size; ++i)
 *       {
 *         (u->data)[i] = dbuffer[i+1];
 *       }
 *     *u_ptr = u;
 *   
 *     return 0;
 *   }
 * @endcode
 
 
<a name="ann-src/BraidFuncs.hh"></a>
<h1>Annotated version of src/BraidFuncs.hh</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #ifndef _BRAIDFUNCS_H_
 *   #define _BRAIDFUNCS_H_
 *   
 *   /**
 *    *   \file BraidFuncs.cc
 *    *   \brief Contains the implementation of the mandatory X-Braid functions
 *    *
 *    *  X-Braid mandates several functions in order to drive the solution.
 *    *  This file contains the implementation of said mandatory functions.
 *    *  See the X-Braid documentation for more information.
 *    *  There are several functions that are optional in X-Braid that may
 *    *  or may not be implemented in here.
 *    *
 *    */
 *   
 *   
 *   /*-------- Third Party --------*/
 *   #include <deal.II/numerics/vector_tools.h>
 *   
 *   #include <braid.h>
 *   #include <braid_test.h>
 *   
 *   /*-------- Project --------*/
 *   #include "HeatEquation.hh"
 *   
 * @endcode
 * 
 * This struct contains all data that changes with time. For now
 * this is just the solution data. When doing AMR this should
 * probably include the triangulization, the sparsity patter,
 * constraints, etc.
 * 
 * @code
 *   /**
 *    * \brief Struct that contains the deal.ii vector.
 *    */
 *   typedef struct _braid_Vector_struct
 *   {
 *     ::Vector<double> data;
 *   } my_Vector;
 *   
 * @endcode
 * 
 * This struct contains all the data that is unchanging with time.
 * 
 * @code
 *   /**
 *    * \brief Struct that contains the HeatEquation and final
 *    * time step number.
 *    */
 *   typedef struct _braid_App_struct
 *   {
 *     HeatEquation<2> eq;
 *     int final_step;
 *   } my_App;
 *   
 *   
 *   /**
 *    * @brief my_Step - Takes a step in time, advancing the u vector
 *    *
 *    * @param app - The braid app struct
 *    * @param ustop - The solution data at the end of this time step
 *    * @param fstop - RHS data (such as forcing function?)
 *    * @param u - The solution data at the beginning of this time step
 *    * @param status - Status structure that contains various info of this time
 *    *
 *    * @return Success (0) or failure (1)
 *    **/
 *   int my_Step(braid_App        app,
 *               braid_Vector     ustop,
 *               braid_Vector     fstop,
 *               braid_Vector     u,
 *               braid_StepStatus status);
 *   
 *   
 *   /**
 *    * @brief my_Init - Initializes a solution data at the given time
 *    * For now, initializes the solution to zero no matter what time we are at
 *    *
 *    * @param app - The braid app struct containing user data
 *    * @param t - Time at which the solution is initialized
 *    * @param u_ptr - The solution data that needs to be filled
 *    *
 *    * @return Success (0) or failure (1)
 *    **/
 *   int
 *   my_Init(braid_App     app,
 *           double        t,
 *           braid_Vector *u_ptr);
 *   
 *   
 *   /**
 *    * @brief my_Clone - Clones a vector into a new vector
 *    *
 *    * @param app - The braid app struct containing user data
 *    * @param u - The existing vector containing data
 *    * @param v_ptr - The empty vector that needs to be filled
 *    *
 *    * @return Success (0) or failure (1)
 *    **/
 *   int
 *   my_Clone(braid_App     app,
 *            braid_Vector  u,
 *            braid_Vector *v_ptr);
 *   
 *   
 *   /**
 *    * @brief my_Free - Deletes a vector
 *    *
 *    * @param app - The braid app struct containing user data
 *    * @param u - The vector that needs to be deleted
 *    *
 *    * @return Success (0) or failure (1)
 *    **/
 *   int
 *   my_Free(braid_App    app,
 *           braid_Vector u);
 *   
 *   
 *   /**
 *    * @brief my_Sum - Sums two vectors in an AXPY operation
 *    * The operation is y = alpha*x + beta*y
 *    *
 *    * @param app - The braid app struct containing user data
 *    * @param alpha - The coefficient in front of x
 *    * @param x - A vector that is multiplied by alpha then added to y
 *    * @param beta - The coefficient of y
 *    * @param y - A vector that is multiplied by beta then summed with x
 *    *
 *    * @return Success (0) or failure (1)
 *    **/
 *   int
 *   my_Sum(braid_App app,
 *          double alpha,
 *          braid_Vector x,
 *          double beta,
 *          braid_Vector y);
 *   
 *   /**
 *    *  \brief Returns the spatial norm of the provided vector
 *    *
 *    *  Calculates and returns the spatial norm of the provided vector.
 *    *  Interestingly enough, X-Braid does not specify a particular norm.
 *    *  to keep things simple, we implement the Euclidean norm.
 *    *
 *    *  \param app - The braid app struct containing user data
 *    *  \param u - The vector we need to take the norm of
 *    *  \param norm_ptr - Pointer to the norm that was calculated, need to modify this
 *    *  \return Success (0) or failure (1)
 *    */
 *   int
 *   my_SpatialNorm(braid_App     app,
 *                  braid_Vector  u,
 *                  double       *norm_ptr);
 *   
 *   /**
 *    *  \brief Allows the user to output details
 *    *
 *    *  The Access function is called at various points to allow the user to output
 *    *  information to the screen or to files.
 *    *  The astatus parameter provides various information about the simulation,
 *    *  see the XBraid documentation for details on what information you can get.
 *    *  Example information is what the current timestep number and current time is.
 *    *  If the access level (in parallel_in_time.cc) is set to 0, this function is
 *    *  never called.
 *    *  If the access level is set to 1, the function is called after the last
 *    *  XBraid cycle.
 *    *  If the access level is set to 2, it is called every XBraid cycle.
 *    *
 *    *  \param app - The braid app struct containing user data
 *    *  \param u - The vector containing the data at the status provided
 *    *  \param astatus - The Braid status structure
 *    *  \return Success (0) or failure (1)
 *    */
 *   int
 *   my_Access(braid_App          app,
 *             braid_Vector       u,
 *             braid_AccessStatus astatus);
 *   
 *   /**
 *    *  \brief Calculates the size of a buffer for MPI data transfer
 *    *
 *    *  Calculates the size of the buffer that is needed to transfer
 *    *  a solution vector to another processor.
 *    *  The bstatus parameter provides various information on the
 *    *  simulation, see the XBraid documentation for all possible
 *    *  fields.
 *    *
 *    *  \param app - The braid app struct containing user data
 *    *  \param size_ptr A pointer to the calculated size
 *    *  \param bstatus The XBraid status structure
 *    *  \return Success (0) or failure (1)
 *    */
 *   int
 *   my_BufSize(braid_App           app,
 *              int                 *size_ptr,
 *              braid_BufferStatus  bstatus);
 *   
 *   /**
 *    *  \brief Linearizes a vector to be sent to another processor
 *    *
 *    *  Linearizes (packs) a data buffer with the contents of
 *    *  some solution state u.
 *    *
 *    *  \param app - The braid app struct containing user data
 *    *  \param u The vector that must be packed into buffer
 *    *  \param buffer The buffer that must be filled with u
 *    *  \param bstatus The XBraid status structure
 *    *  \return Success (0) or failure (1)
 *    */
 *   int
 *   my_BufPack(braid_App           app,
 *              braid_Vector        u,
 *              void               *buffer,
 *              braid_BufferStatus  bstatus);
 *   
 *   /**
 *    *  \brief Unpacks a vector that was sent from another processor
 *    *
 *    *  Unpacks a linear data buffer into the vector pointed to by
 *    *  u_ptr.
 *    *
 *    *  \param app - The braid app struct containing user data
 *    *  \param buffer The buffer that must be unpacked
 *    *  \param u_ptr The pointer to the vector that is filled
 *    *  \param bstatus The XBraid status structure
 *    *  \return Success (0) or failure (1)
 *    */
 *   int
 *   my_BufUnpack(braid_App           app,
 *                void               *buffer,
 *                braid_Vector       *u_ptr,
 *                braid_BufferStatus  bstatus);
 *   
 *   #endif // _BRAIDFUNCS_H_
 * @endcode
 
 
<a name="ann-src/HeatEquation.hh"></a>
<h1>Annotated version of src/HeatEquation.hh</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #ifndef _HEATEQUATION_H_
 *   #define _HEATEQUATION_H_
 *   
 *   #include <deal.II/base/utilities.h>
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/base/logstream.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/solver_cg.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/affine_constraints.h>
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/grid/grid_refinement.h>
 *   #include <deal.II/grid/grid_out.h>
 *   #include <deal.II/grid/tria_accessor.h>
 *   #include <deal.II/grid/tria_iterator.h>
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/dofs/dof_accessor.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/numerics/data_out.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/numerics/error_estimator.h>
 *   #include <deal.II/numerics/solution_transfer.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/base/convergence_table.h>
 *   
 *   #include <fstream>
 *   
 *   using namespace dealii;
 *   
 * @endcode
 * 
 * The HeatEquation class is describes the finite element
 * solver for the heat equation. It contains all the functions
 * needed to define the problem domain and advance the solution
 * in time.
 * 
 * @code
 *   template <int dim>
 *   class HeatEquation
 *   {
 *   public:
 *     HeatEquation();
 *     void define();
 *     void step(Vector<double>& braid_data,
 *               double deltaT,
 *               double a_time,
 *               int a_time_idx);
 *   
 *     int size() const; /// Returns the size of the solution vector
 *   
 *     void output_results(int a_time_idx,
 *                         double a_time,
 *                         Vector<double>& a_solution) const;
 *   
 *     void initialize(double a_time,
 *                     Vector<double>& a_vector) const;
 *   
 *     void process_solution(double a_time,
 *                           int a_index,
 *                           const Vector<double>& a_vector);
 *   
 *   private:
 *     void setup_system();
 *     void solve_time_step(Vector<double>& a_solution);
 *   
 *     Triangulation<dim>   triangulation;
 *     FE_Q<dim>            fe;
 *     DoFHandler<dim>      dof_handler;
 *   
 *     AffineConstraints<double> constraints;
 *   
 *     SparsityPattern      sparsity_pattern;
 *     SparseMatrix<double> mass_matrix;
 *     SparseMatrix<double> laplace_matrix;
 *     SparseMatrix<double> system_matrix;
 *   
 *     Vector<double>       system_rhs;
 *   
 *     std::ofstream        myfile;
 *   
 *     const double         theta;
 *   
 * @endcode
 * 
 * These were originally in the run() function but because
 * I am splitting the run() function up into define and step
 * they need to become member data
 * 
 * @code
 *     Vector<double> tmp;
 *     Vector<double> forcing_terms;
 *   
 *     ConvergenceTable convergence_table;
 *   };
 *   
 * @endcode
 * 
 * The RightHandSide class describes the RHS of the governing
 * equations. In this case, it is the forcing function.
 * 
 * @code
 *   template <int dim>
 *   class RightHandSide : public Function<dim>
 *   {
 *   public:
 *     RightHandSide ()
 *       :
 *       Function<dim>(),
 *       period (0.2)
 *     {}
 *   
 *     virtual double value (const Point<dim> &p,
 *                           const unsigned int component = 0) const override;
 *   
 *   private:
 *     const double period;
 *   };
 *   
 * @endcode
 * 
 * The BoundaryValues class describes the boundary conditions
 * of the governing equations.
 * 
 * @code
 *   template <int dim>
 *   class BoundaryValues : public Function<dim>
 *   {
 *   public:
 *     virtual double value (const Point<dim>  &p,
 *                           const unsigned int component = 0) const override;
 *   };
 *   
 * @endcode
 * 
 * The RightHandSideMFG class describes the right hand side
 * function when doing the method of manufactured solutions.
 * 
 * @code
 *   template <int dim>
 *   class RightHandSideMFG : public Function<dim>
 *   {
 *   public:
 *     virtual double value (const Point<dim>  &p,
 *                           const unsigned int component = 0) const override;
 *   };
 *   
 * @endcode
 * 
 * The InitialValuesMFG class describes the initial values
 * when doing the method of manufactured solutions.
 * 
 * @code
 *   template <int dim>
 *   class InitialValuesMFG : public Function<dim>
 *   {
 *   public:
 *     virtual double value (const Point<dim> &p,
 *                           const unsigned int component = 0) const override;
 *   };
 *   
 * @endcode
 * 
 * Provides the exact value for the manufactured solution. This
 * is used for the boundary conditions as well.
 * 
 * @code
 *   template <int dim>
 *   class ExactValuesMFG : public Function<dim>
 *   {
 *   public:
 *     /**
 *      *  \brief Computes the value at the given point and member data time
 *      *
 *      *  Computes the exact value of the manufactured solution at point p and
 *      *  the member data time. See the class documentation and the design doc
 *      *  for details on what the exact solution is.
 *      *
 *      *  \param p The point that the exact solution is computed at
 *      *  \param component The component of the exact solution (always 0 for now)
 *      *  \return double The exact value that was computed
 *      */
 *     virtual double value (const Point<dim> &p,
 *                           const unsigned int component = 0) const override;
 *   
 *     /**
 *      *  \brief Computes the gradient of the exact solution at the given point
 *      *
 *      *  Computes the gradient of the exact/manufactured solution value at
 *      *  point p and member data time. See the design doc for details on
 *      *  what the gradient of the exact solution is
 *      *
 *      *  \param p The point that the gradient is calculated at
 *      *  \param component The component of the system of equations this gradient is for
 *      *  \return Tensor<1,dim> A rank 1 tensor that contains the gradient
 *      *  in each spatial dimension
 *      */
 *     virtual Tensor<1,dim> gradient (const Point<dim>   &p,
 *                                     const unsigned int  component = 0) const override;
 *   };
 *   
 *   
 *   #include "HeatEquationImplem.hh"
 *   
 *   #endif // _HEATEQUATION_H_
 * @endcode
 
 
<a name="ann-src/HeatEquationImplem.hh"></a>
<h1>Annotated version of src/HeatEquationImplem.hh</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #include "Utilities.hh"
 *   
 *   #include <iomanip>
 *   #include <math.h>
 *   
 * @endcode
 * 
 * Calculates the forcing function for the RightHandSide. See the
 * documentation for the math.
 * 
 * @code
 *   template <int dim>
 *   double RightHandSide<dim>::value (const Point<dim> &p,
 *                                     const unsigned int component) const
 *   {
 *     (void) component;
 *     Assert (component == 0, ExcIndexRange(component, 0, 1));
 *     Assert (dim == 2, ExcNotImplemented());
 *   
 *     double time = this->get_time();
 *   
 *     if ((p[0] > 0.5) && (p[1] > -0.5))
 *       {
 *         return std::exp(-0.5*(time-0.125)*(time-0.125)/(0.005));
 *       }
 *     else if ((p[0] > -0.5) && (p[1] > 0.5))
 *       {
 *         return std::exp(-0.5*(time-0.375)*(time-0.375)/(0.005));
 *       }
 *     else
 *       {
 *         return 0;
 *       }
 *   
 *     return 0; // No forcing function
 *   }
 *   
 * @endcode
 * 
 * Calculates the forcing function for the method of manufactured
 * solutions. See the documentation for the math.
 * 
 * @code
 *   template <int dim>
 *   double RightHandSideMFG<dim>::value (const Point<dim> &p,
 *                                        const unsigned int component) const
 *   {
 *     (void) component;
 *     Assert (component == 0, ExcIndexRange(component, 0, 1));
 *     Assert (dim == 2, ExcNotImplemented());
 *   
 *     double time = this->get_time();
 *   
 *     double pi = numbers::PI;
 *     return 4*pi*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
 *   }
 *   
 * @endcode
 * 
 * Calculates the boundary conditions, essentially zero everywhere.
 * 
 * @code
 *   template <int dim>
 *   double BoundaryValues<dim>::value (const Point<dim> &p,
 *                               const unsigned int component) const
 *   {
 *     UNUSED(p);
 *     (void) component;
 *     Assert (component == 0, ExcIndexRange(component, 0, 1));
 *     return 0;
 *   }
 *   
 * @endcode
 * 
 * Calculates the exact solution (and thus also boundary conditions)
 * for the method of manufactured solutions.
 * 
 * @code
 *   template <int dim>
 *   double ExactValuesMFG<dim>::value (const Point<dim> &p,
 *                                      const unsigned int component) const
 *   {
 *     (void) component;
 *     Assert (component == 0, ExcIndexRange(component, 0, 1));
 *   
 *     double time = this->get_time();
 *     const double pi = numbers::PI;
 *   
 *     return std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
 *   }
 *   
 * @endcode
 * 
 * Calculates the gradient of the exact solution for the method of manufactured
 * solutions. See the documentation for the math.
 * 
 * @code
 *   template <int dim>
 *   Tensor<1,dim> ExactValuesMFG<dim>::gradient (const Point<dim>   &p,
 *                                                const unsigned int) const
 *   {
 *     Assert (dim == 2, ExcNotImplemented());
 *   
 *     Tensor<1,dim> return_value;
 *     const double pi = numbers::PI;
 *     double time = this->get_time();
 *     return_value[0] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[1])*std::sin(2*pi*p[0]);
 *     return_value[1] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::sin(2*pi*p[1]);
 *     return return_value;
 *   }
 *   
 * @endcode
 * 
 * Calculates the initial values for the method of manufactured solutions.
 * See the documentation for the math.
 * 
 * @code
 *   template <int dim>
 *   double InitialValuesMFG<dim>::value (const Point<dim> &p,
 *                                        const unsigned int component) const
 *   {
 *     (void) component;
 *     Assert (component == 0, ExcIndexRange(component, 0, 1));
 *     const double pi = numbers::PI;
 *   
 *     return std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
 *   }
 *   
 *   template <int dim>
 *   HeatEquation<dim>::HeatEquation ()
 *     :
 *     fe(1),
 *     dof_handler(triangulation),
 *     theta(0.5)
 *   {
 *   }
 *   
 *   template <int dim>
 *   void HeatEquation<dim>::initialize(double a_time,
 *                                      Vector<double>& a_vector) const
 *   {
 *   #if DO_MFG
 * @endcode
 * 
 * We only initialize values in the manufactured solution case
 * 
 * @code
 *     InitialValuesMFG<dim> iv_function;
 *     iv_function.set_time(a_time);
 *     VectorTools::project (dof_handler, constraints,
 *                           QGauss<dim>(fe.degree+1), iv_function,
 *                           a_vector);
 *   #else
 *     UNUSED(a_time);
 *     UNUSED(a_vector);
 *   #endif // DO_MFG
 * @endcode
 * 
 * If not the MFG solution case, a_vector is already zero'd so do nothing
 * 
 * @code
 *   }
 *   
 *   template <int dim>
 *   void HeatEquation<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 *   
 *     constraints.clear ();
 *     DoFTools::make_hanging_node_constraints (dof_handler,
 *                                              constraints);
 *     constraints.close();
 *   
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler,
 *                                     dsp,
 *                                     constraints,
 *                                     /*keep_constrained_dofs = */ true);
 *     sparsity_pattern.copy_from(dsp);
 *   
 *     mass_matrix.reinit(sparsity_pattern);
 *     laplace_matrix.reinit(sparsity_pattern);
 *     system_matrix.reinit(sparsity_pattern);
 *   
 *     MatrixCreator::create_mass_matrix(dof_handler,
 *                                       QGauss<dim>(fe.degree+1),
 *                                       mass_matrix);
 *     MatrixCreator::create_laplace_matrix(dof_handler,
 *                                          QGauss<dim>(fe.degree+1),
 *                                          laplace_matrix);
 *   
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   }
 *   
 *   
 *   template <int dim>
 *   void HeatEquation<dim>::solve_time_step(Vector<double>& a_solution)
 *   {
 *     SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
 *     SolverCG<> cg(solver_control);
 *   
 *     PreconditionSSOR<> preconditioner;
 *     preconditioner.initialize(system_matrix, 1.0);
 *   
 *     cg.solve(system_matrix, a_solution, system_rhs,
 *              preconditioner);
 *   
 *     constraints.distribute(a_solution);
 *   }
 *   
 *   
 *   
 *   template <int dim>
 *   void HeatEquation<dim>::output_results(int a_time_idx,
 *                                          double a_time,
 *                                          Vector<double>& a_solution) const
 *   {
 *   
 *     DataOutBase::VtkFlags vtk_flags;
 *     vtk_flags.time = a_time;
 *     vtk_flags.cycle = a_time_idx;
 *   
 *     DataOut<dim> data_out;
 *     data_out.set_flags(vtk_flags);
 *   
 *     data_out.attach_dof_handler(dof_handler);
 *     data_out.add_data_vector(a_solution, "U");
 *   
 *     data_out.build_patches();
 *   
 *     const std::string filename = "solution-"
 *       + Utilities::int_to_string(a_time_idx, 3) +
 *       ".vtk";
 *     std::ofstream output(filename.c_str());
 *     data_out.write_vtk(output);
 *   }
 *   
 * @endcode
 * 
 * We define the geometry here, this is called on each processor
 * and doesn't change in time. Once doing AMR, this won't need
 * to exist anymore.
 * 
 * @code
 *   template <int dim>
 *   void HeatEquation<dim>::define()
 *   {
 *     const unsigned int initial_global_refinement = 6;
 *   
 *     GridGenerator::hyper_L (triangulation);
 *     triangulation.refine_global (initial_global_refinement);
 *   
 *     setup_system();
 *   
 *     tmp.reinit (dof_handler.n_dofs());
 *     forcing_terms.reinit (dof_handler.n_dofs());
 *   }
 *   
 * @endcode
 * 
 * Here we advance the solution forward in time. This is done
 * the same way as in the loop in @ref step_26 "step-26"'s run function.
 * 
 * @code
 *   template<int dim>
 *   void HeatEquation<dim>::step(Vector<double>& braid_data,
 *                                double deltaT,
 *                                double a_time,
 *                                int a_time_idx)
 *   {
 *     a_time += deltaT;
 *     ++a_time_idx;
 *   
 *     mass_matrix.vmult(system_rhs, braid_data);
 *   
 *     laplace_matrix.vmult(tmp, braid_data);
 *   
 *     system_rhs.add(-(1 - theta) * deltaT, tmp);
 *   
 *   #if DO_MFG
 *     RightHandSideMFG<dim> rhs_function;
 *   #else
 *     RightHandSide<dim> rhs_function;
 *   #endif
 *     rhs_function.set_time(a_time);
 *     VectorTools::create_right_hand_side(dof_handler,
 *                                         QGauss<dim>(fe.degree+1),
 *                                         rhs_function,
 *                                         tmp);
 *   
 *     forcing_terms = tmp;
 *     forcing_terms *= deltaT * theta;
 *   
 *     rhs_function.set_time(a_time - deltaT);
 *     VectorTools::create_right_hand_side(dof_handler,
 *                                         QGauss<dim>(fe.degree+1),
 *                                         rhs_function,
 *                                         tmp);
 *   
 *     forcing_terms.add(deltaT * (1 - theta), tmp);
 *     system_rhs += forcing_terms;
 *   
 *     system_matrix.copy_from(mass_matrix);
 *     system_matrix.add(theta * deltaT, laplace_matrix);
 *   
 *     constraints.condense (system_matrix, system_rhs);
 *   
 *     {
 *   #if DO_MFG
 * @endcode
 * 
 * If we are doing the method of manufactured solutions
 * then we set the boundary conditions to the exact solution.
 * Otherwise the boundary conditions are zero.
 * 
 * @code
 *       ExactValuesMFG<dim> boundary_values_function;
 *   #else
 *       BoundaryValues<dim> boundary_values_function;
 *   #endif
 *       boundary_values_function.set_time(a_time);
 *   
 *       std::map<types::global_dof_index, double> boundary_values;
 *       VectorTools::interpolate_boundary_values(dof_handler,
 *                                                0,
 *                                                boundary_values_function,
 *                                                boundary_values);
 *   
 *       MatrixTools::apply_boundary_values(boundary_values,
 *                                          system_matrix,
 *                                          braid_data,
 *                                          system_rhs);
 *     }
 *   
 *     solve_time_step(braid_data);
 *   }
 *   
 *   template<int dim>
 *   int HeatEquation<dim>::size() const
 *   {
 *     return dof_handler.n_dofs();
 *   }
 *   
 * @endcode
 * 
 * This function computes the error for the time step when doing
 * the method of manufactured solutions. First the exact values
 * is calculated, then the difference per cell is computed for
 * the various norms, and the error is computed. This is written
 * out to a pretty table.
 * 
 * @code
 *   template<int dim> void
 *   HeatEquation<dim>::process_solution(double a_time,
 *                                       int a_index,
 *                                       const Vector<double>& a_vector)
 *   {
 * @endcode
 * 
 * Compute the exact value for the manufactured solution case
 * 
 * @code
 *     ExactValuesMFG<dim> exact_function;
 *     exact_function.set_time(a_time);
 *   
 *     Vector<double> difference_per_cell (triangulation.n_active_cells());
 *     VectorTools::integrate_difference(dof_handler,
 *                                       a_vector,
 *                                       exact_function,
 *                                       difference_per_cell,
 *                                       QGauss<dim>(fe.degree+1),
 *                                       VectorTools::L2_norm);
 *   
 *     const double L2_error = VectorTools::compute_global_error(triangulation,
 *                                                               difference_per_cell,
 *                                                               VectorTools::L2_norm);
 *   
 *     VectorTools::integrate_difference(dof_handler,
 *                                       a_vector,
 *                                       exact_function,
 *                                       difference_per_cell,
 *                                       QGauss<dim>(fe.degree+1),
 *                                       VectorTools::H1_seminorm);
 *   
 *     const double H1_error = VectorTools::compute_global_error(triangulation,
 *                                                               difference_per_cell,
 *                                                               VectorTools::H1_seminorm);
 *   
 *     const QTrapezoid<1> q_trapez;
 *     const QIterated<dim> q_iterated (q_trapez, 5);
 *     VectorTools::integrate_difference (dof_handler,
 *                                        a_vector,
 *                                        exact_function,
 *                                        difference_per_cell,
 *                                        q_iterated,
 *                                        VectorTools::Linfty_norm);
 *     const double Linfty_error = VectorTools::compute_global_error(triangulation,
 *                                                                   difference_per_cell,
 *                                                                   VectorTools::Linfty_norm);
 *   
 *     const unsigned int n_active_cells = triangulation.n_active_cells();
 *     const unsigned int n_dofs = dof_handler.n_dofs();
 *   
 *     pout() << "Cycle " << a_index << ':'
 *            << std::endl
 *            << "   Number of active cells:       "
 *            << n_active_cells
 *            << std::endl
 *            << "   Number of degrees of freedom: "
 *            << n_dofs
 *            << std::endl;
 *   
 *     convergence_table.add_value("cycle", a_index);
 *     convergence_table.add_value("cells", n_active_cells);
 *     convergence_table.add_value("dofs", n_dofs);
 *     convergence_table.add_value("L2", L2_error);
 *     convergence_table.add_value("H1", H1_error);
 *     convergence_table.add_value("Linfty", Linfty_error);
 *   
 *     convergence_table.set_precision("L2", 3);
 *     convergence_table.set_precision("H1", 3);
 *     convergence_table.set_precision("Linfty", 3);
 *   
 *     convergence_table.set_scientific("L2", true);
 *     convergence_table.set_scientific("H1", true);
 *     convergence_table.set_scientific("Linfty", true);
 *   
 *     convergence_table.set_tex_caption("cells", "\\# cells");
 *     convergence_table.set_tex_caption("dofs", "\\# dofs");
 *     convergence_table.set_tex_caption("L2", "@fL^2@f-error");
 *     convergence_table.set_tex_caption("H1", "@fH^1@f-error");
 *     convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
 *   
 *     convergence_table.set_tex_format("cells", "r");
 *     convergence_table.set_tex_format("dofs", "r");
 *   
 *     std::cout << std::endl;
 *     convergence_table.write_text(std::cout);
 *   
 *     std::ofstream error_table_file("tex-conv-table.tex");
 *     convergence_table.write_tex(error_table_file);
 *   }
 * @endcode
 
 
<a name="ann-src/Utilities.cc"></a>
<h1>Annotated version of src/Utilities.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #include "Utilities.hh"
 *   
 *   #include <string>
 *   #include <fstream>
 *   
 *   #include <mpi.h>
 *   
 *   int procID = 0;
 *   
 * @endcode
 * 
 * The shared variables
 * 
 
 * 
 * 
 * @code
 *   static std::string       s_pout_filename ;
 *   static std::string       s_pout_basename ;
 *   static std::ofstream     s_pout ;
 *   
 *   static bool              s_pout_init = false ;
 *   static bool              s_pout_open = false ;
 *   
 *   #ifdef USE_MPI
 * @endcode
 * 
 * in parallel, compute the filename give the basename
 * [NOTE: dont call this before MPI is initialized.]
 * 
 * @code
 *   static void setFileName()
 *   {
 *     static const size_t ProcnumSize = 1 + 10 + 1 ;  //'.' + 10digits + '\0'
 *     char procnum[ProcnumSize] ;
 *     snprintf( procnum ,ProcnumSize ,".%d" ,procID);
 *               s_pout_filename = s_pout_basename + procnum ;
 *   }
 *   
 * @endcode
 * 
 * in parallel, close the file if nec., open it and check for success
 * 
 * @code
 *   static void openFile()
 *   {
 *     if ( s_pout_open )
 *     {
 *       s_pout.close();
 *     }
 *     s_pout.open( s_pout_filename.c_str() );
 * @endcode
 * 
 * if open() fails, we have problems, but it's better
 * to try again later than to make believe it succeeded
 * 
 * @code
 *     s_pout_open = (bool)s_pout ;
 *   }
 *   
 *   #else
 * @endcode
 * 
 * in serial, filename is always cout
 * 
 * @code
 *   static void setFileName()
 *   {
 *     s_pout_filename = "cout" ;
 *   }
 *   
 * @endcode
 * 
 * in serial, this does absolutely nothing
 * 
 * @code
 *   static void openFile()
 *   {
 *   }
 *   #endif
 *   
 *   std::ostream& pout()
 *   {
 *   #ifdef USE_MPI
 * @endcode
 * 
 * the common case is _open == true, which just returns s_pout
 * 
 * @code
 *     if ( ! s_pout_open )
 *       {
 * @endcode
 * 
 * the uncommon cae: the file isn't opened, MPI may not be
 * initialized, and the basename may not have been set
 * 
 * @code
 *         int flag_i, flag_f;
 *         MPI_Initialized(&flag_i);
 *         MPI_Finalized(&flag_f);
 * @endcode
 * 
 * app hasn't set a basename yet, so set the default
 * 
 * @code
 *         if ( ! s_pout_init )
 *           {
 *             s_pout_basename = "pout" ;
 *             s_pout_init = true ;
 *           }
 * @endcode
 * 
 * if MPI not initialized, we cant open the file so return cout
 * 
 * @code
 *         if ( ! flag_i || flag_f)
 *           {
 *             return std::cout; // MPI hasn't been started yet, or has ended....
 *           }
 * @endcode
 * 
 * MPI is initialized, so file must not be, so open it
 * 
 * @code
 *         setFileName() ;
 *         openFile() ;
 * @endcode
 * 
 * finally, in case the open failed, return cout
 * 
 * @code
 *         if ( ! s_pout_open )
 *           {
 *             return std::cout ;
 *           }
 *       }
 *     return s_pout ;
 *   #else
 *     return std::cout;
 *   #endif
 *   }
 * @endcode
 
 
<a name="ann-src/Utilities.hh"></a>
<h1>Annotated version of src/Utilities.hh</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #ifndef _UTILITIES_H_
 *   #define _UTILITIES_H_
 *   
 *   #include <iostream>
 *   
 * @endcode
 * 
 * This preprocessor macro is used on function arguments
 * that are not used in the function. It is used to
 * suppress compiler warnings.
 * 
 * @code
 *   #define UNUSED(x) (void)(x)
 *   
 * @endcode
 * 
 * Contains the current MPI processor ID.
 * 
 * @code
 *   extern int procID;
 *   
 * @endcode
 * 
 * Function to return the ostream to write out to. In MPI
 * mode it returns a stream to a file named pout.<#> where
 * <#> is the procID. This allows the user to write output
 * from each processor to a separate file. In serial mode
 * (no MPI), it returns the standard output.
 * 
 * @code
 *   std::ostream& pout();
 *   #endif
 * @endcode
 
 
<a name="ann-src/parallel_in_time.cc"></a>
<h1>Annotated version of src/parallel_in_time.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* ---------------------------------------------------------------------
 *    *
 *    * Copyright (C) 2013 - 2018 by the deal.II authors
 *    *
 *    * This file is part of the deal.II library.
 *    *
 *    * The deal.II library is free software; you can use it, redistribute
 *    * it, and/or modify it under the terms of the GNU Lesser General
 *    * Public License as published by the Free Software Foundation; either
 *    * version 2.1 of the License, or (at your option) any later version.
 *    * The full text of the license can be found in the file LICENSE at
 *    * the top level of the deal.II distribution.
 *    *
 *    * ---------------------------------------------------------------------
 *   
 *    *
 *    * Author: Joshua Christopher, Colorado State University, 2018
 *    */
 *   
 *   #include "BraidFuncs.hh"
 *   #include "HeatEquation.hh"
 *   #include "Utilities.hh"
 *   
 *   #include <fstream>
 *   #include <iostream>
 *   
 *   int main(int argc, char *argv[])
 *   {
 *     try
 *       {
 *         using namespace dealii;
 *   
 *         /* Initialize MPI */
 *         MPI_Comm      comm; //, comm_x, comm_t;
 *         int rank;
 *         MPI_Init(&argc, &argv);
 *         comm   = MPI_COMM_WORLD;
 *         MPI_Comm_rank(comm, &rank);
 *         procID = rank;
 *   
 * @endcode
 * 
 * Set up X-Braid
 * 
 * @code
 *         /* Initialize Braid */
 *         braid_Core core;
 *         double tstart = 0.0;
 *         double tstop = 0.002;
 *         int    ntime = 10;
 *         my_App *app = new(my_App);
 *   
 *         braid_Init(MPI_COMM_WORLD, comm, tstart, tstop, ntime, app,
 *                    my_Step, my_Init, my_Clone, my_Free, my_Sum, my_SpatialNorm,
 *                    my_Access, my_BufSize, my_BufPack, my_BufUnpack, &core);
 *   
 *         /* Define XBraid parameters
 *          * See -help message forf descriptions */
 *         int       max_levels    = 3;
 * @endcode
 * 
 * int       nrelax        = 1;
 * int       skip          = 0;
 * 
 * @code
 *         double    tol           = 1.e-7;
 * @endcode
 * 
 * int       cfactor       = 2;
 * 
 * @code
 *         int       max_iter      = 5;
 * @endcode
 * 
 * int       min_coarse    = 10;
 * int       fmg           = 0;
 * int       scoarsen      = 0;
 * int       res           = 0;
 * int       wrapper_tests = 0;
 * 
 * @code
 *         int       print_level   = 1;
 *         int       access_level  = 1;
 *         int       use_sequential= 0;
 *   
 *         braid_SetPrintLevel( core, print_level);
 *         braid_SetAccessLevel( core, access_level);
 *         braid_SetMaxLevels(core, max_levels);
 * @endcode
 * 
 * braid_SetMinCoarse( core, min_coarse );
 * braid_SetSkip(core, skip);
 * braid_SetNRelax(core, -1, nrelax);
 * 
 * @code
 *         braid_SetAbsTol(core, tol);
 * @endcode
 * 
 * braid_SetCFactor(core, -1, cfactor);
 * 
 * @code
 *         braid_SetMaxIter(core, max_iter);
 *         braid_SetSeqSoln(core, use_sequential);
 *   
 *         app->eq.define();
 *         app->final_step = ntime;
 *   
 *         braid_Drive(core);
 *   
 * @endcode
 * 
 * Free the memory now that we are done
 * 
 * @code
 *         braid_Destroy(core);
 *   
 *         delete app;
 *   
 * @endcode
 * 
 * Clean up MPI
 * MPI_Comm_free(&comm);
 * 
 * @code
 *         MPI_Finalize();
 *       }
 *     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="ann-test/test_braid.cc"></a>
<h1>Annotated version of test/test_braid.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   #include "BraidFuncs.hh"
 *   
 *   #include <braid.h>
 *   #include <braid_test.h>
 *   
 *   #include <iostream>
 *   
 *   int main(int argc, char** argv)
 *   {
 *     MPI_Comm comm;
 *     int rank;
 *     MPI_Init(&argc, &argv);
 *     comm   = MPI_COMM_WORLD;
 *     MPI_Comm_rank(comm, &rank);
 *   
 *     my_App *app = new(my_App);
 *     app->eq.define();
 *   
 *     double time = 0.2;
 *   
 *     braid_Int init_access_result = braid_TestInitAccess(app,
 *                                                         comm,
 *                                                         stdout,
 *                                                         time,
 *                                                         my_Init,
 *                                                         my_Access,
 *                                                         my_Free);
 *     (void)init_access_result;
 *   
 *     braid_Int clone_result = braid_TestClone(app,
 *                                              comm,
 *                                              stdout,
 *                                              time,
 *                                              my_Init,
 *                                              my_Access,
 *                                              my_Free,
 *                                              my_Clone);
 *     (void)clone_result;
 *   
 *     braid_Int sum_result = braid_TestSum(app,
 *                                          comm,
 *                                          stdout,
 *                                          time,
 *                                          my_Init,
 *                                          my_Access,
 *                                          my_Free,
 *                                          my_Clone,
 *                                          my_Sum);
 *     (void)sum_result;
 *   
 *     braid_Int norm_result = braid_TestSpatialNorm(app,
 *                                                   comm,
 *                                                   stdout,
 *                                                   time,
 *                                                   my_Init,
 *                                                   my_Free,
 *                                                   my_Clone,
 *                                                   my_Sum,
 *                                                   my_SpatialNorm);
 *     (void)norm_result;
 *   
 *     braid_Int buf_result = braid_TestBuf(app,
 *                                          comm,
 *                                          stdout,
 *                                          time,
 *                                          my_Init,
 *                                          my_Free,
 *                                          my_Sum,
 *                                          my_SpatialNorm,
 *                                          my_BufSize,
 *                                          my_BufPack,
 *                                          my_BufUnpack);
 *     (void)buf_result;
 *   
 * @endcode
 * 
 * /* Create spatial communicator for wrapper-tests */
 * braid_SplitCommworld(&comm, 1, &comm_x, &comm_t);
 * 
 * 
 * braid_TestAll(app, comm_x, stdout, 0.0, (tstop-tstart)/ntime,
 * 2*(tstop-tstart)/ntime, my_Init, my_Free, my_Clone,
 * my_Sum, my_SpatialNorm, my_BufSize, my_BufPack,
 * my_BufUnpack, my_Coarsen, my_Interp, my_Residual, my_Step);
 * 
 
 * 
 * 
 * @code
 *     /* Finalize MPI */
 *     MPI_Finalize();
 *   }
 * @endcode
 
 
*/