`MatrixFreeTools::compute_diagonal()` and `compute_matrix()` for face and boundary integrals #17029

peterrum · 2024-05-16T12:33:48Z

This feature was already requested multiple times (also on the mailing list).

bangerth · 2024-05-16T12:46:39Z

Did you mean to perhaps merge some of the commits? ;-)

peterrum · 2024-05-16T12:38:14Z

include/deal.II/matrix_free/tools.h

+            for (unsigned int i = 0; i < phi_m.dofs_per_cell; ++i)
+              {
+                helper_m.prepare_basis_vector(i);
+                helper_p.zero_basis_vector();
+                face_operation(phi_m, phi_p);
+                helper_m.submit();
+              }
+
+            helper_m.distribute_local_to_global(diagonal_global_components);
+
+            for (unsigned int i = 0; i < phi_p.dofs_per_cell; ++i)
+              {
+                helper_m.zero_basis_vector();
+                helper_p.prepare_basis_vector(i);
+                face_operation(phi_m, phi_p);
+                helper_p.submit();
+              }
+
+            helper_p.distribute_local_to_global(diagonal_global_components);


This is how the diagonal for the face integral is determined. Currently, there is a limitation that not arbitrary constraints are supported (which should be fine as a first step, since in the DG constraints you normally don't have constraints).

peterrum · 2024-05-16T12:39:15Z

include/deal.II/matrix_free/tools.h

+            for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+              {
+                helper.prepare_basis_vector(i);
+                boundary_operation(phi);
+                helper.submit();
+              }
+
+            helper.distribute_local_to_global(diagonal_global_components);


Similarly, the diagonal entries related to boundary values are determined here (without any limitations).

peterrum · 2024-05-16T12:40:01Z

include/deal.II/matrix_free/tools.h

+          return;
+
+        // shortcut for faces containing purely FE_Nothing
+        if (internal::is_fe_nothing<false>(matrix_free,


We ignore cells/faces with FE_Nothing attached to them.

peterrum · 2024-05-16T12:44:59Z

include/deal.II/matrix_free/tools.h

+        internal::create_new_affine_constraints_if_needed(
+          matrix, constraints_in, constraints_for_matrix);
+
+      const auto batch_operation =


This is the relevant function that computes the stiffness matrices. It is generalized in the sense that it works on a list of FEEvaluationData (blocks). This is useful, since it allows to have one code for computing stiffness matrices related to cell (1 FEEvaluation), face (2 FEFaceEvaluation), and boundary integrals (1 FEFaceEvaluation). In the future, one should also be able to pass in multiple FEEvaluation (as needed in the Navier-Stokes case: velocity and pressure).

peterrum · 2024-05-16T12:45:51Z

include/deal.II/matrix_free/tools.h

+  namespace internal
+  {
+    template <int dim, typename Number, bool is_face_>
+    class ComputeMatrixScratchData
+    {
+    public:
+      using FEEvalType = FEEvaluationData<dim, Number, is_face_>;
+
+      std::vector<unsigned int> dof_numbers;
+      std::vector<unsigned int> quad_numbers;
+      std::vector<unsigned int> first_selected_components;
+      std::vector<unsigned int> batch_type;
+      static const bool         is_face = is_face_;
+
+      std::function<std::vector<std::unique_ptr<FEEvalType>>(
+        const std::pair<unsigned int, unsigned int> &)>
+        op_create;
+      std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
+                         const unsigned int)>
+        op_reinit;
+      std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
+        op_compute;
+    };
+  } // namespace internal


This function is configured by this struct.

peterrum · 2024-05-16T12:47:49Z

include/deal.II/matrix_free/tools.h

+
+    data_face.op_create =
+      [&](const std::pair<unsigned int, unsigned int> &range) {
+        std::vector<
+          std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+          phi;
+
+        if (!internal::is_fe_nothing<true>(matrix_free,
+                                           range,
+                                           dof_no,
+                                           quad_no,
+                                           first_selected_component,
+                                           fe_degree,
+                                           n_q_points_1d,
+                                           true) &&
+            !internal::is_fe_nothing<true>(matrix_free,
+                                           range,
+                                           dof_no,
+                                           quad_no,
+                                           first_selected_component,
+                                           fe_degree,
+                                           n_q_points_1d,
+                                           false))
          {
-            integrator.reinit(cell);
+            phi.emplace_back(
+              std::make_unique<FEFaceEvalType>(matrix_free,
+                                               range,
+                                               true,
+                                               dof_no,
+                                               quad_no,
+                                               first_selected_component));
+            phi.emplace_back(
+              std::make_unique<FEFaceEvalType>(matrix_free,
+                                               range,
+                                               false,
+                                               dof_no,
+                                               quad_no,
+                                               first_selected_component));
+          }

-            const unsigned int n_filled_lanes =
-              matrix_free.n_active_entries_per_cell_batch(cell);
+        return phi;
+      };

-            for (unsigned int v = 0; v < n_filled_lanes; ++v)
-              matrices[v] = 0.0;
+    data_face.op_reinit = [](auto &phi, const unsigned batch) {
+      static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+      static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
+    };

-            for (unsigned int j = 0; j < dofs_per_cell; ++j)
-              {
-                for (unsigned int i = 0; i < dofs_per_cell; ++i)
-                  integrator.begin_dof_values()[i] =
-                    static_cast<Number>(i == j);
+    if (face_operation)
+      data_face.op_compute = [&](auto &phi) {
+        face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
+                       static_cast<FEFaceEvalType &>(*phi[1]));
+      };
+
+    internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+      data_boundary;
+
+    data_boundary.dof_numbers               = {dof_no};
+    data_boundary.quad_numbers              = {dof_no};
+    data_boundary.first_selected_components = {first_selected_component};
+    data_boundary.batch_type                = {1};
+
+    data_boundary
+      .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
+      std::vector<
+        std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+        phi;
+
+      if (!internal::is_fe_nothing<true>(matrix_free,
+                                         range,
+                                         dof_no,
+                                         quad_no,
+                                         first_selected_component,
+                                         fe_degree,
+                                         n_q_points_1d,
+                                         true))
+        phi.emplace_back(std::make_unique<FEFaceEvalType>(
+          matrix_free, range, true, dof_no, quad_no, first_selected_component));
+
+      return phi;
+    };
+
+    data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
+      static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+    };
+
+    if (boundary_operation)
+      data_boundary.op_compute = [&](auto &phi) {
+        boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
+      };


