]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Replace version with timers by a cleaned-up version.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 28 Oct 2010 15:32:51 +0000 (15:32 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 28 Oct 2010 15:32:51 +0000 (15:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@22537 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-40/step-40-commented.cc [deleted file]
deal.II/examples/step-40/step-40.cc

diff --git a/deal.II/examples/step-40/step-40-commented.cc b/deal.II/examples/step-40/step-40-commented.cc
deleted file mode 100644 (file)
index b949234..0000000
+++ /dev/null
@@ -1,406 +0,0 @@
-/* $Id$ */
-/* Author: Wolfgang Bangerth, Texas A&M University, 2009, 2010 */
-/*         Timo Heister, University of Goettingen, 2009, 2010 */
-
-/*    $Id$       */
-/*                                                                */
-/*    Copyright (C) 2009, 2010 by Timo Heister and the deal.II authors */
-/*                                                                */
-/*    This file is subject to QPL and may not be  distributed     */
-/*    without copyright and license information. Please refer     */
-/*    to the file deal.II/doc/license.html for the  text  and     */
-/*    further information on this license.                        */
-
-#include <base/quadrature_lib.h>
-#include <base/function.h>
-#include <base/utilities.h>
-#include <base/conditional_ostream.h>
-#include <base/index_set.h>
-#include <base/timer.h>
-#include <lac/vector.h>
-#include <lac/full_matrix.h>
-#include <lac/solver_cg.h>
-#include <lac/constraint_matrix.h>
-#include <lac/compressed_simple_sparsity_pattern.h>
-#include <lac/sparsity_tools.h>
-
-#include <lac/petsc_parallel_sparse_matrix.h>
-#include <lac/petsc_parallel_vector.h>
-#include <lac/petsc_solver.h>
-#include <lac/petsc_precondition.h>
-
-#include <grid/grid_generator.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe_values.h>
-#include <fe/fe_q.h>
-#include <numerics/vectors.h>
-#include <numerics/data_out.h>
-#include <numerics/error_estimator.h>
-
-#include <distributed/tria.h>
-#include <distributed/grid_refinement.h>
-
-#include <fstream>
-#include <iostream>
-
-using namespace dealii;
-
-
-template <int dim>
-class LaplaceProblem
-{
-  public:
-    LaplaceProblem ();
-    ~LaplaceProblem ();
-
-    void run (const unsigned int initial_global_refine);
-
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    MPI_Comm mpi_communicator;
-
-    parallel::distributed::Triangulation<dim>   triangulation;
-
-    DoFHandler<dim>      dof_handler;
-    FE_Q<dim>            fe;
-
-    IndexSet             locally_owned_dofs;
-    IndexSet             locally_relevant_dofs;
-
-    ConstraintMatrix     constraints;
-
-    PETScWrappers::MPI::SparseMatrix system_matrix;
-    PETScWrappers::MPI::Vector locally_relevant_solution;
-    PETScWrappers::MPI::Vector system_rhs;
-
-    ConditionalOStream                pcout;
-};
-
-
-
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem ()
-               :
-               mpi_communicator (MPI_COMM_WORLD),
-               triangulation (mpi_communicator,
-                              typename Triangulation<dim>::MeshSmoothing
-                              (Triangulation<dim>::smoothing_on_refinement |
-                               Triangulation<dim>::smoothing_on_coarsening)),
-               dof_handler (triangulation),
-                fe (2),
-               pcout (std::cout,
-                      (Utilities::System::
-                       get_this_mpi_process(mpi_communicator)
-                       == 0))
-{}
-
-
-
-template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
-{
-  dof_handler.clear ();
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::setup_system ()
-{
-  locally_relevant_solution.reinit (mpi_communicator,
-                                   locally_owned_dofs,
-                                   locally_relevant_dofs);
-  system_rhs.reinit (mpi_communicator,
-                    dof_handler.n_dofs(),
-                    dof_handler.n_locally_owned_dofs());
-
-  locally_relevant_solution = 0;
-  system_rhs = 0;
-
-  constraints.clear ();
-  constraints.reinit (locally_relevant_dofs);
-
-  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           constraints);
-  constraints.close ();
-
-
-  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                      dof_handler.n_dofs(),
-                                      locally_relevant_dofs);
-  DoFTools::make_sparsity_pattern (dof_handler,
-                                  csp,
-                                  constraints, false);
-  SparsityTools::distribute_sparsity_pattern (csp,
-                                             dof_handler.n_locally_owned_dofs_per_processor(),
-                                             mpi_communicator,
-                                             locally_relevant_dofs);
-  system_matrix.reinit (mpi_communicator,
-                       csp,
-                       dof_handler.n_locally_owned_dofs_per_processor(),
-                       dof_handler.n_locally_owned_dofs_per_processor(),
-                       Utilities::System::get_this_mpi_process(mpi_communicator));
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
-  const QGauss<dim>  quadrature_formula(3);
-
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                          update_values    |  update_gradients |
-                          update_quadrature_points |
-                          update_JxW_values);
-
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
-
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double>       cell_rhs (dofs_per_cell);
-
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
-      {
-       cell_matrix = 0;
-       cell_rhs = 0;
-
-       fe_values.reinit (cell);
-
-       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           {
-             for (unsigned int j=0; j<dofs_per_cell; ++j)
-               cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
-                                    fe_values.shape_grad(j,q_point) *
-                                    fe_values.JxW(q_point));
-
-             cell_rhs(i) += (fe_values.shape_value(i,q_point) *
-                             (fe_values.quadrature_point(q_point)[1]
-                              >
-                              0.5+0.25*sin(4.0*numbers::PI*fe_values.quadrature_point(q_point)[0])
-                              ? 1 : -1) *
-                             fe_values.JxW(q_point));
-           }
-
-       cell->get_dof_indices (local_dof_indices);
-       constraints.distribute_local_to_global (cell_matrix,
-                                               cell_rhs,
-                                               local_dof_indices,
-                                               system_matrix,
-                                               system_rhs);
-      }
-
-  system_matrix.compress ();
-  system_rhs.compress ();
-}
-
-
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
-  PETScWrappers::MPI::Vector
-    completely_distributed_solution (mpi_communicator,
-                                    dof_handler.n_dofs(),
-                                    dof_handler.n_locally_owned_dofs());
-
-  SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
-
-  PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
-  PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
-
-  solver.solve (system_matrix, completely_distributed_solution, system_rhs,
-               preconditioner);
-
-  pcout << "  Solved in " << solver_control.last_step()
-       << " iterations." << std::endl;
-
-  constraints.distribute (completely_distributed_solution);
-
-  locally_relevant_solution = completely_distributed_solution;
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(3),
-                                     typename FunctionMap<dim>::type(),
-                                     locally_relevant_solution,
-                                     estimated_error_per_cell);
-  parallel::distributed::GridRefinement::
-    refine_and_coarsen_fixed_number (triangulation,
-                                    estimated_error_per_cell,
-                                    0.3, 0.03);
-  triangulation.execute_coarsening_and_refinement ();
-}
-
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
-{
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (locally_relevant_solution, "u");
-
-  Vector<float> subdomain (triangulation.n_active_cells());
-                                  // could just fill entire vector with subdomain_id()
-  {
-    unsigned int index = 0;
-    for (typename Triangulation<dim>::active_cell_iterator
-          cell = triangulation.begin_active();
-        cell != triangulation.end(); ++cell, ++index)
-      subdomain(index) = (cell->is_ghost() || cell->is_artificial()
-                         ?
-                         -1
-                         :
-                         cell->subdomain_id());
-  }
-  data_out.add_data_vector (subdomain, "subdomain");
-  data_out.build_patches ();
-
-  const std::string filename = ("solution-" +
-                               Utilities::int_to_string (cycle, 2) +
-                               "." +
-                               Utilities::int_to_string
-                               (triangulation.locally_owned_subdomain(), 4));
-
-  std::ofstream output ((filename + ".vtu").c_str());
-  data_out.write_vtu (output);
-
-  if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
-    {
-      std::vector<std::string> filenames;
-      for (unsigned int i=0;
-          i<Utilities::System::get_n_mpi_processes(mpi_communicator);
-          ++i)
-       filenames.push_back ("solution-" +
-                            Utilities::int_to_string (cycle, 2) +
-                            "." +
-                            Utilities::int_to_string (i, 4) +
-                            ".vtu");
-
-      std::ofstream master_output ((filename + ".pvtu").c_str());
-      data_out.write_pvtu_record (master_output, filenames);
-    }
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::run (const unsigned int initial_global_refine)
-{
-  const unsigned int n_cycles = 12;
-  for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
-    {
-      pcout << "Cycle " << cycle << ':' << std::endl;
-
-      if (cycle == 0)
-       {
-         GridGenerator::hyper_cube (triangulation);
-         triangulation.refine_global (initial_global_refine);
-       }
-      else
-       refine_grid ();
-
-      pcout << "   Number of active cells:       "
-           << triangulation.n_global_active_cells()
-           << std::endl;
-
-      dof_handler.distribute_dofs (fe);
-
-                                      // eliminate
-      locally_owned_dofs = dof_handler.locally_owned_dofs ();
-      DoFTools::extract_locally_relevant_dofs (dof_handler,
-                                              locally_relevant_dofs);
-
-      pcout << "   Number of degrees of freedom: "
-           << dof_handler.n_dofs()
-           << std::endl;
-
-      setup_system ();
-
-      assemble_system ();
-      solve ();
-
-      if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 100)
-       output_results (cycle);
-
-      pcout << std::endl;
-    }
-}
-
-
-
-int main(int argc, char *argv[])
-{
-  try
-    {
-      PetscInitialize(&argc, &argv);
-      deallog.depth_console (0);
-
-      int refine=5;
-      if (argc>1)
-       {
-         refine = (unsigned int)Utilities::string_to_int(argv[1]);
-       }
-
-      {
-       LaplaceProblem<2> laplace_problem_2d;
-       laplace_problem_2d.run (refine);
-      }
-
-      PetscFinalize();
-    }
-  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;
-}
index b8683b8024b38d9e242946e01c49621ca9de56f0..b94923467cb2ebdf63d2078a7047dbe7e1af6a57 100644 (file)
@@ -1,6 +1,6 @@
 /* $Id$ */
