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 6a4c370d305a0cb66851a85e0ac46d3a254007ee..5c3ed6104a27ea71645e72af6a12a4586f4eb547 100644 (file)
--- 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.
 
 
 <a name="step-64-extensions"></a>
 <h3> Possibilities for extensions </h3>
 
-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 a6fd34cde20cbc17b0073a97f1d15294e7ffab5f..e78510062d6ad0ef21f19591b8e471e8da2dbcc9 100644 (file)
--- a/include/deal.II/matrix_free/tools.h
+++ b/include/deal.II/matrix_free/tools.h
@@ -1384,6 +1384,8 @@ namespace MatrixFreeTools
     {
       Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> 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<dim, fe_degree, Number, QuadOperation> cell_action(
       quad_operation, evaluation_flags, integration_flags);
-    // diagonal_global.zero_out_ghost_values();
     LinearAlgebra::distributed::Vector<Number, MemorySpace> 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 12017ebdc6e066e71fd37e0a98c3a968fe9aea98..d4d04f184af191f535ee3e1bf55026a69e026512 100644 (file)
--- 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<dim, Number>::Data &)
+  {}
+
   static const unsigned int n_q_points =
     dealii::Utilities::pow(fe_degree + 1, dim);