From: tcclevenger Date: Mon, 19 Jun 2017 07:56:33 +0000 (+0200) Subject: add periodicity to MGTransfer matrix free and prebuilt X-Git-Tag: v9.0.0-rc1~1433^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4527%2Fhead;p=dealii.git add periodicity to MGTransfer matrix free and prebuilt --- diff --git a/doc/news/changes/minor/20170708ConradClevenger b/doc/news/changes/minor/20170708ConradClevenger new file mode 100644 index 0000000000..6494e6acd7 --- /dev/null +++ b/doc/news/changes/minor/20170708ConradClevenger @@ -0,0 +1,5 @@ +New: MGConstrainedDoFs object now stores information for periodicity +constraints which is used during the multigrid transfer by +MGTransferMatrixFree and MGTransferPrebuilt. +
+(Conrad Clevenger, Martin Kronbichler, 2017/07/08) diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 17c1c84b95..ea8b2e7ced 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -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 void initialize (const DoFHandler &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 refinement_edge_indices; + + /** + * Constraint matrices containing information regarding potential + * periodic boundary conditions for each level . + */ + std::vector level_constraints; + }; @@ -146,13 +169,70 @@ MGConstrainedDoFs::initialize(const DoFHandler &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::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::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 dofs_1(dofs_per_face); + std::vector 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 &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; iget_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::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 &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 template @@ -195,6 +216,8 @@ void MGTransferPrebuilt::build_matrices { cell->get_mg_dof_indices (dof_indices_parent); + replace(this->mg_constrained_dofs, level, dof_indices_parent); + Assert(cell->n_children()==GeometryInfo::max_children_per_cell, ExcNotImplemented()); for (unsigned int child=0; childn_children(); ++child) @@ -208,6 +231,8 @@ void MGTransferPrebuilt::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::build_matrices { cell->get_mg_dof_indices (dof_indices_parent); + replace(this->mg_constrained_dofs, level, dof_indices_parent); + Assert(cell->n_children()==GeometryInfo::max_children_per_cell, ExcNotImplemented()); for (unsigned int child=0; childn_children(); ++child) @@ -259,6 +286,8 @@ void MGTransferPrebuilt::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; iset (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 index 0000000000..9791ec7e09 --- /dev/null +++ b/tests/multigrid/transfer_matrix_free_11.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include + + +template +void check(const unsigned int fe_degree) +{ + deallog.threshold_double(std::max(5e2*(double)std::numeric_limits::epsilon(), + 1e-11)); + FE_Q fe(fe_degree); + deallog << "FE: " << fe.get_name() << std::endl; + + parallel::distributed::Triangulation + tr(MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + GridGenerator::hyper_cube(tr, -1, 1, true); + std::vector::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 mgdof(tr); + mgdof.distribute_dofs(fe); + mgdof.distribute_mg_dofs(fe); + + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(mgdof); + + // build reference + MGTransferPrebuilt > + 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 transfer(mg_constrained_dofs); + transfer.build(mgdof); + + // check prolongation for all levels using random vector + for (unsigned int level=1; level v1, v2; + LinearAlgebra::distributed::Vector 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 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 v1, v2; + LinearAlgebra::distributed::Vector 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(1); +} diff --git a/tests/multigrid/transfer_matrix_free_11.output b/tests/multigrid/transfer_matrix_free_11.output new file mode 100644 index 0000000000..49963ec916 --- /dev/null +++ b/tests/multigrid/transfer_matrix_free_11.output @@ -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::