-/* Author: Wolfgang Bangerth, Texas A&M University, 2009 */
-/*         Timo Heister, University of Goettingen, 2009 */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2009, 2010 */
+/*         Timo Heister, University of Goettingen, 2009, 2010 */
 
 /*    $Id$       */
 /*                                                                */
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
-// comment to use PETSc:
-//#define USE_TRILINOS
-
-
 #include <base/quadrature_lib.h>
 #include <base/function.h>
-#include <base/logstream.h>
 #include <base/utilities.h>
 #include <base/conditional_ostream.h>
 #include <base/index_set.h>
 #include <base/timer.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
-#include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
-#include <lac/precondition.h>
 #include <lac/constraint_matrix.h>
 #include <lac/compressed_simple_sparsity_pattern.h>
 #include <lac/sparsity_tools.h>
 
-#ifdef USE_TRILINOS
-#include <lac/trilinos_precondition.h>
-#include <lac/trilinos_sparse_matrix.h>
-#include <lac/trilinos_sparsity_pattern.h>
-#include <lac/trilinos_vector.h>
-#else
 #include <lac/petsc_parallel_sparse_matrix.h>
 #include <lac/petsc_parallel_vector.h>
 #include <lac/petsc_solver.h>
 #include <lac/petsc_precondition.h>
