]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Provide a function for generating nice output of computing times. Use that in step-32.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 12 Jan 2009 16:38:50 +0000 (16:38 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 12 Jan 2009 16:38:50 +0000 (16:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@18185 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/timer.h
deal.II/base/source/timer.cc
deal.II/examples/step-32/step-32.cc

index 28cf778628658330aebffb064fca737afbe51735..a42f49fdb0abc48ad32a19aec69771a082ca1713 100644 (file)
@@ -14,6 +14,9 @@
 #define __deal2__timer_h
 
 #include <base/config.h>
+#include <base/conditional_ostream.h>
+#include <string>
+#include <vector>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -167,6 +170,80 @@ class Timer
     bool                running;
 };
 
+
+
+/**
+ * This class can be used to generate formatted output from time
+ * measurements of different subsections in a program. It is possible to
+ * create several sections that perform certain aspects of the program. A
+ * section can be entered several times. By changing the options in
+ * OutputFrequency and OutputType, the user can choose whether output should
+ * be generated every time a section is joined or just in the end of the
+ * program. Moreover, it is possible to show CPU times, wall times or both.
+ *
+ * TODO: Write a more extensive documentation.
+ *
+ * @ingroup utilities
+ * @author M. Kronbichler, 2009.
+ */
+class TimerOutput
+{
+  public:
+                                  // Sets whether to generate output every
+                                  // time we exit a section, just in the
+                                  // end, or both.
+    enum OutputFrequency {every_call, summary, every_call_and_summary}
+                                                             output_frequency;
+
+                                  // Sets whether to show CPU times, wall
+                                  // times, or both CPU and wall times.
+    enum OutputType      {cpu_times, wall_times, cpu_and_wall_times} 
+                                                             output_type;
+
+                                  // Constructor that takes std::cout as
+                                  // output stream.
+    TimerOutput (std::ostream              &stream, 
+                const enum OutputFrequency output_frequency,
+                const enum OutputType      output_type);
+
+                                  // Constructor that takes a
+                                  // ConditionalOStream to write output to.
+    TimerOutput (ConditionalOStream        &stream, 
+                const enum OutputFrequency output_frequency,
+                const enum OutputType      output_type);
+
+                                  // Destructor. Calls print_summary() in
+                                  // case the option for writing the
+                                  // summary output is set.
+    ~TimerOutput();
+
+                                  // Open a section by given a string name
+                                  // of it. In case the name already
+                                  // exists, that section is done once
+                                  // again.
+    void enter_section (const std::string &section_name);
+
+                                  // Leave a section. If no name is given,
+                                  // the last section that was entered is
+                                  // left.
+    void exit_section (const std::string &section_name = std::string());
+
+                                  // Print a formatted table that
+                                  // summarizes the time consumed in the
+                                  // various sections.
+    void print_summary ();
+
+  private:
+    Timer              timer_all;
+    std::vector<Timer> section_timers;
+    std::vector<std::string> section_names;
+    std::vector<double> section_total_cpu_times;
+    std::vector<double> section_total_wall_times;
+    std::vector<unsigned int> section_n_calls;
+    ConditionalOStream out_stream;
+    std::vector<unsigned int> active_sections;
+};
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index eda66492835dd622890de805b80e84d674f52e7d..f20c70dfffc32315bb6d0d87df5af4ea1b29cdd7 100644 (file)
 
 
 #include <base/timer.h>
-
+#include <base/exceptions.h>
+#include <iostream>
 
 // these includes should probably be properly
 // ./configure'd using the AC_HEADER_TIME macro:
 #include <sys/time.h>
 #include <sys/resource.h>
 
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+#include <mpi.h>
+#endif
+
 
 // on SunOS 4.x, getrusage is stated in the man pages and exists, but
 // is not declared in resource.h. declare it ourselves
@@ -131,4 +136,347 @@ void Timer::reset ()
   running         = false;
 }
 