Here, the struct is filled for the cell/face/boundary integral case.

peterrum · 2024-05-16T12:50:38Z

include/deal.II/matrix_free/tools.h

+    template <int dim, typename Number, bool is_face_>
+    class ComputeMatrixScratchData


Currently, this class and the corresponding compute_matrix() function are internal. I apply these in different contexts to see the applicability. Once the interface seems to be general enough, I would suggest to make these public so that users have full flexibility.

include/deal.II/matrix_free/tools.h

peterrum · 2024-05-18T17:12:56Z

The failing test is unrelated.

kronbichler

This looks OK, thank you for tackling it. I have a few questions below. If I read the code correctly, there is no change (apart from empty std::functions) in the diagonal code for cells before, and the main extension is in the matrix creation in terms of multiple evaluators (op.create()).

doc/news/changes/minor/20240516Munch

include/deal.II/lac/full_matrix.h

include/deal.II/matrix_free/tools.h

kronbichler · 2024-06-05T11:11:30Z

include/deal.II/matrix_free/tools.h

+      const auto init_data = dealii::internal::
+        extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
+          matrix_free,
+          dof_no,
+          first_selected_component,
+          quad_no,
+          fe_degree,
+          static_n_q_points,
+          active_fe_index,
+          numbers::invalid_unsigned_int /*active_quad_index*/,
+          numbers::invalid_unsigned_int /*face_type*/);


Is this the simplest way to get the case of FE_Nothing? I would have looked up the active fe index and then checked the matrix_free.get_dof_handler(dof_no).get_fe_collection()[active_fe_index]? Is there something not working, e.g. for the faces?

The reason why we do this is that we did not copy the (somewhat complicated) logic in:

dealii/include/deal.II/matrix_free/fe_evaluation.h

Lines 2291 to 2297 in f8df7d6

init_data.active_fe_index =

fe_degree != numbers::invalid_unsigned_int ?

