]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add periodicity to MGTransfer matrix free and prebuilt 4527/head
authortcclevenger <tcleven@clemson.edu>
Mon, 19 Jun 2017 07:56:33 +0000 (09:56 +0200)
committertcclevenger <tcleven@clemson.edu>
Sat, 8 Jul 2017 12:20:34 +0000 (08:20 -0400)
doc/news/changes/minor/20170708ConradClevenger [new file with mode: 0644]
include/deal.II/multigrid/mg_constrained_dofs.h
source/multigrid/mg_transfer_internal.cc
source/multigrid/mg_transfer_prebuilt.cc
tests/multigrid/transfer_matrix_free_11.cc [new file with mode: 0644]
tests/multigrid/transfer_matrix_free_11.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20170708ConradClevenger b/doc/news/changes/minor/20170708ConradClevenger
new file mode 100644 (file)
index 0000000..6494e6a
--- /dev/null
@@ -0,0 +1,5 @@
+New: MGConstrainedDoFs object now stores information for periodicity 
+constraints which is used during the multigrid transfer by 
+MGTransferMatrixFree and MGTransferPrebuilt.
+<br>
+(Conrad Clevenger, Martin Kronbichler, 2017/07/08)
index 17c1c84b95c5fd22f577a0ced47e738f1afcd627..ea8b2e7cedff67a058222db0ee1a3496cfdd2255 100644 (file)
@@ -50,6 +50,15 @@ public:
    * edges of a refinement level are treated corrected but leaves degrees of
    * freedom at the boundary of the domain untouched assuming that no
    * Dirichlet boundary conditions for them exist.
+   *
+   * Furthermore, this call sets up a ConstraintMatrix on each level that
+   * contains possible periodicity constraints in case those have been added to
+   * the underlying triangulation. The constraint matrix can be queried by
+   * get_level_constraint_matrix(level). Note that the current implementation of
+   * periodicity constraints in this class does not support rotation matrices in the
+   * periodicity definition, i.e., the respective argument in the
+   * GridTools::collect_periodic_faces() may not be different from the identity
+   * matrix.
    */
   template <int dim, int spacedim>
   void initialize (const DoFHandler<dim,spacedim> &dof);
@@ -124,6 +133,13 @@ public:
    */
   bool have_boundary_indices () const;
 
+  /**
+   * Return the level constraint matrix for a given level, containing
+   * periodicity constraints (if enabled on the triangulation).
+   */
+  const ConstraintMatrix &
+  get_level_constraint_matrix (const unsigned int level) const;
+
 private:
 
   /**
@@ -136,6 +152,13 @@ private:
    * between the level and cells on a coarser level.
    */
   std::vector<IndexSet> refinement_edge_indices;
+
+  /**
+   * Constraint matrices containing information regarding potential
+   * periodic boundary conditions for each level .
+   */
+  std::vector<ConstraintMatrix> level_constraints;
+
 };
 
 