+
+
+/* ---------------------------- TimerOutput -------------------------- */
+
+TimerOutput::TimerOutput (std::ostream &stream, 
+                         const enum OutputFrequency output_frequency,
+                         const enum OutputType output_type)
+                          :
+                          output_frequency (output_frequency),
+                         output_type (output_type),
+                          out_stream (stream, true)
+{}
+
+
+
+TimerOutput::TimerOutput (ConditionalOStream &stream,
+                         const enum OutputFrequency output_frequency,
+                         const enum OutputType output_type)
+                          :
+                          output_frequency (output_frequency),
+                         output_type (output_type),
+                          out_stream (stream)
+{}
+
+
+
+TimerOutput::~TimerOutput()
+{
+  while (active_sections.size() > 0)
+    exit_section();
+
+  if (output_frequency != every_call)
+    print_summary();
+}
+
+
+
+void 
+TimerOutput::enter_section (const std::string &section_name)
+{
+  Assert (section_name.empty() == false,
+         ExcMessage ("Section string is empty."));
+
+  unsigned int this_section_number = numbers::invalid_unsigned_int;
+
+                                  // check whether the requested section
+                                  // already exists
+  for (unsigned int i=0; i<section_names.size(); ++i)
+    {
+                                  // if we found it, we have to have the
+                                  // index into which we write the
+                                  // data. ensure that the section was not
+                                  // active before to avoid any mess
+      if (section_names[i].compare (section_name) == 0)
+       {
+         this_section_number = i;
+         for (unsigned int j=0; j<active_sections.size(); ++j)
+           {
+             if (active_sections[j] == i)
+               {
+#ifdef DEBUG
+                 std::string exc_text = "Cannot enter the already active section " 
+                   + section_name;
+                 Assert (active_sections[j] != i, ExcMessage (exc_text.data()));
+#endif
+                 exit_section(section_name);
+               }
+           }
+         break;
+       }
+    }
+
+                                  // in case this is a new section, enlarge
+                                  // our data storage by one, otherwise
+                                  // just reset the timer
+  if (this_section_number == numbers::invalid_unsigned_int)
+    {
+      section_names.push_back(section_name);
+      section_total_cpu_times.push_back (0.);
+      section_total_wall_times.push_back (0.);
+      section_timers.push_back (Timer());
+      section_n_calls.push_back (0);
+      this_section_number = section_names.size() - 1;
+    }
+  else
+    {
+      section_timers[this_section_number].reset();
+      section_timers[this_section_number].start();
+    }
+
+                                  // now turn the current section into the
+                                  // list of active sections and increase
+                                  // the counter for that section
+  active_sections.push_back (this_section_number);
+  section_n_calls[this_section_number]++;
+}
+
+
+
+void 
+TimerOutput::exit_section (const std::string &section_name)
+{
+                               // if no string is given, exit the last
+                               // active section.
+  unsigned int active_index_to_delete = numbers::invalid_unsigned_int;
+  if (section_name.empty() == true)
+    active_index_to_delete = active_sections.size() - 1;
+                                  // if we got a string, we need to find
+                                  // the index in the list of active
+                                  // functions
+  else
+    {
+      for (unsigned int i=0; i<active_sections.size(); ++i)
+       if (section_names[active_sections[i]].compare (section_name) == 0)
+         {
+           active_index_to_delete = i;
+           break;
+         }
+    }
+
+  Assert (active_index_to_delete != numbers::invalid_unsigned_int,
+         ExcMessage ("Cannot delete a section that was never created."));
+
+  const unsigned int section_to_exit = active_sections[active_index_to_delete];
+
+  Assert (section_to_exit < section_names.size(),
+         ExcInternalError());
+
+                               // stop the timer for this section.
+  section_timers[section_to_exit].stop();
+  section_total_wall_times[section_to_exit] += section_timers[section_to_exit].wall_time();
+
+                               // get cpu time. on MPI systems, add the
+                               // local contributions.
+                               // 
+                               // TODO: this should rather be in the Timer
+                               // class itself, shouldn't it?
+  double cpu_time = section_timers[section_to_exit]();
+  {
+
+                               // On MPI, sum up all the local CPU
+                               // times.
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    double total_cpu_time = 0.;
+    int mpiInitialized;
+    MPI_Initialized(&mpiInitialized);
+
+    if( mpiInitialized ) 
+      {
+       MPI_Allreduce (&cpu_time, &total_cpu_time, 1, MPI_DOUBLE, MPI_SUM, 
+                      MPI_COMM_WORLD);
+       cpu_time = total_cpu_time;
+      }
+#endif
+    section_total_cpu_times[section_to_exit] += cpu_time;
+  }
+
+
+                                  // in case we have to print out
+                                  // something, do that here...
+  if (output_frequency != summary)
+    {
+      std::string output_time;
+      std::ostringstream cpu;
+      cpu << cpu_time << "s";
+      std::ostringstream wall;
+      wall << section_timers[section_to_exit].wall_time() << "s";
+      if (output_type == cpu_times)
+       output_time = " CPU time: " + cpu.str();
+      else if (output_type == wall_times)
+       output_time = " wall time: " + wall.str() + ".";
+      else
+       output_time = ", CPU/wall time: " + cpu.str() + " / " + wall.str() + ".";
+
+      out_stream << section_names[section_to_exit] << output_time
+                << std::endl;
+    }
+
+                                  // delete the index from the list of
+                                  // active ones
+  std::vector<unsigned int>::iterator position_to_delete 
+    = active_sections.begin() + active_index_to_delete;
+  active_sections.erase(position_to_delete);
+}
+
+
+
+void 
+TimerOutput::print_summary ()
+{
+                               // in case we want to write CPU times
+  if (output_type != wall_times)
+    {
+      double total_cpu_time;
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+      {
+       double my_cpu_time = timer_all();
+       int mpiInitialized;
+       MPI_Initialized(&mpiInitialized);
+
+       if( mpiInitialized ) 
+         {
+           MPI_Allreduce (&my_cpu_time, &total_cpu_time, 1, MPI_DOUBLE, MPI_SUM, 
+                          MPI_COMM_WORLD);
+         }
+       else
+         total_cpu_time = my_cpu_time;
+      }
+#else
+      total_cpu_time = timer_all();
+#endif
+
+                                  // check that the sum of all times is
+                                  // less or equal than the total
+                                  // time. otherwise, we might have
+                                  // generated a lot of overhead in this
+                                  // function.
+      double check_time = 0.;
+      for (unsigned int i=0; i<section_names.size(); ++i)
+       {
+         check_time += section_total_cpu_times[i];
+       }
+      if (check_time > total_cpu_time)
+       {
+         total_cpu_time = check_time;
+         out_stream << std::endl << "Sum of counted times is larger than total time. "
+                    << "Timer function may have introduced too much overhead." << std::endl;
+       }
+
+                                  // now generate a nice table
+      out_stream << "\n\n"
+                << "+---------------------------------------------+------------"
+                << "+------------+\n"
+                << "| Total CPU time elapsed since start          |";
+      std::cout.width(10);
+      std::cout.precision(3);
+      out_stream << total_cpu_time << "s |            |\n";
+      out_stream << "|                                             |            "
+                << "|            |\n";
+      out_stream << "| Section                         | no. calls |";
+      std::cout.width(10);
+      std::cout.precision(3);
+      out_stream << "  CPU time "  << " | % of total |\n";
+      out_stream << "+---------------------------------------------+------------"
+                << "+------------+";
+      for (unsigned int i=0; i<section_names.size(); ++i)
+       {
+         std::string name_out = section_names[i];
+
+                               // resize the array so that it is always
+                               // of the same size
+         unsigned int pos_non_space = name_out.find_first_not_of (" ");
+         name_out.erase(0, pos_non_space);
+         name_out.resize (32, ' ');
+         out_stream << std::endl;
+         out_stream << "| " << name_out;
+         out_stream << "| ";
+         std::cout.width(9);
+         out_stream << section_n_calls[i] << " |";
+         std::cout.width(10);
+         std::cout.precision(3);
+         out_stream << section_total_cpu_times[i] << "s |";
+         std::cout.width(10);
+         std::cout.precision(2);
+         out_stream << section_total_cpu_times[i]/total_cpu_time * 100 << "% |";
+       }
+      out_stream << std::endl
+                << "+---------------------------------------------+"
+                << "------------+------------+\n"
+                << std::endl;
+    }
+
+                               // in case we want to write out wallclock times
+  if (output_type != cpu_times)
+    {
+      double total_wall_time = timer_all.wall_time();
+
+                                  // check that the sum of all times is
+                                  // less or equal than the total
+                                  // time. otherwise, we might have
+                                  // generated a lot of overhead in this
+                                  // function.
+      double check_time = 0.;
+      for (unsigned int i=0; i<section_names.size(); ++i)
+       {
+         check_time += section_total_wall_times[i];
+       }
+      if (check_time > total_wall_time)
+       {
+         total_wall_time = check_time;
+         out_stream << std::endl 
+                    << "Sum of counted times is larger than total time. "
+                    << "Timer function may have introduced too much overhead." 
+                    << std::endl;
+       }
+
+                                  // now generate a nice table
+     out_stream << "\n\n"
+                << "+---------------------------------------------+------------"
+                << "+------------+\n"
+                << "| Total wallclock time elapsed from start     |";
+      std::cout.width(10);
+      std::cout.precision(3);
+      out_stream << total_wall_time << "s |            |\n";
+      out_stream << "|                                             |            "
+                << "|            |\n";
+      out_stream << "| Section                         | no. calls |";
+      std::cout.width(10);
+      std::cout.precision(3);
+      out_stream << "  CPU time "  << " | % of total |\n";
+      out_stream << "+---------------------------------------------+------------"
+                << "+------------+";
+      for (unsigned int i=0; i<section_names.size(); ++i)
+       {
+         std::string name_out = section_names[i];
+
+                               // resize the array so that it is always
+                               // of the same size
+         unsigned int pos_non_space = name_out.find_first_not_of (" ");
+         name_out.erase(0, pos_non_space);
+         name_out.resize (32, ' ');
+         out_stream << std::endl;
+         out_stream << "| " << name_out;
+         out_stream << "| ";
+         std::cout.width(9);
+         out_stream << section_n_calls[i] << " |";
+         std::cout.width(10);
+         std::cout.precision(3);
+         out_stream << section_total_wall_times[i] << "s |";
+         std::cout.width(10);
+         std::cout.precision(2);
+         out_stream << section_total_wall_times[i]/total_wall_time * 100 << "% |";
+       }
+      out_stream << std::endl
+                << "+---------------------------------------------+"
+                << "------------+------------+\n"
+                << std::endl;
+    }
+
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
index 76191ab78cdfb556c785c367235def57d48a5025..77684c6587102485ac7556eb85a6cafeec5ba315 100644 (file)
@@ -56,6 +56,9 @@
 
 #include <Epetra_Map.h>
 
+                                       // Time measurements. 
+#include <base/timer.h>
+
 #include <fstream>
 #include <iostream>
 #include <sstream>
@@ -65,7 +68,6 @@
                                 // names into global namespace
 using namespace dealii;
 
-
                                 // @sect3{Equation data}
 
 
@@ -95,10 +97,11 @@ namespace EquationData
 
   template <int dim>
   double
-  TemperatureInitialValues<dim>::value (const Point<dim>  &p,
+  TemperatureInitialValues<dim>::value (const Point<dim>  &,
                                        const unsigned int) const
   {
-    return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);
+    //return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);
+    return 0.;
   }
 
 
@@ -130,10 +133,32 @@ namespace EquationData
 
   template <int dim>
   double
-  TemperatureRightHandSide<dim>::value (const Point<dim>  &,
-                                       const unsigned int /*component*/) const
+  TemperatureRightHandSide<dim>::value (const Point<dim>  &p,
+                                       const unsigned int component) const
   {
-    return 0;
+    //    return 0;
+    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);
+
   }
 
 
