From 49022f8c3e06fd3bb307aa6f3bbb265df84f27a9 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 21 Aug 2024 14:18:36 -0400 Subject: [PATCH] Fix hanging node constraints --- include/deal.II/matrix_free/tools.h | 43 +++++++++++++++++-- .../matrix_free_kokkos/compute_diagonal_01.cc | 12 ++---- .../compute_diagonal_util.h | 10 ++--- 3 files changed, 47 insertions(+), 18 deletions(-) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index 21531c3785..a6fd34cde2 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -25,6 +25,7 @@ #include #include + DEAL_II_NAMESPACE_OPEN /** @@ -1375,7 +1376,7 @@ namespace MatrixFreeTools {} KOKKOS_FUNCTION void - operator()(const unsigned int /*cell*/, + operator()(const unsigned int cell, const typename Portable::MatrixFree::Data *gpu_data, Portable::SharedData *shared_data, const Number *, @@ -1384,15 +1385,31 @@ namespace MatrixFreeTools Portable::FEEvaluation fe_eval( gpu_data, shared_data); constexpr int dofs_per_cell = decltype(fe_eval)::tensor_dofs_per_cell; - double diagonal[dofs_per_cell] = {}; + Number diagonal[dofs_per_cell] = {}; for (unsigned int i = 0; i < dofs_per_cell; ++i) { Kokkos::parallel_for( Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell), [&](int j) { fe_eval.submit_dof_value(i == j ? 1 : 0, j); }); + + Portable::internal:: + resolve_hanging_nodes( + shared_data->team_member, + gpu_data->constraint_weights, + gpu_data->constraint_mask(cell), + Kokkos::subview(shared_data->values, Kokkos::ALL, 0)); + fe_eval.evaluate(m_evaluation_flags); fe_eval.apply_for_each_quad_point(m_quad_operation); fe_eval.integrate(m_integration_flags); + + Portable::internal:: + resolve_hanging_nodes( + shared_data->team_member, + gpu_data->constraint_weights, + gpu_data->constraint_mask(cell), + Kokkos::subview(shared_data->values, Kokkos::ALL, 0)); + Kokkos::single(Kokkos::PerTeam(shared_data->team_member), [&] { diagonal[i] = fe_eval.get_dof_value(i); }); } @@ -1401,7 +1418,27 @@ namespace MatrixFreeTools for (unsigned int i = 0; i < dofs_per_cell; ++i) fe_eval.submit_dof_value(diagonal[i], i); }); - fe_eval.distribute_local_to_global(dst); + + // We need to do the same as distribute_local_to_global but without + // constraints since we have already taken care of them earlier + if (gpu_data->use_coloring) + { + Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member, + dofs_per_cell), + [&](const int &i) { + dst[gpu_data->local_to_global(cell, i)] += + shared_data->values(i, 0); + }); + } + else + { + Kokkos::parallel_for( + Kokkos::TeamThreadRange(shared_data->team_member, dofs_per_cell), + [&](const int &i) { + Kokkos::atomic_add(&dst[gpu_data->local_to_global(cell, i)], + shared_data->values(i, 0)); + }); + } }; static constexpr unsigned int n_local_dofs = QuadOperation::n_local_dofs; diff --git a/tests/matrix_free_kokkos/compute_diagonal_01.cc b/tests/matrix_free_kokkos/compute_diagonal_01.cc index 2865a5069a..6df68126c4 100644 --- a/tests/matrix_free_kokkos/compute_diagonal_01.cc +++ b/tests/matrix_free_kokkos/compute_diagonal_01.cc @@ -30,12 +30,10 @@ test() GridGenerator::hyper_cube(tria, -numbers::PI / 2, numbers::PI / 2); tria.refine_global(2); - // FIXME - // for (auto &cell : tria.active_cell_iterators()) - // if (cell->is_active() && cell->is_locally_owned() && - // cell->center()[0] < 0.0) - // tria.begin_active()->set_refine_flag(); - // tria.execute_coarsening_and_refinement(); + for (auto &cell : tria.active_cell_iterators()) + if (cell->is_active() && cell->is_locally_owned() && + cell->center()[0] < 0.0) + tria.execute_coarsening_and_refinement(); const FE_Q fe_q(fe_degree); const FESystem fe(fe_q, n_components); @@ -54,8 +52,6 @@ test() constraint); constraint.close(); - // constraint.print(std::cout); - typename Portable::MatrixFree::AdditionalData additional_data; additional_data.mapping_update_flags = update_values | update_gradients; diff --git a/tests/matrix_free_kokkos/compute_diagonal_util.h b/tests/matrix_free_kokkos/compute_diagonal_util.h index ae71b3f1c4..12017ebdc6 100644 --- a/tests/matrix_free_kokkos/compute_diagonal_util.h +++ b/tests/matrix_free_kokkos/compute_diagonal_util.h @@ -120,7 +120,6 @@ public: Assert(diagonal_global_host[i] > 0, ExcMessage("Diagonal non-positive at position " + std::to_string(i))); - // std::cout << i << ": " << diagonal_global_host[i] << '\n'; } } @@ -151,12 +150,9 @@ public: { if (!constraints.is_constrained(i)) { - if (std::abs(A_ref(i, i) - diagonal_global_host(i)) > 1e-6) - // std::cout << i << ": " << A_ref(i,i) << " should be " << - // diagonal_global_host(i) << '\n'; - Assert(std::abs(A_ref(i, i) - diagonal_global_host(i)) < 1e-6, - ExcMessage("Wrong diagonal entry at position " + - std::to_string(i))); + Assert(std::abs(A_ref(i, i) - diagonal_global_host(i)) < 1e-6, + ExcMessage("Wrong diagonal entry at position " + + std::to_string(i))); } } } -- 2.39.5