-#endif
 
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <grid/grid_out.h>
-#include <grid/filtered_iterator.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_values.h>
 #include <fe/fe_q.h>
 #include <numerics/vectors.h>
-#include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 
 #include <distributed/tria.h>
 #include <distributed/grid_refinement.h>
-#include <distributed/solution_transfer.h>
 
 #include <fstream>
 #include <iostream>
 
-#include "timeblock.h"
-
 using namespace dealii;
 
-void print_it(Utilities::System::MinMaxAvg & result)
-{
-  std::cout// << "sum: " << result.sum
-    << " avg: " << (long)result.avg/1024
-         << " min: " << (long)result.min/1024 << " @" << result.min_index
-         << " max: " << (long)result.max/1024 << " @" << result.max_index
-         << std::endl;
-}
-
-void print_memory_stats()
-{
-  int myid = Utilities::System::get_this_mpi_process(MPI_COMM_WORLD);
-  Utilities::System::MemoryStats stats;
-  Utilities::System::get_memory_stats(stats);
-  Utilities::System::MinMaxAvg r;
-  Utilities::System::calculate_collective_mpi_min_max_avg(MPI_COMM_WORLD, stats.VmPeak, r);
-  if (myid==0)
-    {
-      std::cout << "MEM: VmPeak: ";
-      print_it(r);
-    }
-  Utilities::System::calculate_collective_mpi_min_max_avg(MPI_COMM_WORLD, stats.VmSize, r);
-  if (myid==0)
-    {
-      std::cout << "MEM: VmSize: ";
-      print_it(r);
-    }
-  Utilities::System::calculate_collective_mpi_min_max_avg(MPI_COMM_WORLD, stats.VmHWM, r);
-  if (myid==0)
-    {
-      std::cout << "MEM: VmHWM:  ";
-      print_it(r);
-    }
-  Utilities::System::calculate_collective_mpi_min_max_avg(MPI_COMM_WORLD, stats.VmRSS, r);
-  if (myid==0)
-    {
-      std::cout << "MEM: VmRSS:  ";
-      print_it(r);
-    }
-}
-
 
 template <int dim>
 class LaplaceProblem