@@ -278,7 +303,6 @@ class BoussinesqFlowProblem
     void output_results () const;
     void refine_mesh (const unsigned int max_grid_level);
 
-    static
     double
     compute_viscosity(const std::vector<double>          &old_temperature,
                      const std::vector<double>          &old_old_temperature,
@@ -286,13 +310,12 @@ class BoussinesqFlowProblem
                      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> >  &present_velocity_values,
+                     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                        global_Omega_diameter,
-                     const double                        cell_diameter,
-                     const double                        old_time_step);
+                     const double                        cell_diameter);
 
 
     const Epetra_Comm                  &trilinos_communicator;
@@ -300,6 +323,7 @@ class BoussinesqFlowProblem
     ConditionalOStream                  pcout;
 
     Triangulation<dim>                  triangulation;
+    double                              global_Omega_diameter;
 
     const unsigned int                  stokes_degree;
     FESystem<dim>                       stokes_fe;
@@ -311,6 +335,7 @@ class BoussinesqFlowProblem
     TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
 
     TrilinosWrappers::MPI::BlockVector  stokes_solution;
+    TrilinosWrappers::BlockVector       old_stokes_solution;
     TrilinosWrappers::MPI::BlockVector  stokes_rhs;
 
 
@@ -334,12 +359,16 @@ class BoussinesqFlowProblem
     double old_time_step;
     unsigned int timestep_number;
 
