From 0b2f314779f9ea21c867c37d020f6f0cfacbc994 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Mon, 6 Aug 2018 09:42:05 -0400 Subject: [PATCH] Add support for hanging nodes --- .../deal.II/matrix_free/cuda_fe_evaluation.h | 40 +- .../matrix_free/cuda_hanging_nodes_internal.h | 985 ++++++++++++++++++ .../matrix_free/cuda_matrix_free.templates.h | 45 +- source/matrix_free/CMakeLists.txt | 4 - 4 files changed, 1026 insertions(+), 48 deletions(-) create mode 100644 include/deal.II/matrix_free/cuda_hanging_nodes_internal.h diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index 4fd2ab7dad..08d6e1c83f 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -16,14 +16,19 @@ #ifndef dealii_cuda_fe_evaluation_h #define dealii_cuda_fe_evaluation_h -#include -#include +#include -#include +#ifdef DEAL_II_WITH_CUDA -#include -#include -#include +# include +# include + +# include + +# include +# include +# include +# include DEAL_II_NAMESPACE_OPEN @@ -32,18 +37,6 @@ DEAL_II_NAMESPACE_OPEN */ namespace CUDAWrappers { - namespace internal - { - template - __device__ void - resolve_hanging_nodes_shmem(Number *values, const unsigned int constr) - { - // TODO - } - } // namespace internal - - - /** * This class provides all the functions necessary to evaluate functions at * quadrature points and cell integrations. In functionality, this class is @@ -237,8 +230,8 @@ namespace CUDAWrappers values[idx] = __ldg(&src[src_idx]); if (constraint_mask) - internal::resolve_hanging_nodes_shmem( - values, constraint_mask); + internal::resolve_hanging_nodes(constraint_mask, + values); __syncthreads(); } @@ -257,9 +250,8 @@ namespace CUDAWrappers static_assert(n_components_ == 1, "This function only supports FE with one \ components"); if (constraint_mask) - internal::resolve_hanging_nodes_shmem( - values, constraint_mask); - + internal::resolve_hanging_nodes(constraint_mask, + values); const unsigned int idx = (threadIdx.x % n_q_points_1d) + @@ -456,3 +448,5 @@ namespace CUDAWrappers DEAL_II_NAMESPACE_CLOSE #endif + +#endif diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h new file mode 100644 index 0000000000..39a8589535 --- /dev/null +++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h @@ -0,0 +1,985 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__cuda_hanging_nodes_internal_h +#define dealii__cuda_hanging_nodes_internal_h + +#include + +#ifdef DEAL_II_WITH_CUDA + +# include + +# include +# include + +# include +# include + +DEAL_II_NAMESPACE_OPEN + +namespace CUDAWrappers +{ + namespace internal + { + /** + * This class creates the mask used in the treatment of hanging nodes in + * CUDAWrappers::MatrixFree. + * The implementation of this class is explained in Section 3 of + * Matrix-Free Finite-Element Computations On Graphics Processors With + * Adaptively Refined Unstructured Meshes by Karl Ljungkvist, + * SpringSim-HPC, 2017 April 23-26. + */ + template + class HangingNodes + { + public: + /** + * Constructor. + */ + HangingNodes(unsigned int fe_degree, + const DoFHandler & dof_handler, + const std::vector &lexicographic_mapping); + + /** + * Compute the value of the constraint mask for a given cell. + */ + template + void + setup_constraints(std::vector &dof_indices, + const CellIterator & cell, + unsigned int & mask) const; + + private: + /** + * Set up line-to-cell mapping for edge constraints in 3D. + */ + void + setup_line_to_cell(); + + void + rotate_subface_index(int times, unsigned int &subface_index) const; + + void + rotate_face(int times, + unsigned int n_dofs_1d, + std::vector &dofs) const; + + unsigned int + line_dof_idx(int local_line, + unsigned int dof, + unsigned int n_dofs_1d) const; + + void + transpose_face(std::vector &dofs) const; + + void + transpose_subface_index(unsigned int &subface) const; + + typedef typename DoFHandler::cell_iterator cell_iterator; + typedef + typename DoFHandler::active_cell_iterator active_cell_iterator; + const unsigned int n_raw_lines; + std::vector>> + line_to_cells; + const std::vector &lexicographic_mapping; + const unsigned int fe_degree; + const DoFHandler & dof_handler; + }; + + namespace internal + { + // TODO: use a template parameter instead of a macro +# define MAX_ELEM_DEGREE 10 + __constant__ double + constraint_weights[(MAX_ELEM_DEGREE + 1) * (MAX_ELEM_DEGREE + 1)]; + + // Here is the system for how we store constraint types in a binary mask. + // This is not a complete contradiction-free system, i.e., there are + // invalid states that we just assume that we never get. + + // If the mask is zero, there are no constraints. Then, there are three + // different fields with one bit per dimension. The first field determines + // the type, or the position of an element along each direction. The + // second field determines if there is a constrained face with that + // direction as normal. The last field determines if there is a + // constrained edge of a given pair of coordinate planes, but where + // neither of the corresponding faces are constrained (only valid in 3D). + + // The element is placed in the 'first position' along *-axis. These also + // determine which face is constrained. For example, in 2D, if + // constr_face_x and constr_type are set, then x = 0 is contrained. + constexpr unsigned int constr_type_x = 1 << 0; + constexpr unsigned int constr_type_y = 1 << 1; + constexpr unsigned int constr_type_z = 1 << 2; + + // Element has as a constraint at * = 0 or * = fe_degree face + constexpr unsigned int constr_face_x = 1 << 3; + constexpr unsigned int constr_face_y = 1 << 4; + constexpr unsigned int constr_face_z = 1 << 5; + + // Element has as a constraint at * = 0 or * = fe_degree edge + constexpr unsigned int constr_edge_xy = 1 << 6; + constexpr unsigned int constr_edge_yz = 1 << 7; + constexpr unsigned int constr_edge_zx = 1 << 8; + + template + void + setup_constraint_weigths(unsigned int fe_degree) + { + FE_Q<2> fe_q(fe_degree); + FullMatrix interpolation_matrix(fe_q.dofs_per_face, + fe_q.dofs_per_face); + fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix); + + std::vector mapping = + FETools::lexicographic_to_hierarchic_numbering<1>(FE_Q<1>(fe_degree)); + + FullMatrix mapped_matrix(fe_q.dofs_per_face, + fe_q.dofs_per_face); + for (unsigned int i = 0; i < fe_q.dofs_per_face; ++i) + for (unsigned int j = 0; j < fe_q.dofs_per_face; ++j) + mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]); + + cudaError_t error_code = + cudaMemcpyToSymbol(internal::constraint_weights, + &mapped_matrix[0][0], + sizeof(double) * fe_q.dofs_per_face * + fe_q.dofs_per_face); + AssertCuda(error_code); + } + } // namespace internal + + + + template + HangingNodes::HangingNodes( + unsigned int fe_degree, + const DoFHandler & dof_handler, + const std::vector &lexicographic_mapping) + : n_raw_lines(dof_handler.get_triangulation().n_raw_lines()) + , line_to_cells(dim == 3 ? n_raw_lines : 0) + , lexicographic_mapping(lexicographic_mapping) + , fe_degree(fe_degree) + , dof_handler(dof_handler) + { + AssertThrow( + (dim == 3) || ((fe_degree % 2) == 1), + ExcMessage( + "This function is not implemented when dim = 2 and fe_degree is even.")); + + // Set up line-to-cell mapping for edge constraints (only if dim = 3) + setup_line_to_cell(); + + internal::setup_constraint_weigths(fe_degree); + } + + + + template + void + HangingNodes::setup_line_to_cell() + {} + + + + template <> + void + HangingNodes<3>::setup_line_to_cell() + { + // In 3D, we can have DoFs on only an edge being constrained (e.g. in a + // cartesian 2x2x2 grid, where only the upper left 2 cells are refined). + // This sets up a helper data structure in the form of a mapping from + // edges (i.e. lines) to neighboring cells. + + // Mapping from an edge to which children that share that edge. + const unsigned int line_to_children[12][2] = {{0, 2}, + {1, 3}, + {0, 1}, + {2, 3}, + {4, 6}, + {5, 7}, + {4, 5}, + {6, 7}, + {0, 4}, + {1, 5}, + {2, 6}, + {3, 7}}; + + std::vector>> + line_to_inactive_cells(n_raw_lines); + + // First add active and inactive cells to their lines: + for (auto cell : dof_handler.cell_iterators()) + { + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; + ++line) + { + const unsigned int line_idx = cell->line(line)->index(); + if (cell->active()) + line_to_cells[line_idx].push_back(std::make_pair(cell, line)); + else + line_to_inactive_cells[line_idx].push_back( + std::make_pair(cell, line)); + } + } + + // Now, we can access edge-neighboring active cells on same level to also + // access of an edge to the edges "children". These are found from looking + // at the corresponding edge of children of inactive edge neighbors. + for (unsigned int line_idx = 0; line_idx < n_raw_lines; ++line_idx) + { + if ((line_to_cells[line_idx].size() > 0) && + line_to_inactive_cells[line_idx].size() > 0) + { + // We now have cells to add (active ones) and edges to which they + // should be added (inactive cells). + const cell_iterator &inactive_cell = + line_to_inactive_cells[line_idx][0].first; + const unsigned int neighbor_line = + line_to_inactive_cells[line_idx][0].second; + + for (unsigned int c = 0; c < 2; ++c) + { + const cell_iterator &child = + inactive_cell->child(line_to_children[neighbor_line][c]); + const unsigned int child_line_idx = + child->line(neighbor_line)->index(); + + // Now add all active cells + for (auto cl : line_to_cells[line_idx]) + line_to_cells[child_line_idx].push_back(cl); + } + } + } + } + + + + template + template + void + HangingNodes::setup_constraints(std::vector &dof_indices, + const CellIterator & cell, + unsigned int & mask) const + { + mask = 0; + const unsigned int n_dofs_1d = fe_degree + 1; + const unsigned int dofs_per_face = + Utilities::fixed_power(n_dofs_1d); + + std::vector neighbor_dofs(dofs_per_face); + + std::vector lex_face_mapping = + FETools::lexicographic_to_hierarchic_numbering( + FE_Q(fe_degree)); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + if ((!cell->at_boundary(face)) && + (cell->neighbor(face)->has_children() == false)) + { + const active_cell_iterator &neighbor = cell->neighbor(face); + + // Neighbor is coarser than us, i.e., face is constrained + if (neighbor->level() < cell->level()) + { + const unsigned int neighbor_face = + cell->neighbor_face_no(face); + + // Find position of face on neighbor + unsigned int subface = 0; + for (; subface < GeometryInfo::max_children_per_face; + ++subface) + if (neighbor->neighbor_child_on_subface(neighbor_face, + subface) == cell) + break; + + // Get indices to read + neighbor->face(neighbor_face)->get_dof_indices(neighbor_dofs); + + if (dim == 2) + { + if (face < 2) + { + mask |= internal::constr_face_x; + if (face == 0) + mask |= internal::constr_type_x; + if (subface == 0) + mask |= internal::constr_type_y; + } + else + { + mask |= internal::constr_face_y; + if (face == 2) + mask |= internal::constr_type_y; + if (subface == 0) + mask |= internal::constr_type_x; + } + + // Reorder neighbor_dofs and copy into faceth face of + // dof_indices + + // Offset if upper/right face + unsigned int offset = (face % 2 == 1) ? fe_degree : 0; + + for (unsigned int i = 0; i < n_dofs_1d; ++i) + { + unsigned int idx = 0; + // If X-line, i.e., if y = 0 or y = fe_degree + if (face > 1) + idx = n_dofs_1d * offset + i; + // If Y-line, i.e., if x = 0 or x = fe_degree + else + idx = n_dofs_1d * i + offset; + + dof_indices[idx] = neighbor_dofs[lex_face_mapping[i]]; + } + } + else if (dim == 3) + { + const bool transpose = !(cell->face_orientation(face)); + + int rotate = 0; + + if (cell->face_rotation(face)) + rotate -= 1; + if (cell->face_flip(face)) + rotate -= 2; + + rotate_face(rotate, n_dofs_1d, neighbor_dofs); + rotate_subface_index(rotate, subface); + + if (transpose) + { + transpose_face(neighbor_dofs); + transpose_subface_index(subface); + } + + // YZ-plane + if (face < 2) + { + mask |= internal::constr_face_x; + if (face == 0) + mask |= internal::constr_type_x; + if (subface % 2 == 0) + mask |= internal::constr_type_y; + if (subface / 2 == 0) + mask |= internal::constr_type_z; + } + // XZ-plane + else if (face < 4) + { + mask |= internal::constr_face_y; + if (face == 2) + mask |= internal::constr_type_y; + if (subface % 2 == 0) + mask |= internal::constr_type_z; + if (subface / 2 == 0) + mask |= internal::constr_type_x; + } + // XY-plane + else + { + mask |= internal::constr_face_z; + if (face == 4) + mask |= internal::constr_type_z; + if (subface % 2 == 0) + mask |= internal::constr_type_x; + if (subface / 2 == 0) + mask |= internal::constr_type_y; + } + + // Offset if upper/right/back face + unsigned int offset = (face % 2 == 1) ? fe_degree : 0; + + for (unsigned int i = 0; i < n_dofs_1d; ++i) + { + for (unsigned int j = 0; j < n_dofs_1d; ++j) + { + unsigned int idx = 0; + // If YZ-plane, i.e., if x = 0 or x = fe_degree, + // and orientation standard + if (face < 2) + idx = n_dofs_1d * n_dofs_1d * i + + n_dofs_1d * j + offset; + // If XZ-plane, i.e., if y = 0 or y = fe_degree, + // and orientation standard + else if (face < 4) + idx = n_dofs_1d * n_dofs_1d * j + + n_dofs_1d * offset + i; + // If XY-plane, i.e., if z = 0 or z = fe_degree, + // and orientation standard + else + idx = n_dofs_1d * n_dofs_1d * offset + + n_dofs_1d * i + j; + + dof_indices[idx] = + neighbor_dofs[lex_face_mapping[n_dofs_1d * i + + j]]; + } + } + } + else + ExcNotImplemented(); + } + } + } + + // In 3D we can have a situation where only DoFs on an edge are + // constrained. Append these here. + if (dim == 3) + { + // For each line on cell, which faces does it belong to, what is the + // edge mask, what is the types of the faces it belong to, and what is + // the type along the edge. + const unsigned int line_to_edge[12][4] = { + {internal::constr_face_x | internal::constr_face_z, + internal::constr_edge_zx, + internal::constr_type_x | internal::constr_type_z, + internal::constr_type_y}, + {internal::constr_face_x | internal::constr_face_z, + internal::constr_edge_zx, + internal::constr_type_z, + internal::constr_type_y}, + {internal::constr_face_y | internal::constr_face_z, + internal::constr_edge_yz, + internal::constr_type_y | internal::constr_type_z, + internal::constr_type_x}, + {internal::constr_face_y | internal::constr_face_z, + internal::constr_edge_yz, + internal::constr_type_z, + internal::constr_type_x}, + {internal::constr_face_x | internal::constr_face_z, + internal::constr_edge_zx, + internal::constr_type_x, + internal::constr_type_y}, + {internal::constr_face_x | internal::constr_face_z, + internal::constr_edge_zx, + 0, + internal::constr_type_y}, + {internal::constr_face_y | internal::constr_face_z, + internal::constr_edge_yz, + internal::constr_type_y, + internal::constr_type_x}, + {internal::constr_face_y | internal::constr_face_z, + internal::constr_edge_yz, + 0, + internal::constr_type_x}, + {internal::constr_face_x | internal::constr_face_y, + internal::constr_edge_xy, + internal::constr_type_x | internal::constr_type_y, + internal::constr_type_z}, + {internal::constr_face_x | internal::constr_face_y, + internal::constr_edge_xy, + internal::constr_type_y, + internal::constr_type_z}, + {internal::constr_face_x | internal::constr_face_y, + internal::constr_edge_xy, + internal::constr_type_x, + internal::constr_type_z}, + {internal::constr_face_x | internal::constr_face_y, + internal::constr_edge_xy, + 0, + internal::constr_type_z}}; + + for (unsigned int local_line = 0; + local_line < GeometryInfo::lines_per_cell; + ++local_line) + { + // If we don't already have a constraint for as part of a face + if (!(mask & line_to_edge[local_line][0])) + { + // For each cell which share that edge + const unsigned int line = cell->line(local_line)->index(); + for (auto edge_neighbor : line_to_cells[line]) + { + // If one of them is coarser than us + const cell_iterator neighbor_cell = edge_neighbor.first; + if (neighbor_cell->level() < cell->level()) + { + const unsigned int local_line_neighbor = + edge_neighbor.second; + mask |= line_to_edge[local_line][1] | + line_to_edge[local_line][2]; + + bool flipped = false; + if (cell->line(local_line)->vertex_index(0) == + neighbor_cell->line(local_line_neighbor) + ->vertex_index(0)) + { + // Assuming line directions match axes directions, + // we have an unflipped edge of first type + mask |= line_to_edge[local_line][3]; + } + else if (cell->line(local_line)->vertex_index(1) == + neighbor_cell->line(local_line_neighbor) + ->vertex_index(1)) + { + // We have an unflipped edge of second type + } + else if (cell->line(local_line)->vertex_index(1) == + neighbor_cell->line(local_line_neighbor) + ->vertex_index(0)) + { + // We have a flipped edge of second type + flipped = true; + } + else if (cell->line(local_line)->vertex_index(0) == + neighbor_cell->line(local_line_neighbor) + ->vertex_index(1)) + { + // We have a flipped edge of first type + mask |= line_to_edge[local_line][3]; + flipped = true; + } + else + ExcInternalError(); + + // Copy the unconstrained values + neighbor_dofs.resize(n_dofs_1d * n_dofs_1d * + n_dofs_1d); + neighbor_cell->get_dof_indices(neighbor_dofs); + + for (unsigned int i = 0; i < n_dofs_1d; ++i) + { + // Get local dof index along line + const unsigned int idx = + line_dof_idx(local_line, i, n_dofs_1d); + dof_indices[idx] = neighbor_dofs + [lexicographic_mapping[line_dof_idx( + local_line_neighbor, + flipped ? fe_degree - i : i, + n_dofs_1d)]]; + } + + // Stop looping over edge neighbors + break; + } + } + } + } + } + } + + + + template + void + HangingNodes::rotate_subface_index(int times, + unsigned int &subface_index) const + { + const unsigned int rot_mapping[4] = {2, 0, 3, 1}; + + times = times % 4; + times = times < 0 ? times + 4 : times; + for (int t = 0; t < times; ++t) + subface_index = rot_mapping[subface_index]; + } + + + + template + void + HangingNodes::rotate_face( + int times, + unsigned int n_dofs_1d, + std::vector &dofs) const + { + const unsigned int rot_mapping[4] = {2, 0, 3, 1}; + + times = times % 4; + times = times < 0 ? times + 4 : times; + + std::vector copy(dofs.size()); + for (int t = 0; t < times; ++t) + { + std::swap(copy, dofs); + + // Vertices + for (unsigned int i = 0; i < 4; ++i) + dofs[rot_mapping[i]] = copy[i]; + + // Edges + const unsigned int n_int = n_dofs_1d - 2; + unsigned int offset = 4; + for (unsigned int i = 0; i < n_int; ++i) + { + // Left edge + dofs[offset + i] = copy[offset + 2 * n_int + (n_int - 1 - i)]; + // Right edge + dofs[offset + n_int + i] = + copy[offset + 3 * n_int + (n_int - 1 - i)]; + // Bottom edge + dofs[offset + 2 * n_int + i] = copy[offset + n_int + i]; + // Top edge + dofs[offset + 3 * n_int + i] = copy[offset + i]; + } + + // Interior points + offset += 4 * n_int; + + for (unsigned int i = 0; i < n_int; ++i) + for (unsigned int j = 0; j < n_int; ++j) + dofs[offset + i * n_int + j] = + copy[offset + j * n_int + (n_int - 1 - i)]; + } + } + + + + template + unsigned int + HangingNodes::line_dof_idx(int local_line, + unsigned int dof, + unsigned int n_dofs_1d) const + { + unsigned int x, y, z; + + if (local_line < 8) + { + x = + (local_line % 4 == 0) ? 0 : (local_line % 4 == 1) ? fe_degree : dof; + y = + (local_line % 4 == 2) ? 0 : (local_line % 4 == 3) ? fe_degree : dof; + z = (local_line / 4) * fe_degree; + } + else + { + x = ((local_line - 8) % 2) * fe_degree; + y = ((local_line - 8) / 2) * fe_degree; + z = dof; + } + + return n_dofs_1d * n_dofs_1d * z + n_dofs_1d * y + x; + } + + + + template + void + HangingNodes::transpose_face( + std::vector &dofs) const + { + const std::vector copy(dofs); + + // Vertices + dofs[1] = copy[2]; + dofs[2] = copy[1]; + + // Edges + const unsigned int n_int = fe_degree - 1; + unsigned int offset = 4; + for (unsigned int i = 0; i < n_int; ++i) + { + // Right edge + dofs[offset + i] = copy[offset + 2 * n_int + i]; + // Left edge + dofs[offset + n_int + i] = copy[offset + 3 * n_int + i]; + // Bottom edge + dofs[offset + 2 * n_int + i] = copy[offset + i]; + // Top edge + dofs[offset + 3 * n_int + i] = copy[offset + n_int + i]; + } + + // Interior + offset += 4 * n_int; + for (unsigned int i = 0; i < n_int; ++i) + for (unsigned int j = 0; j < n_int; ++j) + dofs[offset + i * n_int + j] = copy[offset + j * n_int + i]; + } + + + + template + void + HangingNodes::transpose_subface_index(unsigned int &subface) const + { + if (subface == 1) + subface = 2; + else if (subface == 2) + subface = 1; + } + + + + //-------------------------------------------------------------------------// + // Functions for resolving the hanging node constraints // + //-------------------------------------------------------------------------// + namespace internal + { + template + __device__ inline unsigned int + index2(unsigned int i, unsigned int j) + { + return i + size * j; + } + + + + template + __device__ inline unsigned int + index3(unsigned int i, unsigned int j, unsigned int k) + { + return i + size * j + size * size * k; + } + + + + template + __device__ inline void + interpolate_boundary_2d(const unsigned int constr, Number *values) + { + const unsigned int x_idx = threadIdx.x % (fe_degree + 1); + const unsigned int y_idx = threadIdx.y; + + const unsigned int this_type = + (direction == 0) ? internal::constr_type_x : internal::constr_type_y; + + if (constr & (((direction == 0) ? internal::constr_face_y : 0) | + ((direction == 1) ? internal::constr_face_x : 0))) + { + const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx; + + __syncthreads(); + + Number t = 0; + const bool flag = + ((direction == 0) && + ((constr & internal::constr_type_y) ? (y_idx == 0) : + (y_idx == fe_degree))) || + ((direction == 1) && + ((constr & internal::constr_type_x) ? (x_idx == 0) : + (x_idx == fe_degree))); + + if (flag) + { + const bool type = constr & this_type; + + if (type) + { + for (unsigned int i = 0; i <= fe_degree; ++i) + { + const unsigned int real_idx = + (direction == 0) ? index2(i, y_idx) : + index2(x_idx, i); + + const Number w = + transpose ? + internal::constraint_weights[i * (fe_degree + 1) + + interp_idx] : + internal::constraint_weights[interp_idx * + (fe_degree + 1) + + i]; + t += w * values[real_idx]; + } + } + else + { + for (unsigned int i = 0; i <= fe_degree; ++i) + { + const unsigned int real_idx = + (direction == 0) ? index2(i, y_idx) : + index2(x_idx, i); + + const Number w = + transpose ? + internal::constraint_weights[(fe_degree - i) * + (fe_degree + 1) + + fe_degree - + interp_idx] : + internal::constraint_weights[(fe_degree - + interp_idx) * + (fe_degree + 1) + + fe_degree - i]; + t += w * values[real_idx]; + } + } + } + + __syncthreads(); + + if (flag) + values[index2(x_idx, y_idx)] = t; + } + } + + + + template + __device__ inline void + interpolate_boundary_3d(const unsigned int constr, Number *values) + { + const unsigned int x_idx = threadIdx.x % (fe_degree + 1); + const unsigned int y_idx = threadIdx.y; + const unsigned int z_idx = threadIdx.z; + + const unsigned int this_type = + (direction == 0) ? internal::constr_type_x : + (direction == 1) ? internal::constr_type_y : + internal::constr_type_z; + const unsigned int face1_type = + (direction == 0) ? internal::constr_type_y : + (direction == 1) ? internal::constr_type_z : + internal::constr_type_x; + const unsigned int face2_type = + (direction == 0) ? internal::constr_type_z : + (direction == 1) ? internal::constr_type_x : + internal::constr_type_y; + + // If computing in x-direction, need to match against constr_face_y or + // constr_face_z + const unsigned int face1 = (direction == 0) ? internal::constr_face_y : + (direction == 1) ? + internal::constr_face_z : + internal::constr_face_x; + const unsigned int face2 = (direction == 0) ? internal::constr_face_z : + (direction == 1) ? + internal::constr_face_x : + internal::constr_face_y; + const unsigned int edge = (direction == 0) ? internal::constr_edge_yz : + (direction == 1) ? + internal::constr_edge_zx : + internal::constr_edge_xy; + + if (constr & (face1 | face2 | edge)) + { + const unsigned int interp_idx = + (direction == 0) ? x_idx : (direction == 1) ? y_idx : z_idx; + const unsigned int face1_idx = + (direction == 0) ? y_idx : (direction == 1) ? z_idx : x_idx; + const unsigned int face2_idx = + (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx; + + __syncthreads(); + + Number t = 0; + const bool on_face1 = (constr & face1_type) ? + (face1_idx == 0) : + (face1_idx == fe_degree); + const bool on_face2 = (constr & face2_type) ? + (face2_idx == 0) : + (face2_idx == fe_degree); + const bool flag = (((constr & face1) && on_face1) || + ((constr & face2) && on_face2) || + ((constr & edge) && on_face1 && on_face2)); + + if (flag) + { + const bool type = constr & this_type; + if (type) + { + for (unsigned int i = 0; i <= fe_degree; ++i) + { + const unsigned int real_idx = + (direction == 0) ? + index3(i, y_idx, z_idx) : + (direction == 1) ? + index3(x_idx, i, z_idx) : + index3(x_idx, y_idx, i); + + const Number w = + transpose ? + internal::constraint_weights[i * (fe_degree + 1) + + interp_idx] : + internal::constraint_weights[interp_idx * + (fe_degree + 1) + + i]; + t += w * values[real_idx]; + } + } + else + { + for (unsigned int i = 0; i <= fe_degree; ++i) + { + const unsigned int real_idx = + (direction == 0) ? + index3(i, y_idx, z_idx) : + (direction == 1) ? + index3(x_idx, i, z_idx) : + index3(x_idx, y_idx, i); + + const Number w = + transpose ? + internal::constraint_weights[(fe_degree - i) * + (fe_degree + 1) + + fe_degree - + interp_idx] : + internal::constraint_weights[(fe_degree - + interp_idx) * + (fe_degree + 1) + + fe_degree - i]; + t += w * values[real_idx]; + } + } + } + + __syncthreads(); + + if (flag) + values[index3(x_idx, y_idx, z_idx)] = t; + } + } + } // namespace internal + + + + /** + * This function resolves the hanging nodes using tensor product. + * + * The implementation of this class is explained in Matrix-Free + * Finite-Element Computations On Graphics Processors With Adaptively + * Refined Unstructured Meshes by Karl Ljungkvist, SpringSim-HPC, 2017 + * April 23-26. + */ + template + __device__ void + resolve_hanging_nodes(const unsigned int constr, Number *values) + { + if (dim == 2) + { + internal::interpolate_boundary_2d(constr, + values); + internal::interpolate_boundary_2d(constr, + values); + } + else if (dim == 3) + { + // Interpolate y and z faces (x-direction) + internal::interpolate_boundary_3d(constr, + values); + // Interpolate x and z faces (y-direction) + internal::interpolate_boundary_3d(constr, + values); + // Interpolate x and y faces (z-direction) + internal::interpolate_boundary_3d(constr, + values); + } + } + } // namespace internal +} // namespace CUDAWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h index c6979dba52..c9615aef3e 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -28,6 +28,7 @@ # include +# include # include # include @@ -114,7 +115,6 @@ namespace CUDAWrappers /** * Helper class to (re)initialize MatrixFree object. */ - // TODO for now does not support hanging_nodes template class ReinitHelper { @@ -125,8 +125,9 @@ namespace CUDAWrappers const FiniteElement &fe, const Quadrature<1> & quad, const ::dealii::internal::MatrixFreeFunctions::ShapeInfo - & shape_info, - const UpdateFlags &update_flags); + & shape_info, + const DoFHandler &dof_handler, + const UpdateFlags & update_flags); void setup_color_arrays(const unsigned int n_colors); @@ -160,6 +161,7 @@ namespace CUDAWrappers const unsigned int q_points_per_cell; const UpdateFlags & update_flags; const unsigned int padding_length; + HangingNodes hanging_nodes; }; @@ -171,8 +173,9 @@ namespace CUDAWrappers const FiniteElement &fe, const Quadrature<1> & quad, const ::dealii::internal::MatrixFreeFunctions::ShapeInfo - & shape_info, - const UpdateFlags &update_flags) + & shape_info, + const DoFHandler &dof_handler, + const UpdateFlags & update_flags) : data(data) , fe_degree(data->fe_degree) , dofs_per_cell(data->dofs_per_cell) @@ -185,6 +188,7 @@ namespace CUDAWrappers , lexicographic_inv(shape_info.lexicographic_numbering) , update_flags(update_flags) , padding_length(data->get_padding_length()) + , hanging_nodes(fe_degree, dof_handler, lexicographic_inv) { local_dof_indices.resize(data->dofs_per_cell); lexicographic_dof_indices.resize(dofs_per_cell); @@ -275,9 +279,14 @@ namespace CUDAWrappers { cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) lexicographic_dof_indices[i] = local_dof_indices[lexicographic_inv[i]]; + hanging_nodes.setup_constraints(lexicographic_dof_indices, + cell, + constraint_mask_host[cell_id]); + memcpy(&local_to_global_host[cell_id * padding_length], lexicographic_dof_indices.data(), dofs_per_cell * sizeof(unsigned int)); @@ -323,9 +332,7 @@ namespace CUDAWrappers // Local-to-global mapping if (data->parallelization_scheme == MatrixFree::parallel_over_elem) - internal::transpose_in_place(local_to_global_host, - n_cells, - padding_length); + transpose_in_place(local_to_global_host, n_cells, padding_length); alloc_and_copy(&data->local_to_global[color], local_to_global_host, @@ -336,9 +343,7 @@ namespace CUDAWrappers { if (data->parallelization_scheme == MatrixFree::parallel_over_elem) - internal::transpose_in_place(q_points_host, - n_cells, - padding_length); + transpose_in_place(q_points_host, n_cells, padding_length); alloc_and_copy(&data->q_points[color], q_points_host, @@ -350,7 +355,7 @@ namespace CUDAWrappers { if (data->parallelization_scheme == MatrixFree::parallel_over_elem) - internal::transpose_in_place(JxW_host, n_cells, padding_length); + transpose_in_place(JxW_host, n_cells, padding_length); alloc_and_copy(&data->JxW[color], JxW_host, n_cells * padding_length); } @@ -362,18 +367,18 @@ namespace CUDAWrappers // are together, etc., i.e., reorder indices from // cell_id*q_points_per_cell*dim*dim + q*dim*dim +i to // i*q_points_per_cell*n_cells + cell_id*q_points_per_cell+q - internal::transpose_in_place(inv_jacobian_host, - padding_length * n_cells, - dim * dim); + transpose_in_place(inv_jacobian_host, + padding_length * n_cells, + dim * dim); // Transpose second time means we get the following index order: // q*n_cells*dim*dim + i*n_cells + cell_id which is good for an // element-level parallelization if (data->parallelization_scheme == MatrixFree::parallel_over_elem) - internal::transpose_in_place(inv_jacobian_host, - n_cells * dim * dim, - padding_length); + transpose_in_place(inv_jacobian_host, + n_cells * dim * dim, + padding_length); alloc_and_copy(&data->inv_jacobian[color], inv_jacobian_host, @@ -546,7 +551,7 @@ namespace CUDAWrappers cells_per_block = cells_per_block_shmem(dim, fe_degree); internal::ReinitHelper helper( - this, mapping, fe, quad, shape_info, update_flags); + this, mapping, fe, quad, shape_info, dof_handler, update_flags); // Create a graph coloring using CellFilter = @@ -554,8 +559,6 @@ namespace CUDAWrappers CellFilter begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); CellFilter end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); - using fun_type = - std::function(CellFilter const &)>; const auto fun = [&](const CellFilter &filter) { return internal::get_conflict_indices(filter, constraints); }; diff --git a/source/matrix_free/CMakeLists.txt b/source/matrix_free/CMakeLists.txt index 6f09ed5790..c1b9423f6e 100644 --- a/source/matrix_free/CMakeLists.txt +++ b/source/matrix_free/CMakeLists.txt @@ -40,10 +40,6 @@ IF(DEAL_II_WITH_CUDA) SET(_inst ${_inst} ) - FILE(GLOB _header - ${CMAKE_SOURCE_DIR}/include/deal.II/matrix_free/*.cuh - ${_header} - ) ENDIF() DEAL_II_ADD_LIBRARY(obj_matrix_free OBJECT ${_src} ${_header} ${_inst}) -- 2.39.5