@@ -118,7 +56,7 @@ class LaplaceProblem
     LaplaceProblem ();
     ~LaplaceProblem ();
 
-    void run (unsigned int refine);
+    void run (const unsigned int initial_global_refine);
 
   private:
     void setup_system ();
@@ -127,6 +65,8 @@ class LaplaceProblem
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
 
+    MPI_Comm mpi_communicator;
+
     parallel::distributed::Triangulation<dim>   triangulation;
 
     DoFHandler<dim>      dof_handler;
@@ -137,21 +77,11 @@ class LaplaceProblem
 
     ConstraintMatrix     constraints;
 
-#ifdef USE_TRILINOS
-    TrilinosWrappers::SparseMatrix    system_matrix;
-    TrilinosWrappers::MPI::Vector     solution;
-    TrilinosWrappers::MPI::Vector     system_rhs;
-    parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::Vector> soltrans;
-#else
     PETScWrappers::MPI::SparseMatrix system_matrix;
-    PETScWrappers::MPI::Vector solution;
+    PETScWrappers::MPI::Vector locally_relevant_solution;
     PETScWrappers::MPI::Vector system_rhs;
-    parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector> soltrans;
-#endif
 
     ConditionalOStream                pcout;
-
-
 };
 
 
@@ -159,15 +89,16 @@ class LaplaceProblem
 template <int dim>
 LaplaceProblem<dim>::LaplaceProblem ()
                :
-               triangulation (MPI_COMM_WORLD
-                              ,typename Triangulation<dim>::MeshSmoothing  (Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening)
-               ),
+               mpi_communicator (MPI_COMM_WORLD),
+               triangulation (mpi_communicator,
+                              typename Triangulation<dim>::MeshSmoothing
+                              (Triangulation<dim>::smoothing_on_refinement |
+                               Triangulation<dim>::smoothing_on_coarsening)),
                dof_handler (triangulation),
                 fe (2),
-               soltrans(dof_handler),
                pcout (std::cout,
                       (Utilities::System::
-                       get_this_mpi_process(MPI_COMM_WORLD)
+                       get_this_mpi_process(mpi_communicator)
                        == 0))
 {}
 
@@ -184,82 +115,42 @@ LaplaceProblem<dim>::~LaplaceProblem ()
 template <int dim>
 void LaplaceProblem<dim>::setup_system ()
 {
-  TimeBlock<ConditionalOStream> t(pcout, "setup_system", false);
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "_init_vectors");
-#ifdef USE_TRILINOS
-    solution.reinit (locally_relevant_dofs, MPI_COMM_WORLD);
-    system_rhs.reinit (locally_owned_dofs, MPI_COMM_WORLD);
-#else
-    solution.reinit ( MPI_COMM_WORLD, locally_owned_dofs, locally_relevant_dofs);
-    system_rhs.reinit (MPI_COMM_WORLD, dof_handler.n_dofs(),
-                      dof_handler.n_locally_owned_dofs());
-
-    solution = 0;
-    system_rhs = 0;
-#endif
-  }
+  locally_relevant_solution.reinit (mpi_communicator,
+                                   locally_owned_dofs,
+                                   locally_relevant_dofs);
+  system_rhs.reinit (mpi_communicator,
+                    dof_handler.n_dofs(),
+                    dof_handler.n_locally_owned_dofs());
+
+  locally_relevant_solution = 0;
+  system_rhs = 0;
 
   constraints.clear ();
   constraints.reinit (locally_relevant_dofs);
 