-    boost::shared_ptr<TrilinosWrappers::PreconditionAMG>  Amg_preconditioner;
-    boost::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner;
+    boost::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
+    boost::shared_ptr<TrilinosWrappers::PreconditionIC>  Mp_preconditioner;
+    boost::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;
 
     bool rebuild_stokes_matrix;
-    bool rebuild_temperature_matrices;
     bool rebuild_stokes_preconditioner;
+    bool rebuild_temperature_matrices;
+    bool rebuild_temperature_preconditioner;
+
+    TimerOutput computing_timer;
 };
 
 
@@ -369,8 +398,11 @@ BoussinesqFlowProblem<dim>::BoussinesqFlowProblem ()
                old_time_step (0),
                timestep_number (0),
                rebuild_stokes_matrix (true),
+               rebuild_stokes_preconditioner (true),
                rebuild_temperature_matrices (true),
-               rebuild_stokes_preconditioner (true)
+               rebuild_temperature_preconditioner (true),
+               computing_timer (pcout, TimerOutput::summary,
+                                TimerOutput::wall_times)
 {}
 
 
@@ -513,35 +545,34 @@ BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
 template <int dim>
 double
 BoussinesqFlowProblem<dim>::
-compute_viscosity(const std::vector<double>          &old_temperature,
-                 const std::vector<double>          &old_old_temperature,
-                 const std::vector<Tensor<1,dim> >  &old_temperature_grads,
-                 const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
-                 const std::vector<double>          &old_temperature_laplacians,
-                 const std::vector<double>          &old_old_temperature_laplacians,
-                 const std::vector<Tensor<1,dim> >  &present_velocity_values,
-                 const std::vector<double>          &gamma_values,
-                 const double                        global_u_infty,
-                 const double                        global_T_variation,
-                 const double                        global_Omega_diameter,
-                 const double                        cell_diameter,
-                 const double                        old_time_step)
+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 double beta = 0.015 * dim;
-  const double alpha = 1;
+  const double beta = 0.04 * dim;
+  const double alpha = 2;
   
   if (global_u_infty == 0)
     return 5e-3 * cell_diameter;
   
   const unsigned int n_q_points = old_temperature.size();
   
-                                  // Stage 1: calculate residual
   double max_residual = 0;
   double max_velocity = 0;
   
   for (unsigned int q=0; q < n_q_points; ++q)
     {
-      const Tensor<1,dim> u = present_velocity_values[q];
+      const Tensor<1,dim> u = (old_velocity_values[q] +
+                              old_old_velocity_values[q]) / 2;
       
       const double dT_dt = (old_temperature[q] - old_old_temperature[q])
                           / old_time_step;
@@ -567,15 +598,18 @@ compute_viscosity(const std::vector<double>          &old_temperature,
   return (beta *
          max_velocity *
          std::min (cell_diameter,
-                   std::pow(cell_diameter,alpha) * max_residual / global_scaling));
+                   std::pow(cell_diameter,alpha) *
+                   max_residual / global_scaling));
 }
 
 
 
+
                                 // @sect4{BoussinesqFlowProblem::setup_dofs}
 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs ()
 {
+  computing_timer.enter_section("Setup dof systems");
   std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
   stokes_sub_blocks[dim] = 1;
 
@@ -664,6 +698,13 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
     sp.compress();
 
     stokes_matrix.reinit (sp);
+    /*std::cout << "Processor " << trilinos_communicator.MyPID() 
+             << " stokes(0,0) rows: " 
+             << stokes_matrix.block(0,0).matrix->NumMyRows()
+             << ", nnz: " 
+             << stokes_matrix.block(0,0).matrix->NumMyNonzeros() 
+             << std::endl;*/
+
   }
 
   {
@@ -696,6 +737,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                  0,
                  trilinos_communicator);
   {
+    T_preconditioner.reset ();
     temperature_mass_matrix.clear ();
     temperature_stiffness_matrix.clear ();
     temperature_matrix.clear ();
@@ -711,12 +753,15 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
   }
 
   stokes_solution.reinit (stokes_partitioner);