@@ -146,13 +169,70 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim,spacedim> &dof)
 {
   boundary_indices.clear();
   refinement_edge_indices.clear();
+  level_constraints.clear();
 
   const unsigned int nlevels = dof.get_triangulation().n_global_levels();
 
-  //At this point refinement_edge_indices is empty.
+  //At this point level_constraint and refinement_edge_indices are empty.
+  level_constraints.resize(nlevels);
   refinement_edge_indices.resize(nlevels);
   for (unsigned int l=0; l<nlevels; ++l)
-    refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+    {
+      IndexSet relevant_dofs;
+      DoFTools::extract_locally_relevant_level_dofs(dof,l,relevant_dofs);
+      level_constraints[l].reinit(relevant_dofs);
+
+      // Loop through relevant cells and faces finding those which are periodic neighbors.
+      typename DoFHandler<dim,spacedim>::cell_iterator
+      cell = dof.begin(l),
+      endc = dof.end(l);
+      for (; cell!=endc; ++cell)
+        if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
+          {
+            for (unsigned int f = 0;
+                 f < GeometryInfo<dim>::faces_per_cell;
+                 ++f)
+              if (cell->has_periodic_neighbor(f))
+                {
+                  if (cell->is_locally_owned_on_level())
+                    {
+                      Assert(cell->periodic_neighbor(f)->level_subdomain_id() != numbers::artificial_subdomain_id,
+                             ExcMessage ("Periodic neighbor of a locally owned cell must either be owned or ghost."));
+                    }
+                  // Cell is a level-ghost and its neighbor is a level-artificial cell
+                  // nothing to do here
+                  else if (cell->periodic_neighbor(f)->level_subdomain_id() == numbers::artificial_subdomain_id)
+                    {
+                      Assert (cell->is_locally_owned_on_level() == false, ExcInternalError());
+                      continue;
+                    }
+
+                  const unsigned int dofs_per_face = cell->face(f)->get_fe(0).dofs_per_face;
+                  std::vector<types::global_dof_index> dofs_1(dofs_per_face);
+                  std::vector<types::global_dof_index> dofs_2(dofs_per_face);
+
+                  cell->periodic_neighbor(f)->face(cell->periodic_neighbor_face_no(f))->get_mg_dof_indices(l, dofs_1, 0);
+                  cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
+                  // Store periodicity information in the level constraint matrix
+                  // Skip DoFs for which we've previously entered periodicity constraints already;
+                  // this can happen, for example, for a vertex dof at a periodic boundary that we
+                  // visit from more than one cell
+                  for (unsigned int i=0; i<dofs_per_face; ++i)
+                    if (level_constraints[l].can_store_line(dofs_2[i]) &&
+                        level_constraints[l].can_store_line(dofs_1[i]) &&
+                        !level_constraints[l].is_constrained(dofs_2[i]) &&
+                        !level_constraints[l].is_constrained(dofs_1[i]))
+                      {
+                        level_constraints[l].add_line(dofs_2[i]);
+                        level_constraints[l].add_entry(dofs_2[i], dofs_1[i], 1.);
+                      }
+                }
+          }
+      level_constraints[l].close();
+
+      // Initialize with empty IndexSet of correct size
+      refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+    }
 
   MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices);
 }
@@ -264,6 +344,18 @@ MGConstrainedDoFs::have_boundary_indices () const
 }
 
 
+
+
+inline
+const ConstraintMatrix &
+MGConstrainedDoFs::get_level_constraint_matrix (const unsigned int level) const
+{
+  AssertIndexRange(level, level_constraints.size());
+  return level_constraints[level];
+}
+
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 56121d9615e302f24e5caa5823db3187d7a85da1..9d157ddf4961d738dbc6f1cf5d7edf59c3f3e9eb 100644 (file)
@@ -443,6 +443,27 @@ namespace internal
               fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]);
     }
 
+    namespace
+    {
+      /**
+       * Helper function for setup_transfer. Checks for identity constrained dofs
+       * and replace with the indices of the dofs to which they are constrained
+       */
+      void replace (const MGConstrainedDoFs              *mg_constrained_dofs,
+                    const unsigned int                   level,
+                    std::vector<types::global_dof_index> &dof_indices)
+      {
+        if (mg_constrained_dofs != nullptr &&
+            mg_constrained_dofs->get_level_constraint_matrix(level).n_constraints() > 0)
+          for (auto &ind : dof_indices)
+            if (mg_constrained_dofs->get_level_constraint_matrix(level).is_identity_constrained(ind))
+              {
+                Assert(mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->size() == 1,
+                       ExcInternalError());
+                ind = mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->front().first;
+              }
+      }
+    }
 
     // Sets up most of the internal data structures of the MGTransferMatrixFree
     // class
@@ -564,6 +585,8 @@ namespace internal
                     continue;
                   cell->child(c)->get_mg_dof_indices(local_dof_indices);
 
+                  replace(mg_constrained_dofs, level, local_dof_indices);
+
                   const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
                   for (unsigned int i=0; i<local_dof_indices.size(); ++i)
                     if (!owned_level_dofs.is_element(local_dof_indices[i]))