-  {
-
-    TimeBlock<ConditionalOStream> t(pcout, "_make_hn");
-
-    DoFTools::make_hanging_node_constraints (static_cast<const DoFHandler<dim>&>(dof_handler),
-                                            constraints);
-  }
-
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "_interpol_bv");
-    VectorTools::interpolate_boundary_values (static_cast<const DoFHandler<dim>&>(dof_handler),
-                                             0,
-                                             ZeroFunction<dim>(),
-                                             constraints);
-    constraints.close ();
-  }
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close ();
 
 
   CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
                                       dof_handler.n_dofs(),
                                       locally_relevant_dofs);
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "_make_sp", false);
-
-    DoFTools::make_sparsity_pattern (static_cast<const DoFHandler<dim>&>(dof_handler),
-                                    csp,
-                                    constraints, false);
-  }
-
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "_init_matrix", false);
-#ifdef USE_TRILINOS
-    system_matrix.reinit (locally_owned_dofs, csp, MPI_COMM_WORLD, true);
-#else
-
-    {
-      TimeBlock<ConditionalOStream> t(pcout, "__prep_csp", false);
-      SparsityTools::distribute_sparsity_pattern<>(csp,
-                                                  dof_handler.n_locally_owned_dofs_per_processor(),
-                                                  MPI_COMM_WORLD,
-                                                  locally_relevant_dofs);
-    }
-
-    system_matrix.reinit (MPI_COMM_WORLD,
-                         csp,
-                         dof_handler.n_locally_owned_dofs_per_processor(),
-                         dof_handler.n_locally_owned_dofs_per_processor(),
-                         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD));
-#if (PETSC_VERSION_MAJOR <= 2)
-//  MatSetOption (system_matrix, MAT_YES_NEW_NONZERO_LOCATIONS);
-#else
-//  MatSetOption (system_matrix, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);
-//  MatSetOption (system_matrix, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE);
-#endif
-
-#endif
-  }
+  DoFTools::make_sparsity_pattern (dof_handler,
+                                  csp,
+                                  constraints, false);
+  SparsityTools::distribute_sparsity_pattern (csp,
+                                             dof_handler.n_locally_owned_dofs_per_processor(),
+                                             mpi_communicator,
+                                             locally_relevant_dofs);
+  system_matrix.reinit (mpi_communicator,
+                       csp,
+                       dof_handler.n_locally_owned_dofs_per_processor(),
+                       dof_handler.n_locally_owned_dofs_per_processor(),
+                       Utilities::System::get_this_mpi_process(mpi_communicator));
 }
 
 
@@ -267,8 +158,6 @@ void LaplaceProblem<dim>::setup_system ()
 template <int dim>
 void LaplaceProblem<dim>::assemble_system ()
 {
-  TimeBlock<ConditionalOStream> t(pcout, "assemble");
-
   const QGauss<dim>  quadrature_formula(3);
 
   FEValues<dim> fe_values (fe, quadrature_formula,
@@ -303,17 +192,15 @@ void LaplaceProblem<dim>::assemble_system ()
                                     fe_values.shape_grad(j,q_point) *
                                     fe_values.JxW(q_point));
 
-             double v = fe_values.quadrature_point(q_point)[0];
-             v = 0.5+0.25*sin(4.0*numbers::PI*v);
-
              cell_rhs(i) += (fe_values.shape_value(i,q_point) *
-                             (v < fe_values.quadrature_point(q_point)[1]
+                             (fe_values.quadrature_point(q_point)[1]
+                              >
+                              0.5+0.25*sin(4.0*numbers::PI*fe_values.quadrature_point(q_point)[0])
                               ? 1 : -1) *
                              fe_values.JxW(q_point));
            }
 
        cell->get_dof_indices (local_dof_indices);
-
        constraints.distribute_local_to_global (cell_matrix,
                                                cell_rhs,
                                                local_dof_indices,
@@ -321,12 +208,8 @@ void LaplaceProblem<dim>::assemble_system ()
                                                system_rhs);
       }
 