+  old_stokes_solution.reinit (stokes_partitioner);
   stokes_rhs.reinit (stokes_partitioner);
 
   temperature_solution.reinit (temperature_partitioner);
   old_temperature_solution.reinit (temperature_partitioner);
   old_old_temperature_solution.reinit (temperature_partitioner);
   temperature_rhs.reinit (temperature_partitioner);
+
+  computing_timer.exit_section();
 }
 
 
@@ -791,6 +836,7 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   if (rebuild_stokes_preconditioner == false)
     return;
 
+  computing_timer.enter_section ("   Build Stokes preconditioner");
   pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
 
   assemble_stokes_preconditioner ();
@@ -806,16 +852,17 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
   
   Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
                                 TrilinosWrappers::PreconditionAMG::AdditionalData
-                                (true, true, 5e-2, null_space, 3, 0, false));
+                                (true, true, 1e-2, null_space, 3, 0, false));
 
   Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionIC>
                                    (new TrilinosWrappers::PreconditionIC());
 
   Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
 
-  pcout << std::endl;
-
   rebuild_stokes_preconditioner = false;
+
+  pcout << std::endl;
+  computing_timer.exit_section();
 }
 
 
@@ -826,6 +873,8 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 {
   pcout << "   Assembling..." << std::flush;
 
+  computing_timer.enter_section ("   Assemble Stokes system");
+
   if (rebuild_stokes_matrix == true)
     stokes_matrix=0;
 
@@ -860,7 +909,7 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
   std::vector<SymmetricTensor<2,dim> > grads_phi_u (dofs_per_cell);
   std::vector<double>                  div_phi_u   (dofs_per_cell);
   std::vector<double>                  phi_p       (dofs_per_cell);
-  
+
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
@@ -934,6 +983,7 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
   rebuild_stokes_matrix = false;
 
   pcout << std::endl;
+  computing_timer.exit_section();
 }
 
 
@@ -947,7 +997,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
 {
   if (rebuild_temperature_matrices == false)
     return;
-  
+
+  computing_timer.enter_section ("   Assemble temperature matrices");
   temperature_mass_matrix = 0;
   temperature_stiffness_matrix = 0;
   
@@ -1015,6 +1066,9 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
   temperature_stiffness_matrix.compress();
 
   rebuild_temperature_matrices = false;
+  rebuild_temperature_preconditioner = true;
+
+  computing_timer.exit_section();
 }
 
 
@@ -1038,6 +1092,15 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
       temperature_matrix.add (time_step, temperature_stiffness_matrix);
     }
   temperature_matrix.compress();
+
+  if (rebuild_temperature_preconditioner == true)
+    {
+      T_preconditioner =  boost::shared_ptr<TrilinosWrappers::PreconditionIC>
+                                   (new TrilinosWrappers::PreconditionIC());
+      T_preconditioner->initialize (temperature_matrix);
+
+      rebuild_temperature_preconditioner = false;
+    }
   
   temperature_rhs = 0;
   
@@ -1059,8 +1122,8 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-  std::vector<Tensor<1,dim> > present_velocity_values (n_q_points);
-
+  std::vector<Tensor<1,dim> > old_velocity_values (n_q_points);
+  std::vector<Tensor<1,dim> > old_old_velocity_values (n_q_points);
   
   std::vector<double>         old_temperature_values (n_q_points);
   std::vector<double>         old_old_temperature_values(n_q_points);
@@ -1081,7 +1144,10 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
     global_T_range = get_extrapolated_temperature_range();
   const double global_Omega_diameter = GridTools::diameter (triangulation);
 
-  const TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution);
+  const TrilinosWrappers::BlockVector 
+    localized_stokes_solution (stokes_solution);
+  const TrilinosWrappers::BlockVector 
+    localized_old_stokes_solution (old_stokes_solution);
 
   const FEValuesExtractors::Vector velocities (0);
   
@@ -1117,9 +1183,11 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
        
        temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
                                                gamma_values);
-  
+
        stokes_fe_values[velocities].get_function_values (localized_stokes_solution,
-                                                         present_velocity_values);
+                                                         old_velocity_values);
+       stokes_fe_values[velocities].get_function_values (localized_old_stokes_solution,
+                                                         old_old_velocity_values);
        
        const double nu
          = compute_viscosity (old_temperature_values,
@@ -1128,12 +1196,12 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                               old_old_temperature_grads,
                               old_temperature_laplacians,
                               old_old_temperature_laplacians,
-                              present_velocity_values,
+                              old_velocity_values,
+                              old_old_velocity_values,
                               gamma_values,
                               global_u_infty,
                               global_T_range.second - global_T_range.first,
-                              global_Omega_diameter, cell->diameter(),
-                              old_time_step);
+                              cell->diameter());
        
        for (unsigned int q=0; q<n_q_points; ++q)
          {
@@ -1143,61 +1211,48 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                phi_T[k]      = temperature_fe_values.shape_value (k, q);
              }
   
