]> https://gitweb.dealii.org/ - dealii.git/commitdiff
VectorTools::interpolate should not call compress on parallel::distributed::Vector...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 23 Nov 2012 10:04:01 +0000 (10:04 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 23 Nov 2012 10:04:01 +0000 (10:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@27672 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/parallel_vector.templates.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
tests/matrix_free/step-37.cc [new file with mode: 0644]
tests/matrix_free/step-37/cmp/generic [new file with mode: 0644]
tests/mpi/step-48.cc [new file with mode: 0644]
tests/mpi/step-48/ncpu_1/cmp/generic [new file with mode: 0644]
tests/mpi/step-48/ncpu_10/cmp/generic [new file with mode: 0644]
tests/mpi/step-48/ncpu_4/cmp/generic [new file with mode: 0644]

index cadadfcaca55373ce70c7e3b26bb683315d94dff..d841e0a1f5813d4a7f3f4f10fe7ae3e547185ca0 100644 (file)
@@ -190,6 +190,7 @@ namespace parallel
     {
       AssertDimension (local_range().first, c.local_range().first);
       AssertDimension (local_range().second, c.local_range().second);
+      AssertDimension (vector_view.size(), c.vector_view.size());
       vector_view = c.vector_view;
       if (call_update_ghost_values == true)
         update_ghost_values();
index 8546036d534f33e086df90516878eba252fac351..6936970500cfaea745d1c26a8ebb1a41c674742c 100644 (file)
@@ -66,6 +66,30 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace VectorTools
 {
+  namespace
+  {
+    template <typename VEC>
+    void perform_compress_insert (VEC &vec)
+    {
+      vec.compress(::dealii::VectorOperation::insert);
+    }
+
+    template <typename Number>
+    void perform_compress_insert (::dealii::parallel::distributed::Vector<Number> &vec)
+    {
+                                // should not do compress when inserting
+                                // elements as the ghosts are fixed and some
+                                // will not be set at all. Instead, zero the
+                                // ghost data afterwards to get consistent
+                                // data. Otherwise, at least with the layout
+                                // in r27671 the tests/mpi/step-48/ncpu_10
+                                // will signal incorrect results because of
+                                // incorrect interpolation.
+      vec.zero_out_ghosts();
+    }
+
+  }
+
 
   template <class VECTOR, class DH>
   void interpolate (const Mapping<DH::dimension,DH::space_dimension>    &mapping,
@@ -259,7 +283,7 @@ namespace VectorTools
                   = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
             }
         }
-    vec.compress(::dealii::VectorOperation::insert);
+    perform_compress_insert(vec);
   }
 
 
diff --git a/tests/matrix_free/step-37.cc b/tests/matrix_free/step-37.cc
new file mode 100644 (file)
index 0000000..005c227
--- /dev/null
@@ -0,0 +1,660 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------------------------
+
+// test for correctness of step-37 (without output and only small sizes)
+
+
+#include "../tests.h"
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_dof_handler.h>
+#include <deal.II/multigrid/mg_dof_accessor.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <fstream>
+#include <sstream>
+
+
+namespace Step37
+{
+  using namespace dealii;
+
+  const unsigned int degree_finite_element = 2;
+
+
+
+  template <int dim>
+  class Coefficient : public Function<dim>
+  {
+  public:
+    Coefficient ()  : Function<dim>() {}
+
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    template <typename number>
+    number value (const Point<dim,number> &p,
+                  const unsigned int component = 0) const;
+
+    virtual void value_list (const std::vector<Point<dim> > &points,
+                             std::vector<double>            &values,
+                             const unsigned int              component = 0) const;
+  };
+
+
+
+  template <int dim>
+  template <typename number>
+  number Coefficient<dim>::value (const Point<dim,number> &p,
+                                  const unsigned int /*component*/) const
+  {
+    return 1. / (0.05 + 2.*p.square());
+  }
+
+
+
+  template <int dim>
+  double Coefficient<dim>::value (const Point<dim> &p,
+                                  const unsigned int component) const
+  {
+    return value<double>(p,component);
+  }
+
+
+
+  template <int dim>
+  void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+                                     std::vector<double>            &values,
+                                     const unsigned int              component) const
+  {
+    Assert (values.size() == points.size(),
+            ExcDimensionMismatch (values.size(), points.size()));
+    Assert (component == 0,
+            ExcIndexRange (component, 0, 1));
+
+    const unsigned int n_points = points.size();
+    for (unsigned int i=0; i<n_points; ++i)
+      values[i] = value<double>(points[i],component);
+  }
+
+
+
+
+
+  template <int dim, int fe_degree, typename number>
+  class LaplaceOperator : public Subscriptor
+  {
+  public:
+    LaplaceOperator ();
+
+    void clear();
+
+    void reinit (const MGDoFHandler<dim> &dof_handler,
+                 const ConstraintMatrix &constraints,
+                 const unsigned int       level = numbers::invalid_unsigned_int);
+
+    unsigned int m () const;
+    unsigned int n () const;
+
+    void vmult (Vector<double> &dst,
+                const Vector<double> &src) const;
+    void Tvmult (Vector<double> &dst,
+                 const Vector<double> &src) const;
+    void vmult_add (Vector<double> &dst,
+                    const Vector<double> &src) const;
+    void Tvmult_add (Vector<double> &dst,
+                     const Vector<double> &src) const;
+
+    number el (const unsigned int row,
+               const unsigned int col) const;
+    void set_diagonal (const Vector<number> &diagonal);
+
+    std::size_t memory_consumption () const;
+
+  private:
+    void local_apply (const MatrixFree<dim,number>    &data,
+                      Vector<double>                      &dst,
+                      const Vector<double>                &src,
+                      const std::pair<unsigned int,unsigned int> &cell_range) const;
+
+    void evaluate_coefficient(const Coefficient<dim> &function);
+
+    MatrixFree<dim,number>      data;
+    AlignedVector<VectorizedArray<number> > coefficient;
+
+    Vector<number>  diagonal_values;
+    bool            diagonal_is_available;
+  };
+
+
+
+  template <int dim, int fe_degree, typename number>
+  LaplaceOperator<dim,fe_degree,number>::LaplaceOperator ()
+    :
+    Subscriptor()
+  {}
+
+
+
+  template <int dim, int fe_degree, typename number>
+  unsigned int
+  LaplaceOperator<dim,fe_degree,number>::m () const
+  {
+    return data.get_vector_partitioner()->size();
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  unsigned int
+  LaplaceOperator<dim,fe_degree,number>::n () const
+  {
+    return data.get_vector_partitioner()->size();
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::clear ()
+  {
+    data.clear();
+    diagonal_is_available = false;
+    diagonal_values.reinit(0);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::reinit (const MGDoFHandler<dim> &dof_handler,
+                                                 const ConstraintMatrix &constraints,
+                                                 const unsigned int      level)
+  {
+    typename MatrixFree<dim,number>::AdditionalData additional_data;
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim,number>::AdditionalData::partition_color;
+    additional_data.level_mg_handler = level;
+    additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
+                                            update_quadrature_points);
+    data.reinit (dof_handler, constraints, QGauss<1>(fe_degree+1),
+                 additional_data);
+    evaluate_coefficient(Coefficient<dim>());
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::
+  evaluate_coefficient (const Coefficient<dim> &coefficient_function)
+  {
+    const unsigned int n_cells = data.get_size_info().n_macro_cells;
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data);
+    coefficient.resize (n_cells * phi.n_q_points);
+    for (unsigned int cell=0; cell<n_cells; ++cell)
+      {
+        phi.reinit (cell);
+        for (unsigned int q=0; q<phi.n_q_points; ++q)
+          coefficient[cell*phi.n_q_points+q] =
+            coefficient_function.value(phi.quadrature_point(q));
+      }
+  }
+
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::
+  local_apply (const MatrixFree<dim,number>         &data,
+               Vector<double>                       &dst,
+               const Vector<double>                 &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,number> phi (data);
+    AssertDimension (coefficient.size(),
+                     data.get_size_info().n_macro_cells * phi.n_q_points);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        phi.reinit (cell);
+        phi.read_dof_values(src);
+        phi.evaluate (false,true,false);
+        for (unsigned int q=0; q<phi.n_q_points; ++q)
+          phi.submit_gradient (coefficient[cell*phi.n_q_points+q] *
+                               phi.get_gradient(q), q);
+        phi.integrate (false,true);
+        phi.distribute_local_to_global (dst);
+      }
+  }
+
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::vmult (Vector<double>       &dst,
+                                                const Vector<double> &src) const
+  {
+    dst = 0;
+    vmult_add (dst, src);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::Tvmult (Vector<double>       &dst,
+                                                 const Vector<double> &src) const
+  {
+    dst = 0;
+    vmult_add (dst,src);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::Tvmult_add (Vector<double>       &dst,
+                                                     const Vector<double> &src) const
+  {
+    vmult_add (dst,src);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::vmult_add (Vector<double>       &dst,
+                                                    const Vector<double> &src) const
+  {
+    data.cell_loop (&LaplaceOperator::local_apply, this, dst, src);
+
+    const std::vector<unsigned int> &
+    constrained_dofs = data.get_constrained_dofs();
+    for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+      dst(constrained_dofs[i]) += src(constrained_dofs[i]);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  number
+  LaplaceOperator<dim,fe_degree,number>::el (const unsigned int row,
+                                             const unsigned int col) const
+  {
+    Assert (row == col, ExcNotImplemented());
+    Assert (diagonal_is_available == true, ExcNotInitialized());
+    return diagonal_values(row);
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  void
+  LaplaceOperator<dim,fe_degree,number>::set_diagonal(const Vector<number> &diagonal)
+  {
+    AssertDimension (m(), diagonal.size());
+
+    diagonal_values = diagonal;
+
+    const std::vector<unsigned int> &
+    constrained_dofs = data.get_constrained_dofs();
+    for (unsigned int i=0; i<constrained_dofs.size(); ++i)
+      diagonal_values(constrained_dofs[i]) = 1.0;
+
+    diagonal_is_available = true;
+  }
+
+
+
+  template <int dim, int fe_degree, typename number>
+  std::size_t
+  LaplaceOperator<dim,fe_degree,number>::memory_consumption () const
+  {
+    return (data.memory_consumption () +
+            coefficient.memory_consumption() +
+            diagonal_values.memory_consumption() +
+            MemoryConsumption::memory_consumption(diagonal_is_available));
+  }
+
+
+
+
+  template <int dim>
+  class LaplaceProblem
+  {
+  public:
+    LaplaceProblem ();
+    void run ();
+
+  private:
+    void setup_system ();
+    void assemble_system ();
+    void assemble_multigrid ();
+    void solve ();
+    void output_results (const unsigned int cycle) const;
+
+    typedef LaplaceOperator<dim,degree_finite_element,double> SystemMatrixType;
+    typedef LaplaceOperator<dim,degree_finite_element,float>  LevelMatrixType;
+
+    Triangulation<dim>               triangulation;
+    FE_Q<dim>                        fe;
+    MGDoFHandler<dim>                mg_dof_handler;
+    ConstraintMatrix                 constraints;
+
+    SystemMatrixType                 system_matrix;
+    MGLevelObject<LevelMatrixType>   mg_matrices;
+    FullMatrix<float>                coarse_matrix;
+    MGLevelObject<ConstraintMatrix>  mg_constraints;
+
+    Vector<double>                   solution;
+    Vector<double>                   system_rhs;
+  };
+
+
+
+  template <int dim>
+  LaplaceProblem<dim>::LaplaceProblem ()
+    :
+    fe (degree_finite_element),
+    mg_dof_handler (triangulation)
+  {}
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::setup_system ()
+  {
+    system_matrix.clear();
+    mg_matrices.clear();
+    mg_constraints.clear();
+
+    mg_dof_handler.distribute_dofs (fe);
+
+    deallog << "Number of degrees of freedom: "
+            << mg_dof_handler.n_dofs()
+            << std::endl;
+
+    constraints.clear();
+    VectorTools::interpolate_boundary_values (mg_dof_handler,
+                                              0,
+                                              ZeroFunction<dim>(),
+                                              constraints);
+    constraints.close();
+
+    system_matrix.reinit (mg_dof_handler, constraints);
+    solution.reinit (mg_dof_handler.n_dofs());
+    system_rhs.reinit (mg_dof_handler.n_dofs());
+
+    const unsigned int nlevels = triangulation.n_levels();
+    mg_matrices.resize(0, nlevels-1);
+    mg_constraints.resize (0, nlevels-1);
+
+    typename FunctionMap<dim>::type dirichlet_boundary;
+    ZeroFunction<dim>               homogeneous_dirichlet_bc (1);
+    dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
+    std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
+    MGTools::make_boundary_list (mg_dof_handler,
+                                 dirichlet_boundary,
+                                 boundary_indices);
+    for (unsigned int level=0; level<nlevels; ++level)
+      {
+        std::set<unsigned int>::iterator bc_it = boundary_indices[level].begin();
+        for ( ; bc_it != boundary_indices[level].end(); ++bc_it)
+          mg_constraints[level].add_line(*bc_it);
+
+        mg_constraints[level].close();
+        mg_matrices[level].reinit(mg_dof_handler,
+                                  mg_constraints[level],
+                                  level);
+      }
+    coarse_matrix.reinit (mg_dof_handler.n_dofs(0),
+                          mg_dof_handler.n_dofs(0));
+  }
+
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_system ()
+  {
+    QGauss<dim>  quadrature_formula(fe.degree+1);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_values   | update_JxW_values);
+
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    const Coefficient<dim> coefficient;
+    std::vector<double>    coefficient_values (n_q_points);
+
+    typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(),
+                                                   endc = mg_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell->get_dof_indices (local_dof_indices);
+        fe_values.reinit (cell);
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            double rhs_val = 0;
+            for (unsigned int q=0; q<n_q_points; ++q)
+              rhs_val += (fe_values.shape_value(i,q) * 1.0 *
+                          fe_values.JxW(q));
+            system_rhs(local_dof_indices[i]) += rhs_val;
+          }
+      }
+    constraints.condense(system_rhs);
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::assemble_multigrid ()
+  {
+    coarse_matrix = 0;
+    QGauss<dim>  quadrature_formula(fe.degree+1);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                             update_gradients  | update_inverse_jacobians |
+                             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();
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    const Coefficient<dim>    coefficient;
+    std::vector<double>       coefficient_values (n_q_points);
+    FullMatrix<double>        local_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>            local_diagonal (dofs_per_cell);
+
+    const unsigned int n_levels = triangulation.n_levels();
+    std::vector<Vector<float> > diagonals (n_levels);
+    for (unsigned int level=0; level<n_levels; ++level)
+      diagonals[level].reinit (mg_dof_handler.n_dofs(level));
+
+    std::vector<unsigned int> cell_no(triangulation.n_levels());
+    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+                                              endc = mg_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        const unsigned int level = cell->level();
+        cell->get_mg_dof_indices (local_dof_indices);
+        fe_values.reinit (cell);
+        coefficient.value_list (fe_values.get_quadrature_points(),
+                                coefficient_values);
+
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            double local_diag = 0;
+            for (unsigned int q=0; q<n_q_points; ++q)
+              local_diag += ((fe_values.shape_grad(i,q) *
+                              fe_values.shape_grad(i,q)) *
+                             coefficient_values[q] * fe_values.JxW(q));
+            local_diagonal(i) = local_diag;
+          }
+        mg_constraints[level].distribute_local_to_global(local_diagonal,
+                                                         local_dof_indices,
+                                                         diagonals[level]);
+
+        if (level == 0)
+          {
+            local_matrix = 0;
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                {
+                  double add_value = 0;
+                  for (unsigned int q=0; q<n_q_points; ++q)
+                    add_value += (fe_values.shape_grad(i,q) *
+                                  fe_values.shape_grad(j,q) *
+                                  coefficient_values[q] *
+                                  fe_values.JxW(q));
+                  local_matrix(i,j) = add_value;
+                }
+            mg_constraints[0].distribute_local_to_global (local_matrix,
+                                                          local_dof_indices,
+                                                          coarse_matrix);
+          }
+      }
+
+    for (unsigned int level=0; level<n_levels; ++level)
+      mg_matrices[level].set_diagonal (diagonals[level]);
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::solve ()
+  {
+    GrowingVectorMemory<>   vector_memory;
+
+    MGTransferPrebuilt<Vector<double> > mg_transfer;
+    mg_transfer.build_matrices(mg_dof_handler);
+
+    MGCoarseGridHouseholder<float, Vector<double> > mg_coarse;
+    mg_coarse.initialize(coarse_matrix);
+
+    typedef PreconditionChebyshev<LevelMatrixType,Vector<double> > SMOOTHER;
+    MGSmootherPrecondition<LevelMatrixType, SMOOTHER, Vector<double> >
+    mg_smoother(vector_memory);
+
+    typename SMOOTHER::AdditionalData smoother_data;
+    smoother_data.smoothing_range = 10.;
+    smoother_data.degree = 6;
+    smoother_data.eig_cg_n_iterations = 10;
+    mg_smoother.initialize(mg_matrices, smoother_data);
+
+    MGMatrix<LevelMatrixType, Vector<double> >
+    mg_matrix(&mg_matrices);
+
+    Multigrid<Vector<double> > mg(mg_dof_handler,
+                                  mg_matrix,
+                                  mg_coarse,
+                                  mg_transfer,
+                                  mg_smoother,
+                                  mg_smoother);
+    PreconditionMG<dim, Vector<double>,
+                   MGTransferPrebuilt<Vector<double> > >
+                   preconditioner(mg_dof_handler, mg, mg_transfer);
+
+    const std::size_t multigrid_memory
+      = (mg_matrices.memory_consumption() +
+         mg_transfer.memory_consumption() +
+         coarse_matrix.memory_consumption());
+
+    SolverControl           solver_control (1000, 1e-12*system_rhs.l2_norm());
+    SolverCG<>              cg (solver_control);
+
+    cg.solve (system_matrix, solution, system_rhs,
+              preconditioner);
+  }
+
+
+
+  template <int dim>
+  void LaplaceProblem<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<3; ++cycle)
+      {
+        deallog << "Cycle " << cycle << std::endl;
+
+        if (cycle == 0)
+          {
+            GridGenerator::hyper_cube (triangulation, 0., 1.);
+            triangulation.refine_global (3-dim);
+          }
+        triangulation.refine_global (1);
+        setup_system ();
+        assemble_system ();
+        assemble_multigrid ();
+        solve ();
+        deallog << std::endl;
+      };
+  }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("step-37/output");
+  deallog.attach(logfile);
+  deallog << std::setprecision (3);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console(0);
+
+  {
+    deallog.push("2d");
+    Step37::LaplaceProblem<2> laplace_problem;
+    laplace_problem.run();
+    deallog.pop();
+  }
+  {
+    deallog.push("3d");
+    Step37::LaplaceProblem<3> laplace_problem;
+    laplace_problem.run();
+    deallog.pop();
+  }
+
+
+  return 0;
+}
+
diff --git a/tests/matrix_free/step-37/cmp/generic b/tests/matrix_free/step-37/cmp/generic
new file mode 100644 (file)
index 0000000..4ac0a49
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL:2d::Cycle 0
+DEAL:2d::Number of degrees of freedom: 81
+DEAL:2d:cg::Starting value 0.132
+DEAL:2d:cg::Convergence step 5 value 0
+DEAL:2d::
+DEAL:2d::Cycle 1
+DEAL:2d::Number of degrees of freedom: 289
+DEAL:2d:cg::Starting value 0.0677
+DEAL:2d:cg::Convergence step 5 value 0
+DEAL:2d::
+DEAL:2d::Cycle 2
+DEAL:2d::Number of degrees of freedom: 1089
+DEAL:2d:cg::Starting value 0.0343
+DEAL:2d:cg::Convergence step 6 value 0
+DEAL:2d::
+DEAL:3d::Cycle 0
+DEAL:3d::Number of degrees of freedom: 125
+DEAL:3d:cg::Starting value 0.125
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d::
+DEAL:3d::Cycle 1
+DEAL:3d::Number of degrees of freedom: 729
+DEAL:3d:cg::Starting value 0.0479
+DEAL:3d:cg::Convergence step 5 value 0
+DEAL:3d::
+DEAL:3d::Cycle 2
+DEAL:3d::Number of degrees of freedom: 4913
+DEAL:3d:cg::Starting value 0.0176
+DEAL:3d:cg::Convergence step 4 value 0
+DEAL:3d::
diff --git a/tests/mpi/step-48.cc b/tests/mpi/step-48.cc
new file mode 100644 (file)
index 0000000..d8092d4
--- /dev/null
@@ -0,0 +1,426 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------------------------
+
+// test for correctness of step-48 (without output and only small sizes)
+
+
+#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();
+    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 oldstyle_operation ();
+    void assemble_system ();
+    void output_results (const unsigned int timestep_number) const;
+
+    parallel::distributed::Triangulation<dim>   triangulation;
+    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 ()
+    :
+    triangulation (MPI_COMM_WORLD),
+    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);
+    {
+      typename Triangulation<dim>::active_cell_iterator
+      cell = triangulation.begin_active(),
+      end_cell = triangulation.end();
+      for ( ; cell != end_cell; ++cell)
+        if (cell->is_locally_owned())
+          if (cell->center().norm() < 11)
+            cell->set_refine_flag();
+      triangulation.execute_coarsening_and_refinement();
+
+      cell = triangulation.begin_active();
+      end_cell = triangulation.end();
+      for ( ; cell != end_cell; ++cell)
+        if (cell->is_locally_owned())
+          if (cell->center().norm() < 6)
+            cell->set_refine_flag();
+      triangulation.execute_coarsening_and_refinement();
+    }
+
+    deallog << "   Number of global active cells: "
+            << triangulation.n_global_active_cells()
+            << 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) const
+  {
+    parallel::distributed::Vector<double> locally_relevant_solution;
+    locally_relevant_solution.reinit (dof_handler.locally_owned_dofs(),
+                                      locally_relevant_dofs,
+                                      MPI_COMM_WORLD);
+    locally_relevant_solution.copy_from (solution);
+    locally_relevant_solution.update_ghost_values ();
+    constraints.distribute (locally_relevant_solution);
+
+    Vector<float> norm_per_cell (triangulation.n_active_cells());
+    VectorTools::integrate_difference (dof_handler,
+                                       locally_relevant_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::System::get_this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile(output_file_for_mpi("step-48").c_str());
+      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/mpi/step-48/ncpu_1/cmp/generic b/tests/mpi/step-48/ncpu_1/cmp/generic
new file mode 100644 (file)
index 0000000..f0b2c39
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:0:2d::   Number of global active cells: 964
+DEAL:0:2d::   Number of degrees of freedom: 16033
+DEAL:0:2d::   Time step size: 0.01176, finest cell: 0.4688
+DEAL:0:2d::
+DEAL:0:2d::   Time:-10.0   , solution norm: 9.5599 
+DEAL:0:2d::   Time:-8.99   , solution norm: 21.965 
+DEAL:0:2d::
+DEAL:0:2d::   Performed 86 time steps.
+DEAL:0:2d::
+DEAL:0:3d::   Number of global active cells: 176
+DEAL:0:3d::   Number of degrees of freedom: 14001
+DEAL:0:3d::   Time step size: 0.047619, finest cell: 1.8750
+DEAL:0:3d::
+DEAL:0:3d::   Time:-10.0   , solution norm: 29.559 
+DEAL:0:3d::   Time:-8.95   , solution norm: 90.729 
+DEAL:0:3d::
+DEAL:0:3d::   Performed 22 time steps.
+DEAL:0:3d::
diff --git a/tests/mpi/step-48/ncpu_10/cmp/generic b/tests/mpi/step-48/ncpu_10/cmp/generic
new file mode 100644 (file)
index 0000000..f0b2c39
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:0:2d::   Number of global active cells: 964
+DEAL:0:2d::   Number of degrees of freedom: 16033
+DEAL:0:2d::   Time step size: 0.01176, finest cell: 0.4688
+DEAL:0:2d::
+DEAL:0:2d::   Time:-10.0   , solution norm: 9.5599 
+DEAL:0:2d::   Time:-8.99   , solution norm: 21.965 
+DEAL:0:2d::
+DEAL:0:2d::   Performed 86 time steps.
+DEAL:0:2d::
+DEAL:0:3d::   Number of global active cells: 176
+DEAL:0:3d::   Number of degrees of freedom: 14001
+DEAL:0:3d::   Time step size: 0.047619, finest cell: 1.8750
+DEAL:0:3d::
+DEAL:0:3d::   Time:-10.0   , solution norm: 29.559 
+DEAL:0:3d::   Time:-8.95   , solution norm: 90.729 
+DEAL:0:3d::
+DEAL:0:3d::   Performed 22 time steps.
+DEAL:0:3d::
diff --git a/tests/mpi/step-48/ncpu_4/cmp/generic b/tests/mpi/step-48/ncpu_4/cmp/generic
new file mode 100644 (file)
index 0000000..f0b2c39
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:0:2d::   Number of global active cells: 964
+DEAL:0:2d::   Number of degrees of freedom: 16033
+DEAL:0:2d::   Time step size: 0.01176, finest cell: 0.4688
+DEAL:0:2d::
+DEAL:0:2d::   Time:-10.0   , solution norm: 9.5599 
+DEAL:0:2d::   Time:-8.99   , solution norm: 21.965 
+DEAL:0:2d::
+DEAL:0:2d::   Performed 86 time steps.
+DEAL:0:2d::
+DEAL:0:3d::   Number of global active cells: 176
+DEAL:0:3d::   Number of degrees of freedom: 14001
+DEAL:0:3d::   Time step size: 0.047619, finest cell: 1.8750
+DEAL:0:3d::
+DEAL:0:3d::   Time:-10.0   , solution norm: 29.559 
+DEAL:0:3d::   Time:-8.95   , solution norm: 90.729 
+DEAL:0:3d::
+DEAL:0:3d::   Performed 22 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.