-#ifdef USE_TRILINOS
-  system_rhs.compress (Add);
-#else
-  system_rhs.compress ();
-#endif
   system_matrix.compress ();
+  system_rhs.compress ();
 }
 
 
@@ -336,83 +219,25 @@ void LaplaceProblem<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::solve ()
 {
-                                  // we do not want to get spammed with solver info.
-  deallog.depth_console (0);
-
-  TimeBlock<ConditionalOStream> t(pcout, "solve", false);
-
-#ifdef USE_TRILINOS
-  TrilinosWrappers::MPI::Vector distributed_solution (locally_owned_dofs,
-                                                     MPI_COMM_WORLD);
-
-  distributed_solution.reinit(solution, false, true);
-#else
+  PETScWrappers::MPI::Vector
+    completely_distributed_solution (mpi_communicator,
+                                    dof_handler.n_dofs(),
+                                    dof_handler.n_locally_owned_dofs());
 
-  PETScWrappers::MPI::Vector distributed_solution(MPI_COMM_WORLD,
-                                                 dof_handler.n_dofs(),
-                                                 dof_handler.n_locally_owned_dofs());
+  SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
 
-  distributed_solution = solution;
-#endif
-
-  print_memory_stats();
-
-  SolverControl solver_control (solution.size(), 1e-12);
-#ifdef USE_TRILINOS
-  SolverCG<TrilinosWrappers::MPI::Vector> solver (solver_control);
-
-  TrilinosWrappers::PreconditionAMG preconditioner;
-  preconditioner.initialize(system_matrix);
-#else
-  PETScWrappers::SolverCG solver(solver_control, MPI_COMM_WORLD);
+  PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
   PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
 
-#endif
-
-  { // memory consumption
-    pcout << "Mem: Tria (p4est) DofH Constraints Mat X Rhs I_lo I_lr"
-         << std::endl;
-//    MPI_Barrier(MPI_COMM_WORLD);
-
-    pcout << "MEM: "
-         << " " << triangulation.memory_consumption()
-         << " " << triangulation.memory_consumption_p4est()
-         << " " << dof_handler.memory_consumption()
-         << " " << constraints.memory_consumption()
-         << " " << system_matrix.memory_consumption()
-         << " " << solution.memory_consumption()
-         << " " << system_rhs.memory_consumption()
-         << " " << locally_owned_dofs.memory_consumption()
-         << " " << locally_relevant_dofs.memory_consumption()
-//       << " " << preconditioner.memory_consumption()
-         << " ~sum=" << (triangulation.memory_consumption() + dof_handler.memory_consumption()
-                         + constraints.memory_consumption() + system_matrix.memory_consumption()
-                         + solution.memory_consumption()  + system_rhs.memory_consumption()
-                         + locally_owned_dofs.memory_consumption()
-                         + locally_relevant_dofs.memory_consumption() )/1024
-         << " kB"
-         << std::endl;
-
-//    MPI_Barrier(MPI_COMM_WORLD);
-  }
-
-
-  solver.solve (system_matrix, distributed_solution, system_rhs,
-           preconditioner);
-
-  print_memory_stats();
+  solver.solve (system_matrix, completely_distributed_solution, system_rhs,
+               preconditioner);
 
   pcout << "  Solved in " << solver_control.last_step()
        << " iterations." << std::endl;
 
-  constraints.distribute (distributed_solution);
-  solution = distributed_solution;
-
-#ifndef USE_TRILINOS
-  solution.update_ghost_values();
-#endif
+  constraints.distribute (completely_distributed_solution);
 
-  deallog.depth_console (1);
+  locally_relevant_solution = completely_distributed_solution;
 }
 
 