-           const double        old_T      = old_temperature_values[q];
-           const double        old_old_T  = old_old_temperature_values[q];
-  
-           const Tensor<1,dim> old_grad_T     = old_temperature_grads[q];
-           const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q];
-  
-           
-           const Tensor<1,dim> present_u = present_velocity_values[q];
-  
-           if (use_bdf2_scheme == true)
-             {
-               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                 local_rhs(i) += ((time_step + old_time_step) / old_time_step *
-                                  old_T * phi_T[i]
-                                  -
-                                  (time_step * time_step) /
-                                  (old_time_step * (time_step + old_time_step)) *
-                                  old_old_T * phi_T[i]
-                                  -
-                                  time_step *
-                                  present_u *
-                                  ((1+time_step/old_time_step) * old_grad_T
-                                   -
-                                   time_step / old_time_step * old_old_grad_T) *
-                                  phi_T[i]
-                                  -
-                                  time_step *
-                                  nu *
-                                  ((1+time_step/old_time_step) * old_grad_T
-                                   -
-                                   time_step / old_time_step * old_old_grad_T) *
-                                  grad_phi_T[i]
-                                  +
-                                  time_step *
-                                  gamma_values[q] * phi_T[i])
-                                 * 
-                                 temperature_fe_values.JxW(q);
-             }
-           else
-             {
-               for (unsigned int i=0; i<dofs_per_cell; ++i)
-                 local_rhs(i) += (old_T * phi_T[i]
-                                  -
-                                  time_step *
-                                  present_u * old_grad_T * phi_T[i]
-                                  -
-                                  time_step *
-                                  nu *
-                                  old_grad_T * grad_phi_T[i]
-                                  +
-                                  time_step *
-                                  gamma_values[q] * phi_T[i])
-                                 *
-                                 temperature_fe_values.JxW(q);
-             }
+
+         const double old_Ts
+           = (use_bdf2_scheme ?
+              (old_temperature_values[q] *
+               (time_step + old_time_step) / old_time_step
+               -
+               old_old_temperature_values[q] *
+               (time_step * time_step) /
+               (old_time_step * (time_step + old_time_step)))
+              :
+              old_temperature_values[q]);
+
+         const Tensor<1,dim> ext_grad_T
+           = (use_bdf2_scheme ? 
+              (old_temperature_grads[q] *
+               (1+time_step/old_time_step) 
+               -
+               old_old_temperature_grads[q] *
+               time_step / old_time_step) 
+              :
+              old_temperature_grads[q]);
+         
+         const Tensor<1,dim> extrapolated_u
+           = (use_bdf2_scheme ? 
+              (old_velocity_values[q] * (1+time_step/old_time_step) - 
+               old_old_velocity_values[q] * time_step/old_time_step)
+              :
+              old_velocity_values[q]);
+
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           local_rhs(i) += (old_Ts * phi_T[i]
+                            -
+                            time_step *
+                            extrapolated_u * ext_grad_T * phi_T[i]
+                            -
+                            time_step *
+                            nu * ext_grad_T * grad_phi_T[i]
+                            +
+                            time_step *
+                            gamma_values[q] * phi_T[i])
+                           *
+                           temperature_fe_values.JxW(q);
          }
        
        cell->get_dof_indices (local_dof_indices);