@@ -620,6 +643,8 @@ namespace internal
                 {
                   cell->get_mg_dof_indices(local_dof_indices);
 
+                  replace(mg_constrained_dofs, level-1, local_dof_indices);
+
                   const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
                   for (unsigned int i=0; i<local_dof_indices.size(); ++i)
                     if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
index 69e7f14b94d9891c37e5a2b355249493f0bb58c1..6c0a4db28ebd84a1ba2dffcda8bf58314ef65696 100644 (file)
@@ -123,6 +123,27 @@ void MGTransferPrebuilt<VectorType>::restrict_and_add (const unsigned int from_l
 }
 
 
+namespace
+{
+  /**
+   * Helper function for build_matrices. Checks for identity constrained dofs
+   * and replace with the indices of the dofs to which they are constrained
+   */
+  void replace (const MGConstrainedDoFs              *mg_constrained_dofs,
+                const unsigned int                   level,
+                std::vector<types::global_dof_index> &dof_indices)
+  {
+    if (mg_constrained_dofs != nullptr &&
+        mg_constrained_dofs->get_level_constraint_matrix(level).n_constraints() > 0)
+      for (auto &ind : dof_indices)
+        if (mg_constrained_dofs->get_level_constraint_matrix(level).is_identity_constrained(ind))
+          {
+            Assert(mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->size() == 1,
+                   ExcInternalError());
+            ind = mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->front().first;
+          }
+  }
+}
 
 template <typename VectorType>
 template <int dim, int spacedim>
@@ -195,6 +216,8 @@ void MGTransferPrebuilt<VectorType>::build_matrices
           {
             cell->get_mg_dof_indices (dof_indices_parent);
 
+            replace(this->mg_constrained_dofs, level, dof_indices_parent);
+
             Assert(cell->n_children()==GeometryInfo<dim>::max_children_per_cell,
                    ExcNotImplemented());
             for (unsigned int child=0; child<cell->n_children(); ++child)
@@ -208,6 +231,8 @@ void MGTransferPrebuilt<VectorType>::build_matrices
 
                 cell->child(child)->get_mg_dof_indices (dof_indices_child);
 
+                replace(this->mg_constrained_dofs, level+1, dof_indices_child);
+
                 // now tag the entries in the
                 // matrix which will be used
                 // for this pair of parent/child
@@ -241,6 +266,8 @@ void MGTransferPrebuilt<VectorType>::build_matrices
           {
             cell->get_mg_dof_indices (dof_indices_parent);
 
+            replace(this->mg_constrained_dofs, level, dof_indices_parent);
+
             Assert(cell->n_children()==GeometryInfo<dim>::max_children_per_cell,
                    ExcNotImplemented());
             for (unsigned int child=0; child<cell->n_children(); ++child)
@@ -259,6 +286,8 @@ void MGTransferPrebuilt<VectorType>::build_matrices
 
                 cell->child(child)->get_mg_dof_indices (dof_indices_child);
 
+                replace(this->mg_constrained_dofs, level+1, dof_indices_child);
+
                 // now set the entries in the matrix
                 for (unsigned int i=0; i<dofs_per_cell; ++i)
                   prolongation_matrices[level]->set (dof_indices_child[i],
diff --git a/tests/multigrid/transfer_matrix_free_11.cc b/tests/multigrid/transfer_matrix_free_11.cc
new file mode 100644 (file)
index 0000000..9791ec7
--- /dev/null
@@ -0,0 +1,135 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check MGTransferMatrixFree and MGTransferPrebuilt for periodic boundary
+// conditions
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+
+template <int dim, typename Number>
+void check(const unsigned int fe_degree)
+{
+  deallog.threshold_double(std::max(5e2*(double)std::numeric_limits<Number>::epsilon(),
+                                    1e-11));
+  FE_Q<dim> fe(fe_degree);
+  deallog << "FE: " << fe.get_name() << std::endl;
+
+  parallel::distributed::Triangulation<dim>
+  tr(MPI_COMM_WORLD,
+     Triangulation<dim>::limit_level_difference_at_vertices,
+     parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+  GridGenerator::hyper_cube(tr, -1, 1, true);
+  std::vector<GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator> >
+  periodicity_vector;
+  GridTools::collect_periodic_faces(tr, 0, 1, 0,
+                                    periodicity_vector);
+  GridTools::collect_periodic_faces(tr, 2, 3, 1,
+                                    periodicity_vector);
+  tr.add_periodicity(periodicity_vector);
+  tr.refine_global(2);
+  deallog << "no. cells: " << tr.n_global_active_cells() << std::endl;
+
+  DoFHandler<dim> mgdof(tr);
+  mgdof.distribute_dofs(fe);
+  mgdof.distribute_mg_dofs(fe);
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(mgdof);
+
+  // build reference
+  MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >
+  transfer_ref(mg_constrained_dofs);
+  transfer_ref.build_matrices(mgdof);
+  deallog << "Transfer matrices: " << std::endl;
+  transfer_ref.print_matrices(deallog.get_file_stream());
+  deallog << std::endl;
+
+  // build matrix-free transfer
+  MGTransferMatrixFree<dim, Number> transfer(mg_constrained_dofs);
+  transfer.build(mgdof);
+
+  // check prolongation for all levels using random vector
+  for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+    {
+      LinearAlgebra::distributed::Vector<Number> v1, v2;
+      LinearAlgebra::distributed::Vector<double> v1_cpy, v2_cpy, v3;
+      v1.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+      v2.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+      v3.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+      for (unsigned int i=0; i<v1.local_size(); ++i)
+        v1.local_element(i) = (double)Testing::rand()/RAND_MAX;
+      v1_cpy = v1;
+      transfer.prolongate(level, v2, v1);
+      transfer_ref.prolongate(level, v3, v1_cpy);
+      v2_cpy = v2;
+      v3 -= v2_cpy;
+      deallog << "Diff prolongate   l" << level << ": " << v3.l2_norm() << std::endl;
+      if (v3.l2_norm() > 1e-12)
+        {
+          // On level 0, we expect the matrix-based constraints to be wrong
+          // because it cannot capture the periodicity connections with a
+          // single cell
+          v2_cpy.print(deallog.get_file_stream());
+          v3.print(deallog.get_file_stream());
+        }
+    }
+
+  // check restriction for all levels using random vector
+  for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+    {
+      LinearAlgebra::distributed::Vector<Number> v1, v2;
+      LinearAlgebra::distributed::Vector<double> v1_cpy, v2_cpy, v3;
+      v1.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+      v2.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+      v3.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+      for (unsigned int i=0; i<v1.local_size(); ++i)
+        v1.local_element(i) = (double)Testing::rand()/RAND_MAX;
+      v1_cpy = v1;
+      transfer.restrict_and_add(level, v2, v1);
+      transfer_ref.restrict_and_add(level, v3, v1_cpy);
+      v2_cpy = v2;
+      v3 -= v2_cpy;
+      deallog << "Diff restrict     l" << level << ": " << v3.l2_norm() << std::endl;
+
+      v2 = 1.;
+      v3 = 1.;
+      transfer.restrict_and_add(level, v2, v1);
+      transfer_ref.restrict_and_add(level, v3, v1_cpy);
+      v2_cpy = v2;
+      v3 -= v2_cpy;
+      deallog << "Diff restrict add l" << level << ": " << v3.l2_norm() << std::endl;
+    }
+  deallog << std::endl;
+}
+
+
+int main(int argc, char **argv)
+{
+  // no threading in this test...
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  mpi_initlog();
+
+  check<2,double>(1);
+}
diff --git a/tests/multigrid/transfer_matrix_free_11.output b/tests/multigrid/transfer_matrix_free_11.output
new file mode 100644 (file)
index 0000000..49963ec
--- /dev/null
@@ -0,0 +1,64 @@
+
+DEAL::FE: FE_Q<2>(1)
+DEAL::no. cells: 16
+DEAL::Transfer matrices: 
+Level 0
+(3,3) 0.250000
+(5,3) 0.500000
+(7,3) 0.500000
+(8,3) 1.00000
+
+Level 1
+(3,3) 0.250000
+(3,5) 0.250000
+(3,7) 0.250000
+(3,8) 0.250000
+(5,3) 0.500000
+(5,7) 0.500000
+(7,3) 0.500000
+(7,5) 0.500000
+(8,3) 1.00000
+(10,3) 0.250000
+(10,5) 0.250000
+(10,7) 0.250000
+(10,8) 0.250000
+(12,5) 0.500000
+(12,8) 0.500000
+(13,3) 0.500000
+(13,5) 0.500000
+(14,5) 1.00000
+(16,3) 0.250000
+(16,5) 0.250000
+(16,7) 0.250000
+(16,8) 0.250000
+(17,3) 0.500000
+(17,7) 0.500000
+(19,7) 0.500000
+(19,8) 0.500000
+(20,7) 1.00000
+(21,3) 0.250000
+(21,5) 0.250000
+(21,7) 0.250000
+(21,8) 0.250000
+(22,5) 0.500000
+(22,8) 0.500000
+(23,7) 0.500000
+(23,8) 0.500000
+(24,8) 1.00000
+
+DEAL::
+DEAL::Diff prolongate   l1: 0.823013
+Process #0
+Local range: [0, 9), global size: 9
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 7.984e-01 0.000e+00 7.984e-01 0.000e+00 7.984e-01 7.984e-01 
+Process #0
+Local range: [0, 9), global size: 9
+Vector data:
+0.000e+00 0.000e+00 0.000e+00 -5.988e-01 0.000e+00 -3.992e-01 0.000e+00 -3.992e-01 0.000e+00 
+DEAL::Diff prolongate   l2: 0
+DEAL::Diff restrict     l1: 0.555735
+DEAL::Diff restrict add l1: 0.555735
+DEAL::Diff restrict     l2: 0
+DEAL::Diff restrict add l2: 0
+DEAL::

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.