From: Daniel Arndt Date: Fri, 23 Aug 2024 20:35:28 +0000 (-0400) Subject: Update interface to allow setting cell and data X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7713cb5061db2cf5ae975bd83aa355e3bef724d;p=dealii.git Update interface to allow setting cell and data --- diff --git a/examples/step-64/doc/results.dox b/examples/step-64/doc/results.dox index 6a4c370d30..5c3ed6104a 100644 --- a/examples/step-64/doc/results.dox +++ b/examples/step-64/doc/results.dox @@ -7,48 +7,37 @@ itself, we just show the expected output here: Cycle 0 Number of active cells: 8 Number of degrees of freedom: 343 - Solved in 27 iterations. + Solved in 10 iterations. solution norm: 0.0205439 Cycle 1 Number of active cells: 64 Number of degrees of freedom: 2197 - Solved in 60 iterations. + Solved in 14 iterations. solution norm: 0.0205269 Cycle 2 Number of active cells: 512 Number of degrees of freedom: 15625 - Solved in 114 iterations. + Solved in 29 iterations. solution norm: 0.0205261 Cycle 3 Number of active cells: 4096 Number of degrees of freedom: 117649 - Solved in 227 iterations. + Solved in 58 iterations. solution norm: 0.0205261 @endcode One can make two observations here: First, the norm of the numerical solution converges, presumably to the norm of the exact (but unknown) -solution. And second, the number of iterations roughly doubles with -each refinement of the mesh. (This is in keeping with the expectation -that the number of CG iterations grows with the square root of the -condition number of the matrix; and that we know that the condition -number of the matrix of a second-order differential operation grows -like ${\cal O}(h^{-2})$.) This is of course rather inefficient, as an -optimal solver would have a number of iterations that is independent -of the size of the problem. But having such a solver would require -using a better preconditioner than the identity matrix we have used here. +solution. And second, the number of iterations keep increasing with +each refinement of the mesh since we are only using a preconditioner based on +the diagonal of the operator.

Possibilities for extensions

-Currently, this program uses no preconditioner at all. This is mainly -since constructing an efficient matrix-free preconditioner is -non-trivial. However, simple choices just requiring the diagonal of -the corresponding matrix are good candidates and these can be computed -in a matrix-free way as well. Alternatively, and maybe even better, -one could extend the tutorial to use multigrid with Chebyshev +One could extend the tutorial to use multigrid with Chebyshev smoothers similar to step-37. diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index a6fd34cde2..e78510062d 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -1384,6 +1384,8 @@ namespace MatrixFreeTools { Portable::FEEvaluation fe_eval( gpu_data, shared_data); + m_quad_operation.set_matrix_free_data(*gpu_data); + m_quad_operation.set_cell(cell); constexpr int dofs_per_cell = decltype(fe_eval)::tensor_dofs_per_cell; Number diagonal[dofs_per_cell] = {}; for (unsigned int i = 0; i < dofs_per_cell; ++i) @@ -1444,7 +1446,7 @@ namespace MatrixFreeTools static constexpr unsigned int n_local_dofs = QuadOperation::n_local_dofs; private: - const QuadOperation &m_quad_operation; + mutable QuadOperation m_quad_operation; const EvaluationFlags::EvaluationFlags m_evaluation_flags; const EvaluationFlags::EvaluationFlags m_integration_flags; }; @@ -1483,7 +1485,6 @@ namespace MatrixFreeTools CellAction cell_action( quad_operation, evaluation_flags, integration_flags); - // diagonal_global.zero_out_ghost_values(); LinearAlgebra::distributed::Vector dummy; matrix_free.cell_loop(cell_action, dummy, diagonal_global); diff --git a/tests/matrix_free_kokkos/compute_diagonal_util.h b/tests/matrix_free_kokkos/compute_diagonal_util.h index 12017ebdc6..d4d04f184a 100644 --- a/tests/matrix_free_kokkos/compute_diagonal_util.h +++ b/tests/matrix_free_kokkos/compute_diagonal_util.h @@ -58,6 +58,15 @@ public: fe_eval->submit_gradient(fe_eval->get_gradient(q_point), q_point); } + DEAL_II_HOST_DEVICE + void + set_cell(int) + {} + + DEAL_II_HOST_DEVICE void + set_matrix_free_data(const Portable::MatrixFree::Data &) + {} + static const unsigned int n_q_points = dealii::Utilities::pow(fe_degree + 1, dim);