init_data.dof_info->fe_index_from_degree(first_selected_component,

fe_degree) :

(active_fe_index_given != numbers::invalid_unsigned_int ?

active_fe_index_given :

0);

If we would copy these lines, I am afraid that the codes will diverge sometime. @kronbichler What is your opinion on this?

So in other words, the active_fe_index variable you have at your disposal here is not the index we need to check within the actual FECollection class in the DoFHandler, right? In that case, I think I agree with the code here, even though it is not pretty.

include/deal.II/matrix_free/tools.h

peterrum · 2024-06-07T10:35:24Z

@kronbichler Thank you for the review.

This looks OK, thank you for tackling it. I have a few questions below. If I read the code correctly, there is no change (apart from empty std::functions) in the diagonal code for cells before, and the main extension is in the matrix creation in terms of multiple evaluators (op.create()).

The main new feature of this PR is that is to compute the face/boundary contributions for the computation of the diagonal and matrix.

That one can handle now multiple evaluators in the case of the matrix is a side effect of the extensive refactoring (now cell/face/boundary integrals take the same code path). One should do similar changes to the computation of the diagonal (and I am planing to do this). But, there it is more complex, since we resolve the constraints there on-the-fly (with the main limitation that the same constraints are applied from left and right) without using AffineConstraints.

…ompute_matrix()

peterrum added Matrix-free ready to test labels May 16, 2024

peterrum commented May 16, 2024

View reviewed changes

peterrum force-pushed the mf_tools_new branch from 8064fd4 to 076f7e7 Compare May 16, 2024 14:13

peterrum mentioned this pull request May 16, 2024

Use new dealii::MatrixFreeTools functions exadg/exadg#661

Open

peterrum force-pushed the mf_tools_new branch from 076f7e7 to ca75db7 Compare May 16, 2024 15:40

mschreter reviewed May 16, 2024

View reviewed changes

include/deal.II/matrix_free/tools.h Outdated Show resolved Hide resolved

include/deal.II/matrix_free/tools.h Outdated Show resolved Hide resolved

include/deal.II/matrix_free/tools.h Outdated Show resolved Hide resolved

kronbichler reviewed Jun 5, 2024

View reviewed changes

peterrum force-pushed the mf_tools_new branch from eee87f3 to 9652c5f Compare June 7, 2024 10:27

kronbichler approved these changes Jun 7, 2024

View reviewed changes

Generalize MatrixFreeTools::compute_diagonal() and MatrixFreeTools::c…

94027f4

…ompute_matrix()

peterrum force-pushed the mf_tools_new branch from 9652c5f to 94027f4 Compare June 7, 2024 12:26

kronbichler merged commit 57738c4 into dealii:master Jun 7, 2024
16 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

`MatrixFreeTools::compute_diagonal()` and `compute_matrix()` for face and boundary integrals #17029

`MatrixFreeTools::compute_diagonal()` and `compute_matrix()` for face and boundary integrals #17029

peterrum commented May 16, 2024

bangerth commented May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum May 16, 2024

peterrum commented May 18, 2024

kronbichler left a comment

kronbichler Jun 5, 2024

peterrum Jun 7, 2024

kronbichler Jun 7, 2024

peterrum commented Jun 7, 2024

		template <int dim, typename Number, bool is_face_>
		class ComputeMatrixScratchData

	init_data.active_fe_index =
	fe_degree != numbers::invalid_unsigned_int ?
	init_data.dof_info->fe_index_from_degree(first_selected_component,
	fe_degree) :
	(active_fe_index_given != numbers::invalid_unsigned_int ?
	active_fe_index_given :
	0);

MatrixFreeTools::compute_diagonal() and compute_matrix() for face and boundary integrals #17029

MatrixFreeTools::compute_diagonal() and compute_matrix() for face and boundary integrals #17029

Conversation

peterrum commented May 16, 2024

bangerth commented May 16, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum commented May 18, 2024

kronbichler left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum commented Jun 7, 2024

`MatrixFreeTools::compute_diagonal()` and `compute_matrix()` for face and boundary integrals #17029

`MatrixFreeTools::compute_diagonal()` and `compute_matrix()` for face and boundary integrals #17029