@@ -420,104 +245,17 @@ void LaplaceProblem<dim>::solve ()
 template <int dim>
 void LaplaceProblem<dim>::refine_grid ()
 {
-  TimerOutput computing_timer(pcout, TimerOutput::summary,
-                             TimerOutput::wall_times);
-
-  computing_timer.enter_section ("Estimating");
-
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  {
-
-    TimeBlock<ConditionalOStream> t(pcout, "kelly");
-
-  KellyErrorEstimator<dim>::estimate (static_cast<const DoFHandler<dim>&>(dof_handler),
+  KellyErrorEstimator<dim>::estimate (dof_handler,
                                      QGauss<dim-1>(3),
                                      typename FunctionMap<dim>::type(),
-                                     solution,
+                                     locally_relevant_solution,
                                      estimated_error_per_cell);
-  }
-
-  computing_timer.exit_section();
-
-  computing_timer.enter_section ("Marking");
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "marking");
-
   parallel::distributed::GridRefinement::
     refine_and_coarsen_fixed_number (triangulation,
                                     estimated_error_per_cell,
                                     0.3, 0.03);
-  }
-
-  computing_timer.exit_section();
-
-  computing_timer.enter_section ("Refining");
-
-
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "prep_c&r");
-
-    triangulation.prepare_coarsening_and_refinement();
-  }
-
-  if (0)
-    {
-                                      //output refinement information for
-                                      //debugging
-
-     const std::string filename = ("ref." +
-                                  Utilities::int_to_string
-                                  (triangulation.locally_owned_subdomain(), 4) +
-                                  ".vtk");
-     std::ofstream output (filename.c_str());
-
-     DataOut<dim> data_out;
-     data_out.attach_dof_handler (dof_handler);
-
-     unsigned int n_coarse=0;
-     Vector<float> subdomain (triangulation.n_active_cells());
-     {
-       unsigned int index = 0;
-
-       for (typename Triangulation<dim>::active_cell_iterator
-             cell = triangulation.begin_active();
-           cell != triangulation.end(); ++cell, ++index)
-        {
-          subdomain(index)=0;
-
-          if (cell->is_ghost() || cell->is_artificial())
-            subdomain(index)=-4;
-
-          if (cell->refine_flag_set())
-            subdomain(index)+=1;
-          if (cell->coarsen_flag_set())
-            {
-              subdomain(index)+=2;
-              ++n_coarse;
-            }
-        }
-     }
-     std::cout << "id=" << triangulation.locally_owned_subdomain() << " n_coarsen=" << n_coarse << std::endl;
-
-     data_out.add_data_vector (subdomain, "info");
-     data_out.build_patches ();
-     data_out.write_vtk (output);
-    }
-
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "soltrans.prep");
-
-  soltrans.prepare_for_coarsening_and_refinement(solution);
-  }
-
-  {
-    TimeBlock<ConditionalOStream> t(pcout, "c&r", false);
-
   triangulation.execute_coarsening_and_refinement ();
-  }
-
-  computing_timer.exit_section();
 }
 
 
@@ -526,19 +264,12 @@ void LaplaceProblem<dim>::refine_grid ()
 template <int dim>
 void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 {
-  const std::string filename = ("solution-" +
-                               Utilities::int_to_string (cycle, 2) +
-                               "." +
-                               Utilities::int_to_string
-                               (triangulation.locally_owned_subdomain(), 4) +
-                               ".d2");
-  std::ofstream output (filename.c_str());
-
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution, "u");
+  data_out.add_data_vector (locally_relevant_solution, "u");
 
   Vector<float> subdomain (triangulation.n_active_cells());
+                                  // could just fill entire vector with subdomain_id()
   {
     unsigned int index = 0;
     for (typename Triangulation<dim>::active_cell_iterator
@@ -551,67 +282,59 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
                          cell->subdomain_id());
   }
   data_out.add_data_vector (subdomain, "subdomain");
-
   data_out.build_patches ();
 
-  data_out.write_deal_II_intermediate (output);
+  const std::string filename = ("solution-" +
+                               Utilities::int_to_string (cycle, 2) +
+                               "." +
+                               Utilities::int_to_string
+                               (triangulation.locally_owned_subdomain(), 4));
+
+  std::ofstream output ((filename + ".vtu").c_str());
+  data_out.write_vtu (output);
+
+  if (Utilities::System::get_this_mpi_process(mpi_communicator) == 0)
+    {
+      std::vector<std::string> filenames;
+      for (unsigned int i=0;
+          i<Utilities::System::get_n_mpi_processes(mpi_communicator);
+          ++i)
+       filenames.push_back ("solution-" +
+                            Utilities::int_to_string (cycle, 2) +
+                            "." +
+                            Utilities::int_to_string (i, 4) +
+                            ".vtu");
+
+      std::ofstream master_output ((filename + ".pvtu").c_str());
+      data_out.write_pvtu_record (master_output, filenames);
+    }
 }
 
 
 
 template <int dim>
