From: Daniel Arndt <arndtd@ornl.gov>
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.
 
 
 <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 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<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 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<dim, Number>::Data &)
+  {}
+
   static const unsigned int n_q_points =
     dealii::Utilities::pow(fe_degree + 1, dim);