]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add additional tests (simplied version of matrix_free/step-48) to see what is going on
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 16 Apr 2014 19:30:38 +0000 (19:30 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 16 Apr 2014 19:30:38 +0000 (19:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@32783 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-48/step-48.cc
tests/matrix_free/step-48.cc
tests/matrix_free/step-48b.cc [new file with mode: 0644]
tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output [new file with mode: 0644]
tests/matrix_free/step-48c.cc [new file with mode: 0644]
tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output [new file with mode: 0644]

index 95d9d67d9d171c448944f86373ddfd51d32dae68..77e1f719703b4bdb3b93cec89b6ac6d8e5ac8952 100644 (file)
@@ -286,8 +286,6 @@ namespace Step48
     ConditionalOStream pcout;
 
     void make_grid_and_dofs ();
-    void oldstyle_operation ();
-    void assemble_system ();
     void output_results (const unsigned int timestep_number);
 
 #ifdef DEAL_II_WITH_P4EST
index 3add3fe5d90c79ec5fef54f6cbae643d465acd82..b920f8b52c49d19621ff524a07dfc629805d532f 100644 (file)
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2012 - 2013 by the deal.II authors
+// Copyright (C) 2012 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -207,8 +207,6 @@ namespace Step48
 
   private:
     void make_grid_and_dofs ();
-    void oldstyle_operation ();
-    void assemble_system ();
     void output_results (const unsigned int timestep_number);
 
 #ifdef DEAL_II_WITH_P4EST
diff --git a/tests/matrix_free/step-48b.cc b/tests/matrix_free/step-48b.cc
new file mode 100644 (file)
index 0000000..11cb3c2
--- /dev/null
@@ -0,0 +1,412 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2012 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// simplified version of step-48, no sine term (i.e., usual wave equation)
+
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+
+namespace Step48
+{
+  using namespace dealii;
+
+  const unsigned int dimension = 2;
+  const unsigned int fe_degree = 4;
+
+
+
+  template <int dim, int fe_degree>
+  class SineGordonOperation
+  {
+  public:
+    SineGordonOperation(const MatrixFree<dim,double> &data_in,
+                        const double                      time_step);
+
+    void apply (parallel::distributed::Vector<double>                        &dst,
+                const std::vector<parallel::distributed::Vector<double>*> &src) const;
+
+  private:
+    const MatrixFree<dim,double>         &data;
+    const VectorizedArray<double>         delta_t_sqr;
+    parallel::distributed::Vector<double> inv_mass_matrix;
+
+    void local_apply (const MatrixFree<dim,double>               &data,
+                      parallel::distributed::Vector<double>      &dst,
+                      const std::vector<parallel::distributed::Vector<double>*> &src,
+                      const std::pair<unsigned int,unsigned int> &cell_range) const;
+  };
+
+
+
+
+  template <int dim, int fe_degree>
+  SineGordonOperation<dim,fe_degree>::
+  SineGordonOperation(const MatrixFree<dim,double> &data_in,
+                      const double                  time_step)
+    :
+    data(data_in),
+    delta_t_sqr(make_vectorized_array(time_step *time_step))
+  {
+    VectorizedArray<double> one = make_vectorized_array (1.);
+
+    data.initialize_dof_vector (inv_mass_matrix);
+
+    FEEvaluationGL<dim,fe_degree> fe_eval(data);
+    const unsigned int            n_q_points = fe_eval.n_q_points;
+
+    for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
+      {
+        fe_eval.reinit(cell);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          fe_eval.submit_value(one,q);
+        fe_eval.integrate (true,false);
+        fe_eval.distribute_local_to_global (inv_mass_matrix);
+      }
+
+    inv_mass_matrix.compress(VectorOperation::add);
+    for (unsigned int k=0; k<inv_mass_matrix.local_size(); ++k)
+      if (inv_mass_matrix.local_element(k)>1e-15)
+        inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k);
+      else
+        inv_mass_matrix.local_element(k) = 0;
+  }
+
+
+
+
+
+  template <int dim, int fe_degree>
+  void SineGordonOperation<dim, fe_degree>::
+  local_apply (const MatrixFree<dim>                 &data,
+               parallel::distributed::Vector<double>     &dst,
+               const std::vector<parallel::distributed::Vector<double>*> &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    AssertDimension (src.size(), 2);
+    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        current.reinit (cell);
+        old.reinit (cell);
+
+        current.read_dof_values (*src[0]);
+        old.read_dof_values     (*src[1]);
+
+        current.evaluate (true, true, false);
+        old.evaluate (true, false, false);
+
+        for (unsigned int q=0; q<current.n_q_points; ++q)
+          {
+            const VectorizedArray<double> current_value = current.get_value(q);
+            const VectorizedArray<double> old_value     = old.get_value(q);
+
+            current.submit_value (2.*current_value - old_value -
+                                  delta_t_sqr * std::sin(current_value),q);
+            current.submit_gradient (- delta_t_sqr *
+                                     current.get_gradient(q), q);
+          }
+
+        current.integrate (true,true);
+        current.distribute_local_to_global (dst);
+      }
+  }
+
+
+
+
+  template <int dim, int fe_degree>
+  void SineGordonOperation<dim, fe_degree>::
+  apply (parallel::distributed::Vector<double>                        &dst,
+         const std::vector<parallel::distributed::Vector<double>*> &src) const
+  {
+    dst = 0;
+    data.cell_loop (&SineGordonOperation<dim,fe_degree>::local_apply,
+                    this, dst, src);
+    dst.scale(inv_mass_matrix);
+  }
+
+
+
+  template <int dim>
+  class ExactSolution : public Function<dim>
+  {
+  public:
+    ExactSolution (const unsigned int n_components = 1,
+                   const double time = 0.) : Function<dim>(n_components, time) {}
+    virtual double value (const Point<dim> &p,
+                          const unsigned int component = 0) const;
+  };
+
+  template <int dim>
+  double ExactSolution<dim>::value (const Point<dim> &p,
+                                    const unsigned int /* component */) const
+  {
+    double t = this->get_time ();
+
+    const double m = 0.5;
+    const double c1 = 0.;
+    const double c2 = 0.;
+    const double factor = (m / std::sqrt(1.-m*m) *
+                           std::sin(std::sqrt(1.-m*m)*t+c2));
+    double result = 1.;
+    for (unsigned int d=0; d<dim; ++d)
+      result *= -4. * std::atan (factor / std::cosh(m*p[d]+c1));
+    return result;
+  }
+
+
+
+
+  template <int dim>
+  class SineGordonProblem
+  {
+  public:
+    SineGordonProblem ();
+    void run ();
+
+  private:
+    void make_grid_and_dofs ();
+    void output_results (const unsigned int timestep_number);
+
+#ifdef DEAL_II_WITH_P4EST
+    parallel::distributed::Triangulation<dim>   triangulation;
+#else
+    Triangulation<dim>   triangulation;
+#endif
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+    ConstraintMatrix     constraints;
+    IndexSet             locally_relevant_dofs;
+
+    MatrixFree<dim,double> matrix_free_data;
+
+    parallel::distributed::Vector<double> solution, old_solution, old_old_solution;
+
+    const unsigned int n_global_refinements;
+    double time, time_step;
+    const double final_time;
+    const double cfl_number;
+    const unsigned int output_timestep_skip;
+  };
+
+
+
+  template <int dim>
+  SineGordonProblem<dim>::SineGordonProblem ()
+    :
+#ifdef DEAL_II_WITH_P4EST
+    triangulation (MPI_COMM_WORLD),
+#endif
+    fe (QGaussLobatto<1>(fe_degree+1)),
+    dof_handler (triangulation),
+    n_global_refinements (8-2*dim),
+    time (-10),
+    final_time (-9),
+    cfl_number (.1/fe_degree),
+    output_timestep_skip (200)
+  {}
+
+
+  template <int dim>
+  void SineGordonProblem<dim>::make_grid_and_dofs ()
+  {
+    GridGenerator::hyper_cube (triangulation, -15, 15);
+    triangulation.refine_global (n_global_refinements);
+
+    deallog << "   Number of global active cells: "
+#ifdef DEAL_II_WITH_P4EST
+            << triangulation.n_global_active_cells()
+#else
+            << triangulation.n_active_cells()
+#endif
+            << std::endl;
+
+    dof_handler.distribute_dofs (fe);
+
+    deallog << "   Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << std::endl;
+
+
+    DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                             locally_relevant_dofs);
+    constraints.close();
+
+    QGaussLobatto<1> quadrature (fe_degree+1);
+    typename MatrixFree<dim>::AdditionalData additional_data;
+    additional_data.mpi_communicator = MPI_COMM_WORLD;
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim>::AdditionalData::partition_partition;
+
+    matrix_free_data.reinit (dof_handler, constraints,
+                             quadrature, additional_data);
+
+    matrix_free_data.initialize_dof_vector (solution);
+    old_solution.reinit (solution);
+    old_old_solution.reinit (solution);
+  }
+
+
+
+
+  template <int dim>
+  void
+  SineGordonProblem<dim>::output_results (const unsigned int timestep_number)
+  {
+    constraints.distribute (solution);
+    solution.update_ghost_values();
+
+    Vector<float> norm_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (dof_handler,
+                                       solution,
+                                       ZeroFunction<dim>(),
+                                       norm_per_cell,
+                                       QGauss<dim>(fe_degree+1),
+                                       VectorTools::L2_norm);
+    const double solution_norm =
+      std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD));
+
+    deallog << "   Time:"
+            << std::setw(8) << std::setprecision(3) << time
+            << ", solution norm: "
+            << std::setprecision(5) << std::setw(7) << solution_norm
+            << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  SineGordonProblem<dim>::run ()
+  {
+    make_grid_and_dofs();
+
+    const double local_min_cell_diameter =
+      triangulation.last()->diameter()/std::sqrt(dim);
+    const double global_min_cell_diameter
+      = -Utilities::MPI::max(-local_min_cell_diameter, MPI_COMM_WORLD);
+    time_step = cfl_number * global_min_cell_diameter;
+    time_step = (final_time-time)/(int((final_time-time)/time_step));
+    deallog << "   Time step size: " << time_step << ", finest cell: "
+            << global_min_cell_diameter << std::endl << std::endl;
+
+
+    VectorTools::interpolate (dof_handler,
+                              ExactSolution<dim> (1, time),
+                              solution);
+    VectorTools::interpolate (dof_handler,
+                              ExactSolution<dim> (1, time-time_step),
+                              old_solution);
+    output_results (0);
+
+    std::vector<parallel::distributed::Vector<double>*> previous_solutions;
+    previous_solutions.push_back(&old_solution);
+    previous_solutions.push_back(&old_old_solution);
+
+    SineGordonOperation<dim,fe_degree> sine_gordon_op (matrix_free_data,
+                                                       time_step);
+
+    unsigned int timestep_number = 1;
+
+    for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number)
+      {
+        old_old_solution.swap (old_solution);
+        old_solution.swap (solution);
+        sine_gordon_op.apply (solution, previous_solutions);
+
+        if (timestep_number % output_timestep_skip == 0)
+          output_results(timestep_number / output_timestep_skip);
+      }
+    output_results(timestep_number / output_timestep_skip + 1);
+
+    deallog << std::endl
+            << "   Performed " << timestep_number << " time steps."
+            << std::endl << std::endl;
+  }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      {
+        deallog.push("2d");
+        Step48::SineGordonProblem<2> sg_problem;
+        sg_problem.run ();
+        deallog.pop();
+      }
+      {
+        deallog.push("3d");
+        Step48::SineGordonProblem<3> sg_problem;
+        sg_problem.run ();
+        deallog.pop();
+      }
+    }
+  else
+    {
+      deallog.depth_console(0);
+      {
+        Step48::SineGordonProblem<2> sg_problem;
+        sg_problem.run ();
+      }
+      {
+        Step48::SineGordonProblem<3> sg_problem;
+        sg_problem.run ();
+      }
+    }
+}
+
diff --git a/tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output b/tests/matrix_free/step-48b.with_mpi=true.mpirun=4.output
new file mode 100644 (file)
index 0000000..5cf259b
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:0:2d::   Number of global active cells: 256
+DEAL:0:2d::   Number of degrees of freedom: 4225
+DEAL:0:2d::   Time step size: 0.04762, finest cell: 1.875
+DEAL:0:2d::
+DEAL:0:2d::   Time:-10.0   , solution norm: 9.5599 
+DEAL:0:2d::   Time:-8.95   , solution norm: 21.848 
+DEAL:0:2d::
+DEAL:0:2d::   Performed 22 time steps.
+DEAL:0:2d::
+DEAL:0:3d::   Number of global active cells: 64
+DEAL:0:3d::   Number of degrees of freedom: 4913
+DEAL:0:3d::   Time step size: 0.20000, finest cell: 7.5000
+DEAL:0:3d::
+DEAL:0:3d::   Time:-10.0   , solution norm: 29.779 
+DEAL:0:3d::   Time:-8.80   , solution norm: 84.725 
+DEAL:0:3d::
+DEAL:0:3d::   Performed 6 time steps.
+DEAL:0:3d::
diff --git a/tests/matrix_free/step-48c.cc b/tests/matrix_free/step-48c.cc
new file mode 100644 (file)
index 0000000..c14b9eb
--- /dev/null
@@ -0,0 +1,404 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2012 - 2014 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.
+//
+// ---------------------------------------------------------------------
+
+
+// simplified form for step-48 test 
+
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/lac/parallel_vector.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+
+namespace Step48
+{
+  using namespace dealii;
+
+  const unsigned int dimension = 2;
+  const unsigned int fe_degree = 4;
+
+
+
+  template <int dim, int fe_degree>
+  class SineGordonOperation
+  {
+  public:
+    SineGordonOperation(const MatrixFree<dim,double> &data_in,
+                        const double                      time_step);
+
+    void apply (parallel::distributed::Vector<double>                        &dst,
+                const std::vector<parallel::distributed::Vector<double>*> &src) const;
+
+  private:
+    const MatrixFree<dim,double>         &data;
+    const VectorizedArray<double>         delta_t_sqr;
+    parallel::distributed::Vector<double> inv_mass_matrix;
+
+    void local_apply (const MatrixFree<dim,double>               &data,
+                      parallel::distributed::Vector<double>      &dst,
+                      const std::vector<parallel::distributed::Vector<double>*> &src,
+                      const std::pair<unsigned int,unsigned int> &cell_range) const;
+  };
+
+
+
+
+  template <int dim, int fe_degree>
+  SineGordonOperation<dim,fe_degree>::
+  SineGordonOperation(const MatrixFree<dim,double> &data_in,
+                      const double                  time_step)
+    :
+    data(data_in),
+    delta_t_sqr(make_vectorized_array(time_step *time_step))
+  {
+    VectorizedArray<double> one = make_vectorized_array (1.);
+
+    data.initialize_dof_vector (inv_mass_matrix);
+
+    FEEvaluationGL<dim,fe_degree> fe_eval(data);
+    const unsigned int            n_q_points = fe_eval.n_q_points;
+
+    for (unsigned int cell=0; cell<data.get_size_info().n_macro_cells; ++cell)
+      {
+        fe_eval.reinit(cell);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          fe_eval.submit_value(one,q);
+        fe_eval.integrate (true,false);
+        fe_eval.distribute_local_to_global (inv_mass_matrix);
+      }
+
+    inv_mass_matrix.compress(VectorOperation::add);
+    for (unsigned int k=0; k<inv_mass_matrix.local_size(); ++k)
+      if (inv_mass_matrix.local_element(k)>1e-15)
+        inv_mass_matrix.local_element(k) = 1./inv_mass_matrix.local_element(k);
+      else
+        inv_mass_matrix.local_element(k) = 0;
+  }
+
+
+
+
+
+  template <int dim, int fe_degree>
+  void SineGordonOperation<dim, fe_degree>::
+  local_apply (const MatrixFree<dim>                 &data,
+               parallel::distributed::Vector<double>     &dst,
+               const std::vector<parallel::distributed::Vector<double>*> &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    AssertDimension (src.size(), 2);
+    FEEvaluationGL<dim,fe_degree> current (data), old (data);
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        current.reinit (cell);
+        old.reinit (cell);
+
+        current.read_dof_values (*src[0]);
+        old.read_dof_values     (*src[1]);
+
+        current.evaluate (true, true, false);
+        old.evaluate (true, false, false);
+
+        for (unsigned int q=0; q<current.n_q_points; ++q)
+          {
+            const VectorizedArray<double> current_value = current.get_value(q);
+            const VectorizedArray<double> old_value     = old.get_value(q);
+
+            current.submit_value (2.*current_value - old_value,q);
+            current.submit_gradient (- delta_t_sqr *
+                                     current.get_gradient(q), q);
+          }
+
+        current.integrate (true,true);
+        current.distribute_local_to_global (dst);
+      }
+  }
+
+
+
+
+  template <int dim, int fe_degree>
+  void SineGordonOperation<dim, fe_degree>::
+  apply (parallel::distributed::Vector<double>                        &dst,
+         const std::vector<parallel::distributed::Vector<double>*> &src) const
+  {
+    dst = 0;
+    data.cell_loop (&SineGordonOperation<dim,fe_degree>::local_apply,
+                    this, dst, src);
+    dst.scale(inv_mass_matrix);
+  }
+
+
+
+  template <int dim>
+  class InitialSolution : public Function<dim>
+  {
+  public:
+    InitialSolution (const unsigned int n_components = 1,
+                   const double time = 0.) : Function<dim>(n_components, time) {}
+    virtual double value (const Point<dim> &p,
+                          const unsigned int component = 0) const;
+  };
+
+  template <int dim>
+  double InitialSolution<dim>::value (const Point<dim> &p,
+                                      const unsigned int /* component */) const
+  {
+    return 4.*std::exp(-p.square()*10);
+  }
+
+
+
+
+  template <int dim>
+  class SineGordonProblem
+  {
+  public:
+    SineGordonProblem ();
+    void run ();
+
+  private:
+    void make_grid_and_dofs ();
+    void output_results (const unsigned int timestep_number);
+
+#ifdef DEAL_II_WITH_P4EST
+    parallel::distributed::Triangulation<dim>   triangulation;
+#else
+    Triangulation<dim>   triangulation;
+#endif
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+    ConstraintMatrix     constraints;
+    IndexSet             locally_relevant_dofs;
+
+    MatrixFree<dim,double> matrix_free_data;
+
+    parallel::distributed::Vector<double> solution, old_solution, old_old_solution;
+
+    const unsigned int n_global_refinements;
+    double time, time_step;
+    const double final_time;
+    const double cfl_number;
+    const unsigned int output_timestep_skip;
+  };
+
+
+
+  template <int dim>
+  SineGordonProblem<dim>::SineGordonProblem ()
+    :
+#ifdef DEAL_II_WITH_P4EST
+    triangulation (MPI_COMM_WORLD),
+#endif
+    fe (QGaussLobatto<1>(fe_degree+1)),
+    dof_handler (triangulation),
+    n_global_refinements (9-2*dim),
+    time (-10),
+    final_time (-9),
+    cfl_number (.1/fe_degree),
+    output_timestep_skip (200)
+  {}
+
+
+  template <int dim>
+  void SineGordonProblem<dim>::make_grid_and_dofs ()
+  {
+    GridGenerator::hyper_cube (triangulation, -15, 15);
+    triangulation.refine_global (n_global_refinements);
+
+    deallog << "   Number of global active cells: "
+#ifdef DEAL_II_WITH_P4EST
+            << triangulation.n_global_active_cells()
+#else
+            << triangulation.n_active_cells()
+#endif
+            << std::endl;
+
+    dof_handler.distribute_dofs (fe);
+
+    deallog << "   Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << std::endl;
+
+
+    DoFTools::extract_locally_relevant_dofs (dof_handler,
+                                             locally_relevant_dofs);
+    constraints.clear();
+    constraints.reinit (locally_relevant_dofs);
+    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+    constraints.close();
+
+    QGaussLobatto<1> quadrature (fe_degree+1);
+    typename MatrixFree<dim>::AdditionalData additional_data;
+    additional_data.mpi_communicator = MPI_COMM_WORLD;
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim>::AdditionalData::partition_partition;
+
+    matrix_free_data.reinit (dof_handler, constraints,
+                             quadrature, additional_data);
+
+    matrix_free_data.initialize_dof_vector (solution);
+    old_solution.reinit (solution);
+    old_old_solution.reinit (solution);
+  }
+
+
+
+
+  template <int dim>
+  void
+  SineGordonProblem<dim>::output_results (const unsigned int timestep_number)
+  {
+    constraints.distribute (solution);
+    solution.update_ghost_values();
+
+    Vector<float> norm_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (dof_handler,
+                                       solution,
+                                       ZeroFunction<dim>(),
+                                       norm_per_cell,
+                                       QGauss<dim>(fe_degree+1),
+                                       VectorTools::L2_norm);
+    const double solution_norm =
+      std::sqrt(Utilities::MPI::sum (norm_per_cell.norm_sqr(), MPI_COMM_WORLD));
+
+    deallog << "   Time:"
+            << std::setw(8) << std::setprecision(3) << time
+            << ", solution norm: "
+            << std::setprecision(5) << std::setw(7) << solution_norm
+            << std::endl;
+  }
+
+
+
+  template <int dim>
+  void
+  SineGordonProblem<dim>::run ()
+  {
+    make_grid_and_dofs();
+
+    const double local_min_cell_diameter =
+      triangulation.last()->diameter()/std::sqrt(dim);
+    const double global_min_cell_diameter
+      = -Utilities::MPI::max(-local_min_cell_diameter, MPI_COMM_WORLD);
+    time_step = cfl_number * global_min_cell_diameter;
+    time_step = (final_time-time)/(int((final_time-time)/time_step));
+    deallog << "   Time step size: " << time_step << ", finest cell: "
+            << global_min_cell_diameter << std::endl << std::endl;
+
+
+    VectorTools::interpolate (dof_handler,
+                              InitialSolution<dim> (1, time),
+                              solution);
+    VectorTools::interpolate (dof_handler,
+                              InitialSolution<dim> (1, time-time_step),
+                              old_solution);
+    output_results (0);
+
+    std::vector<parallel::distributed::Vector<double>*> previous_solutions;
+    previous_solutions.push_back(&old_solution);
+    previous_solutions.push_back(&old_old_solution);
+
+    SineGordonOperation<dim,fe_degree> sine_gordon_op (matrix_free_data,
+                                                       time_step);
+
+    unsigned int timestep_number = 1;
+
+    for (time+=time_step; time<=final_time; time+=time_step, ++timestep_number)
+      {
+        old_old_solution.swap (old_solution);
+        old_solution.swap (solution);
+        sine_gordon_op.apply (solution, previous_solutions);
+
+        if (timestep_number % output_timestep_skip == 0)
+          output_results(timestep_number / output_timestep_skip);
+      }
+    output_results(timestep_number / output_timestep_skip + 1);
+
+    deallog << std::endl
+            << "   Performed " << timestep_number << " time steps."
+            << std::endl << std::endl;
+  }
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      {
+        deallog.push("2d");
+        Step48::SineGordonProblem<2> sg_problem;
+        sg_problem.run ();
+        deallog.pop();
+      }
+      {
+        deallog.push("3d");
+        Step48::SineGordonProblem<3> sg_problem;
+        sg_problem.run ();
+        deallog.pop();
+      }
+    }
+  else
+    {
+      deallog.depth_console(0);
+      {
+        Step48::SineGordonProblem<2> sg_problem;
+        sg_problem.run ();
+      }
+      {
+        Step48::SineGordonProblem<3> sg_problem;
+        sg_problem.run ();
+      }
+    }
+}
+
diff --git a/tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output b/tests/matrix_free/step-48c.with_mpi=true.mpirun=4.output
new file mode 100644 (file)
index 0000000..9639684
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:0:2d::   Number of global active cells: 1024
+DEAL:0:2d::   Number of degrees of freedom: 16641
+DEAL:0:2d::   Time step size: 0.02381, finest cell: 0.9375
+DEAL:0:2d::
+DEAL:0:2d::   Time:-10.0   , solution norm: 1.6062 
+DEAL:0:2d::   Time:-8.98   , solution norm: 1.1332 
+DEAL:0:2d::
+DEAL:0:2d::   Performed 43 time steps.
+DEAL:0:2d::
+DEAL:0:3d::   Number of global active cells: 512
+DEAL:0:3d::   Number of degrees of freedom: 35937
+DEAL:0:3d::   Time step size: 0.10000, finest cell: 3.7500
+DEAL:0:3d::
+DEAL:0:3d::   Time:-10.0   , solution norm: 0.78144
+DEAL:0:3d::   Time:-8.90   , solution norm: 0.57949
+DEAL:0:3d::
+DEAL:0:3d::   Performed 11 time steps.
+DEAL:0:3d::

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.