-void LaplaceProblem<dim>::run (unsigned int refine)
+void LaplaceProblem<dim>::run (const unsigned int initial_global_refine)
 {
-      {
-       TimeBlock<ConditionalOStream> t(pcout, "empty");
-      }
-
-      unsigned int cycles = 12;
-
-      pcout << "running " <<  cycles << " cycles with refine= " << refine <<
-#ifdef USE_TRILINOS
-           " with Trilinos"
-#else
-           " with PETSc"
-#endif
-           << std::endl;
-
-  for (unsigned int cycle=0; cycle<cycles; ++cycle)
+  const unsigned int n_cycles = 12;
+  for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
     {
-      TimerOutput computing_timer(pcout, TimerOutput::summary,
-                                 TimerOutput::wall_times);
-
       pcout << "Cycle " << cycle << ':' << std::endl;
 
-      computing_timer.enter_section ("Mesh handling");
       if (cycle == 0)
        {
          GridGenerator::hyper_cube (triangulation);
-         triangulation.refine_global (refine);
+         triangulation.refine_global (initial_global_refine);
        }
       else
        refine_grid ();
 
-
-
       pcout << "   Number of active cells:       "
            << triangulation.n_global_active_cells()
            << std::endl;
-/*      pcout << "                                 (";
-      for (unsigned int i=0; i<triangulation.n_active_cells().size(); ++i)
-       pcout << triangulation.n_active_cells()[i]
-             << (i != triangulation.n_active_cells().size()-1 ? "+" : "");
-             pcout << ")" << std::endl;*/
-      computing_timer.exit_section();
-
-      computing_timer.enter_section ("DoFs handling");
-
-      {
-       TimeBlock<ConditionalOStream> t(pcout, "dist_dofs");
-       dof_handler.distribute_dofs (fe);
-      }
 
+      dof_handler.distribute_dofs (fe);
 
+                                      // eliminate
       locally_owned_dofs = dof_handler.locally_owned_dofs ();
       DoFTools::extract_locally_relevant_dofs (dof_handler,
                                               locally_relevant_dofs);
@@ -619,35 +342,14 @@ void LaplaceProblem<dim>::run (unsigned int refine)
       pcout << "   Number of degrees of freedom: "
            << dof_handler.n_dofs()
            << std::endl;
-/*      pcout << "                                 (";
-      for (unsigned int i=0; i<dof_handler.n_dofs().size(); ++i)
-       pcout << dof_handler.n_dofs()[i]
-             << (i != dof_handler.n_dofs().size()-1 ? "+" : "");
-             pcout << ")" << std::endl;*/
-      computing_timer.exit_section();
-
-      computing_timer.enter_section ("Setting up system");
-      setup_system ();
-
-      if (cycle>0)
-       {
-         TimeBlock<ConditionalOStream> t(pcout, "soltrans.interp");
 
-       soltrans.interpolate(solution);
-       }
-
-
-      computing_timer.exit_section();
+      setup_system ();
 
-      computing_timer.enter_section ("Assembling system");
       assemble_system ();
-      computing_timer.exit_section();
-
-      computing_timer.enter_section ("Solving system");
       solve ();
-      computing_timer.exit_section();
 
-//      output_results (cycle);
+      if (Utilities::System::get_n_mpi_processes(mpi_communicator) <= 100)
+       output_results (cycle);
 
       pcout << std::endl;
     }
@@ -659,13 +361,8 @@ int main(int argc, char *argv[])
 {
   try
     {
-#ifdef USE_TRILINOS
-      Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
-#else
       PetscInitialize(&argc, &argv);
-#endif
-
-      deallog.depth_console (1);
+      deallog.depth_console (0);
 
       int refine=5;
       if (argc>1)
@@ -678,13 +375,7 @@ int main(int argc, char *argv[])
        laplace_problem_2d.run (refine);
       }
 
-      print_memory_stats();
-
-#ifndef USE_TRILINOS
-      PetscLogPrintSummary(MPI_COMM_WORLD, "petsc.log" );
       PetscFinalize();
-#endif
-
     }
   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.