@@ -1216,6 +1271,7 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
 template <int dim>
 void BoussinesqFlowProblem<dim>::solve ()
 {
+  computing_timer.enter_section ("   Solve Stokes system");
   pcout << "   Solving..." << std::endl;
 
   {
@@ -1245,31 +1301,36 @@ void BoussinesqFlowProblem<dim>::solve ()
     stokes_constraints.distribute (localized_stokes_solution);
     stokes_solution = localized_stokes_solution;
   }
+  computing_timer.exit_section();
+
+
+  computing_timer.enter_section ("   Assemble temperature rhs");
 
   old_time_step = time_step;    
-  time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
+  time_step = 1./(1.9*dim*std::sqrt(1.*dim)) /
              temperature_degree *
              GridTools::minimal_cell_diameter(triangulation) /
-              std::max (get_maximal_velocity(), 1.e-5);
+              std::max (get_maximal_velocity(), 0.01);
   
   pcout << "   " << "Time step: " << time_step
-       << std::endl;
-  
+       << std::endl;  
+
   temperature_solution = old_temperature_solution;
 
 
   assemble_temperature_system ();
-  {
 
+  computing_timer.exit_section ();
+
+  computing_timer.enter_section ("   Solve temperature system");
+
+  {
     SolverControl solver_control (temperature_matrix.m(),
                                  1e-8*temperature_rhs.l2_norm());
     SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
 
-    TrilinosWrappers::PreconditionIC preconditioner;
-    preconditioner.initialize (temperature_matrix);
-
     cg.solve (temperature_matrix, temperature_solution,
-             temperature_rhs, preconditioner);
+             temperature_rhs, *T_preconditioner);
 
     TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
     temperature_constraints.distribute (localized_temperature_solution);
@@ -1277,12 +1338,12 @@ void BoussinesqFlowProblem<dim>::solve ()
 
     pcout << "   "
          << solver_control.last_step()
-         << " CG iterations for temperature."
-         << std::endl;
+         << " CG iterations for temperature" << std::endl;
+    computing_timer.exit_section();
 
     double min_temperature = localized_temperature_solution(0),
           max_temperature = localized_temperature_solution(0);
-    for (unsigned int i=0; i<temperature_solution.size(); ++i)
+    for (unsigned int i=1; i<temperature_solution.size(); ++i)
       {
        min_temperature = std::min<double> (min_temperature,
                                            localized_temperature_solution(i));
@@ -1312,60 +1373,60 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
   Assert (joint_dof_handler.n_dofs() ==
          stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(),
          ExcInternalError());
-  
+
   Vector<double> joint_solution (joint_dof_handler.n_dofs());
   TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution);
   TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
 
-  {
-    std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
-    std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
-    std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
-    
-    typename DoFHandler<dim>::active_cell_iterator
-      joint_cell       = joint_dof_handler.begin_active(),
-      joint_endc       = joint_dof_handler.end(),
-      stokes_cell      = stokes_dof_handler.begin_active(),
-      temperature_cell = temperature_dof_handler.begin_active();
-    for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
+  if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0)
+    {
+
       {
-       joint_cell->get_dof_indices (local_joint_dof_indices);
-       stokes_cell->get_dof_indices (local_stokes_dof_indices);
-       temperature_cell->get_dof_indices (local_temperature_dof_indices);
-
-       for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
-         if (joint_fe.system_to_base_index(i).first.first == 0)
-           {
-             Assert (joint_fe.system_to_base_index(i).second
-                     <
-                     local_stokes_dof_indices.size(),
-                     ExcInternalError());
-             joint_solution(local_joint_dof_indices[i])
-               = localized_stokes_solution(local_stokes_dof_indices[joint_fe.system_to_base_index(i).second]);
-           }
-         else
-           {
-             Assert (joint_fe.system_to_base_index(i).first.first == 1,
-                     ExcInternalError());
-             Assert (joint_fe.system_to_base_index(i).second
-                     <
-                     local_stokes_dof_indices.size(),
-                     ExcInternalError());
-             joint_solution(local_joint_dof_indices[i])
-               = localized_temperature_solution(local_temperature_dof_indices[joint_fe.system_to_base_index(i).second]);
-           }
+       std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
+       std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
+       std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
+    
+       typename DoFHandler<dim>::active_cell_iterator
+         joint_cell       = joint_dof_handler.begin_active(),
+         joint_endc       = joint_dof_handler.end(),
+         stokes_cell      = stokes_dof_handler.begin_active(),
+         temperature_cell = temperature_dof_handler.begin_active();
+       for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
+         {
+           joint_cell->get_dof_indices (local_joint_dof_indices);
+           stokes_cell->get_dof_indices (local_stokes_dof_indices);
+           temperature_cell->get_dof_indices (local_temperature_dof_indices);
+
+           for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
+             if (joint_fe.system_to_base_index(i).first.first == 0)
+               {
+                 Assert (joint_fe.system_to_base_index(i).second
+                         <
+                         local_stokes_dof_indices.size(),
+                         ExcInternalError());
+                 joint_solution(local_joint_dof_indices[i])
+                   = localized_stokes_solution(local_stokes_dof_indices[joint_fe.system_to_base_index(i).second]);
+               }
+             else
+               {
+                 Assert (joint_fe.system_to_base_index(i).first.first == 1,
+                         ExcInternalError());
+                 Assert (joint_fe.system_to_base_index(i).second
+                         <
+                         local_stokes_dof_indices.size(),
+                         ExcInternalError());
+                 joint_solution(local_joint_dof_indices[i])
+                   = localized_temperature_solution(local_temperature_dof_indices[joint_fe.system_to_base_index(i).second]);
+               }
+         }
       }
-  }
-  
-  
-  std::vector<std::string> joint_solution_names (dim, "velocity");
-  joint_solution_names.push_back ("p");
-  joint_solution_names.push_back ("T");
 
-  DataOut<dim> data_out;
+      std::vector<std::string> joint_solution_names (dim, "velocity");
+      joint_solution_names.push_back ("p");
+      joint_solution_names.push_back ("T");
+
+      DataOut<dim> data_out;
 
-  if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0)
-    {
       data_out.attach_dof_handler (joint_dof_handler);
 
       std::vector<DataComponentInterpretation::DataComponentInterpretation>
@@ -1394,6 +1455,7 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
 template <int dim>
 void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
 {
+  computing_timer.enter_section ("Refine mesh structure, part 1");
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
   TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
@@ -1413,29 +1475,45 @@ void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
         cell != triangulation.end(); ++cell)
       cell->clear_refine_flag ();
   
-  std::vector<TrilinosWrappers::Vector > x_solution (2);
-  x_solution[0] = temperature_solution;
-  x_solution[1] = old_temperature_solution;
+  std::vector<TrilinosWrappers::Vector> x_temperature (2);
+  x_temperature[0] = temperature_solution;
+  x_temperature[1] = old_temperature_solution;
+  TrilinosWrappers::BlockVector x_stokes = stokes_solution;
 
-  SolutionTransfer<dim,TrilinosWrappers::Vector > soltrans(temperature_dof_handler);
+  SolutionTransfer<dim,TrilinosWrappers::Vector>
+    temperature_trans(temperature_dof_handler);
+  SolutionTransfer<dim,TrilinosWrappers::BlockVector>
+    stokes_trans(stokes_dof_handler);
 
   triangulation.prepare_coarsening_and_refinement();
-  soltrans.prepare_for_coarsening_and_refinement(x_solution);
+  temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
+  stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
 
   triangulation.execute_coarsening_and_refinement ();
+  computing_timer.exit_section();
+
   setup_dofs ();
 
-  std::vector<TrilinosWrappers::Vector > tmp (2);
-  tmp[0] = temperature_solution;
-  tmp[1] = temperature_solution;
-  soltrans.interpolate(x_solution, tmp);
+  computing_timer.enter_section ("Refine mesh structure, part 2");
+
+  std::vector<TrilinosWrappers::Vector> tmp (2);
+  tmp[0].reinit (temperature_solution);
+  tmp[1].reinit (temperature_solution);
+  temperature_trans.interpolate(x_temperature, tmp);
 
   temperature_solution = tmp[0];
   old_temperature_solution = tmp[1];
 
-  rebuild_stokes_matrix         = true;
-  rebuild_temperature_matrices  = true;
-  rebuild_stokes_preconditioner = true;
+  TrilinosWrappers::BlockVector x_stokes_new = stokes_solution;
+  stokes_trans.interpolate (x_stokes, x_stokes_new);
+  stokes_solution = x_stokes_new;
+
+  rebuild_stokes_matrix              = true;
+  rebuild_stokes_preconditioner      = true;
+  rebuild_temperature_matrices       = true;
+  rebuild_temperature_preconditioner = true;
+
+  computing_timer.exit_section();
 }
 
 
@@ -1447,11 +1525,13 @@ void BoussinesqFlowProblem<dim>::run ()
   const unsigned int initial_refinement = (dim == 2 ? 3 : 2);
   const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3);
 
-  GridGenerator::half_hyper_shell (triangulation,
-                                  Point<dim>(), 0.5, 1.0);
+  //GridGenerator::half_hyper_shell (triangulation,
+  //                              Point<dim>(), 0.5, 1.0);
 
-  static HalfHyperShellBoundary<dim> boundary;
-  triangulation.set_boundary (0, boundary);
+  //static HyperShellBoundary<dim> boundary;
+  //triangulation.set_boundary (0, boundary);
+  GridGenerator::hyper_cube (triangulation);
+  global_Omega_diameter = GridTools::diameter (triangulation);
 
   triangulation.refine_global (initial_refinement);
 
@@ -1461,11 +1541,72 @@ void BoussinesqFlowProblem<dim>::run ()
   
   start_time_iteration:
 
-  VectorTools::project (temperature_dof_handler,
-                       temperature_constraints,
-                       QGauss<dim>(temperature_degree+2),
-                       EquationData::TemperatureInitialValues<dim>(),
-                       old_temperature_solution);
+  //VectorTools::project (temperature_dof_handler,
+  //           temperature_constraints,
+  //           QGauss<dim>(temperature_degree+2),
+  //           EquationData::TemperatureInitialValues<dim>(),
+  //           old_temperature_solution);
+  assemble_temperature_matrix ();
+                               // Create right hand side for projection
+  {
+    QGauss<dim> quadrature(temperature_degree+2);
+    UpdateFlags update_flags = UpdateFlags(update_values   |
+                                          update_quadrature_points |
+                                          update_JxW_values);
+    FEValues<dim> fe_values (temperature_fe, quadrature, update_flags);
+
+    const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+      n_q_points    = fe_values.n_quadrature_points;
+  
+    std::vector<unsigned int> dofs (dofs_per_cell);
+    Vector<double> cell_vector (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = temperature_dof_handler.begin_active(),
+      endc = temperature_dof_handler.end();
+
+    std::vector<double> rhs_values(n_q_points);
+
+    TrilinosWrappers::MPI::Vector rhs (temperature_mass_matrix.matrix->RowMap()), 
+      sol (temperature_mass_matrix.matrix->RowMap());
+      
+    for (; cell!=endc; ++cell) 
+      if (cell->subdomain_id() == 
+         Utilities::Trilinos::get_this_mpi_process(trilinos_communicator))
+       {
+         fe_values.reinit(cell);
+         
+         const std::vector<double> &weights   = fe_values.get_JxW_values ();
+         EquationData::TemperatureInitialValues<dim>().value_list 
+           (fe_values.get_quadrature_points(), rhs_values);
+         
+         cell_vector = 0;
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i) 
+             cell_vector(i) += rhs_values[point] *
+                               fe_values.shape_value(i,point) *
+                               weights[point];
+
+         cell->get_dof_indices (dofs);
+        
+         temperature_constraints.distribute_local_to_global (cell_vector,
+                                                             dofs,
+                                                             rhs);
+       }
+    
+    ReductionControl        control(5*rhs.size(), 0., 1e-12, false, false);
+    GrowingVectorMemory<TrilinosWrappers::MPI::Vector> memory;
+    SolverCG<TrilinosWrappers::MPI::Vector> cg(control,memory);
+
+    TrilinosWrappers::PreconditionIC prec;
+    prec.initialize(temperature_mass_matrix);
+                                  // solve
+    cg.solve (temperature_mass_matrix, sol, rhs, prec);
+  
+                                  // distribute solution
+    old_temperature_solution = sol;
+    temperature_constraints.distribute (old_temperature_solution);
+  }
   
   timestep_number           = 0;
   time_step = old_time_step = 0;
@@ -1502,10 +1643,12 @@ void BoussinesqFlowProblem<dim>::run ()
       time += time_step;
       ++timestep_number;
 
+      old_stokes_solution          = stokes_solution;
       old_old_temperature_solution = old_temperature_solution;
       old_temperature_solution     = temperature_solution;      
     }
-  while (time <= 100);
+  while (time <= 0);
+
 }
 
 
@@ -1519,7 +1662,7 @@ int main (int argc, char *argv[])
       
       Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
 
-      BoussinesqFlowProblem<2> flow_problem;
+      BoussinesqFlowProblem<3> flow_problem;
       flow_problem.run ();
     }
   catch (std::exception &exc)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.