Optimize hanging-node constraints in MatrixFree #12560

peterrum · 2021-07-10T17:03:41Z

~~depends on #12577~~

kronbichler

I looked through the whole code and it looks fine. I have some minor comments.

include/deal.II/matrix_free/dof_info.h

include/deal.II/matrix_free/dof_info.templates.h

include/deal.II/matrix_free/evaluation_kernels.h

peterrum · 2021-07-17T14:14:47Z

~~@kronbichler I have addressed 3 out 4 of your comments. There are still two topics I need to address~~:

~~compute diagonal~~
~~storage of plain indices~~

kronbichler

Excellent - I have a few small requests regarding variable names and documentation. Then this should be in a state to give it some try.

include/deal.II/matrix_free/evaluation_kernels.h

kronbichler · 2021-07-23T10:58:24Z

include/deal.II/matrix_free/evaluation_kernels.h

+    template <int fe_degree_, unsigned int side, bool transpose>
+    static void
+    interpolate_2D(const unsigned int fe_degree,
+                   bool               type,


Can you give this a more descriptive name, e.g. flip_direction? Below, you use not_flipped, but I would also prefer the non-negated statement.

I have renamed it to is_subface_0.

include/deal.II/matrix_free/evaluation_kernels.h

kronbichler · 2021-07-23T11:53:37Z

include/deal.II/matrix_free/evaluation_template_factory.templates.h

+    Assert(false, ExcNotImplemented());
+
+    (void)n_components;
+    (void)fe_degree;
+    (void)fe_eval;
+    (void)transpose;
+    (void)c_mask;
+    (void)values;


Does this mean that face integrals computed on FE_Q with hanging nodes (e.g. for inhomogeneous Neumann conditions or weak Dirichlet conditions) end up here? But we should be able to use the same code here, shouldn't we? While there could be a case where FE_Q goes into the contiguous unknown section, I think we should be able to rule that part out. It seems the matrix_vector_faces_xx tests miss the part of FE_Q combined with hanging nodes. We do not have to do it in this PR, but it should happen soon. I can do that next week.

Does this mean that face integrals computed on FE_Q with hanging nodes (e.g. for inhomogeneous Neumann conditions or weak Dirichlet conditions) end up here?

At the moment, this function is not called for faces, since

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

Lines 5307 to 5312 in 85f67b5

FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::

apply_hanging_node_constraints() const

{

if (is_face || this->dof_info == nullptr ||

this->dof_info->hanging_node_constraint_masks.size() == 0)

return; // nothing to do with faces

It is more a way to please the compiler.

But we should be able to use the same code here, shouldn't we? While there could be a case where FE_Q goes into the contiguous unknown section, I think we should be able to rule that part out. It seems the matrix_vector_faces_xx tests miss the part of FE_Q combined with hanging nodes. We do not have to do it in this PR, but it should happen soon. I can do that next week.

I guess the code might work if the face dofs are embedded into the cell dofs. Of course, the mask would need to be updated.

include/deal.II/matrix_free/matrix_free.templates.h

include/deal.II/matrix_free/tools.h

tests/matrix_free/hanging_node_kernels_01.cc

peterrum · 2021-07-24T11:40:28Z

Excellent - I have a few small requests regarding variable names and documentation. Then this should be in a state to give it some try.

@kronbichler I'll address the comments later today. There are two open question worth discussing:

how to allocate the memory for the temporal storage for the 1D interpolations

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

Line 4616 in 85f67b5

typename Number::value_type temp[10 /*TODO*/];
When to set up the optimized hanging-node data structures? I am not really in favor of introducing a new flag in AdditionalData: as it is now the new code is well tested (since it is chosen always). I guess one could alternatively create internally a AffineConstraint object with make_hanging_node_constraints() and check compatibility with the user provided one. But I am not sure the check is always possible!?

peterrum · 2021-07-24T20:34:39Z

@kronbichler I have done the requested changes.

kronbichler · 2021-07-25T15:13:30Z

how to allocate the memory for the temporal storage for the 1D interpolations

Why not choose the length 40 or so, and then assert that we don't exceed it? We won't exceed that length in any realistic scenario, and it is still less than 64 where the compiler might start to touch every page. Since all types are basic types without constructor, we don't pay for the additional elements.

kronbichler

I think this is good to merge now. I have some comments on the generated machine code below.

include/deal.II/matrix_free/evaluation_kernels.h

peterrum · 2021-07-25T19:58:11Z

@kronbichler I have made the changes. Once you give the final go, I will squash the commits!

peterrum · 2021-07-25T20:01:53Z

/rebuild

peterrum · 2021-07-26T20:15:19Z

@kronbichler May I ask you to take look at the last commit:

Does this mean that face integrals computed on FE_Q with hanging nodes (e.g. for inhomogeneous Neumann conditions or weak Dirichlet conditions) end up here? But we should be able to use the same code here, shouldn't we? While there could be a case where FE_Q goes into the contiguous unknown section, I think we should be able to rule that part out. It seems the matrix_vector_faces_xx tests miss the part of FE_Q combined with hanging nodes. We do not have to do it in this PR, but it should happen soon. I can do that next week.

I had to implement this in this PR. Even though this feature is not tested here, it is used in the MeltPoolDG project.

I furthermore have run successfully all tests in the MeltPoolDG project, which extensively uses AMR and is a good indication that the new implementation behaves as expected.

FYI @mschreter

kronbichler

The second commit looks exactly as I thought it should, 👍

peterrum · 2021-07-27T18:24:49Z

include/deal.II/matrix_free/hanging_nodes_internal.h

+      default_value = 0,
+


@kronbichler I have introduced here a default value. Else I got some warnings in the tests in structures like this:

(direction == 0) ? dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y : 0;

Is the name fine or do you prefer another name?

I guess this does not necessarily mean that we have no constraints at all (in which case I would have chosen unconstrained), but only that we do not add new content through a binary 'or' operation, right? In that case, I am fine with the name.

I must admit I like unconstrained, since this is the way we would use it in most cases (and we could replace all the comparisons with 0). The concrete instance that caused trouble looks like this:

const bool constrained_face = (constraint_mask[v] & (((direction == 0) ? dealii::internal::MatrixFreeFunctions:: ConstraintTypes::face_y : dealii::internal::MatrixFreeFunctions:: ConstraintTypes::default_value) | ((direction == 1) ? dealii::internal::MatrixFreeFunctions:: ConstraintTypes::face_x : dealii::internal::MatrixFreeFunctions:: ConstraintTypes::default_value)));

(i.e., we add new content through a binary 'or' operation)

But then we go with unconstrained?

Yup. I have renamed it!

kronbichler

Can you rebase this (to be sure regarding #11108), then we should merge.

peterrum · 2021-07-28T09:33:10Z

Can you rebase this

Done!

jodlbauer · 2021-08-02T11:49:26Z

include/deal.II/matrix_free/tools.h

+                            for (unsigned int j = 0; j < dofs_per_component;
+                                 ++j)
+                              if (1e-10 <
+                                  std::abs(values_dofs[j][0] &&


I think there is something wrong with the parentheses here.

I just noticed the same thing and put up a patch in #12619.

This is very weird. The OSX testers both fail on other branches but they both got green checkmarks here.

I only realised because I accidentally built the examples and clang raised a warning (-Wabsolute-value).

peterrum added the Matrix-free label Jul 10, 2021

peterrum mentioned this pull request Jul 15, 2021

Rename cuda_hanging_nodes_internal.h to hanging_nodes_internal.h #12571

Merged

kronbichler reviewed Jul 16, 2021

View reviewed changes

peterrum added the ready for review label Jul 19, 2021

peterrum mentioned this pull request Jul 19, 2021

DoFInfo::read_dof_indices(): add additional index vector #12581

Merged

peterrum force-pushed the mf_hn branch 2 times, most recently from 78feffe to 6b4cad4 Compare July 20, 2021 05:34

peterrum mentioned this pull request Jul 20, 2021

Add new compute_diagonal() test #12585

Merged

peterrum force-pushed the mf_hn branch from 1602931 to 85f67b5 Compare July 20, 2021 19:19

peterrum changed the title ~~[WIP] Optimize hanging-node constraints in MatrixFree~~ Optimize hanging-node constraints in MatrixFree Jul 20, 2021

kronbichler reviewed Jul 23, 2021

View reviewed changes

peterrum force-pushed the mf_hn branch from 85f67b5 to 984e419 Compare July 24, 2021 20:33

kronbichler approved these changes Jul 25, 2021

View reviewed changes

include/deal.II/matrix_free/evaluation_kernels.h Show resolved Hide resolved

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

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

peterrum added the ready to test label Jul 25, 2021

kronbichler approved these changes Jul 25, 2021

View reviewed changes

peterrum force-pushed the mf_hn branch from 972d1b6 to 29f7601 Compare July 26, 2021 05:32

peterrum added the Do not merge ☠️ label Jul 26, 2021

kronbichler approved these changes Jul 26, 2021

View reviewed changes

peterrum removed the Do not merge ☠️ label Jul 26, 2021

peterrum force-pushed the mf_hn branch 3 times, most recently from a0bfb9a to 62d421b Compare July 27, 2021 18:21

peterrum commented Jul 27, 2021

View reviewed changes

peterrum force-pushed the mf_hn branch from 62d421b to 9f585c0 Compare July 27, 2021 18:54

kronbichler approved these changes Jul 28, 2021

View reviewed changes

peterrum force-pushed the mf_hn branch from 9f585c0 to 91cd3d0 Compare July 28, 2021 09:31

peterrum force-pushed the mf_hn branch from 91cd3d0 to dc79555 Compare July 28, 2021 09:45

Generalize HangingNodes

17f877e

peterrum force-pushed the mf_hn branch from dc79555 to 17f877e Compare July 28, 2021 18:26

kronbichler merged commit 8c776e1 into dealii:master Jul 30, 2021

jodlbauer reviewed Aug 2, 2021

View reviewed changes

kronbichler mentioned this pull request Apr 27, 2023

MatrixFree diagonal: Avoid repeated computations and binary search #15147

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Optimize hanging-node constraints in MatrixFree #12560

Optimize hanging-node constraints in MatrixFree #12560

peterrum commented Jul 10, 2021 •

edited

kronbichler left a comment

peterrum commented Jul 17, 2021 •

edited

kronbichler left a comment

kronbichler Jul 23, 2021

peterrum Jul 24, 2021

kronbichler Jul 23, 2021

peterrum Jul 24, 2021

peterrum commented Jul 24, 2021

peterrum commented Jul 24, 2021

kronbichler commented Jul 25, 2021 •

edited

kronbichler left a comment

peterrum commented Jul 25, 2021

peterrum commented Jul 25, 2021

peterrum commented Jul 26, 2021

kronbichler left a comment

peterrum Jul 27, 2021

kronbichler Jul 27, 2021

peterrum Jul 27, 2021 •

edited

kronbichler Jul 27, 2021

peterrum Jul 27, 2021

kronbichler left a comment

peterrum commented Jul 28, 2021

jodlbauer Aug 2, 2021

drwells Aug 2, 2021

jodlbauer Aug 2, 2021

	FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
	apply_hanging_node_constraints() const
	{
	if (is_face \|\| this->dof_info == nullptr \|\|
	this->dof_info->hanging_node_constraint_masks.size() == 0)
	return; // nothing to do with faces

Optimize hanging-node constraints in MatrixFree #12560

Optimize hanging-node constraints in MatrixFree #12560

Conversation

peterrum commented Jul 10, 2021 • edited

kronbichler left a comment

Choose a reason for hiding this comment

peterrum commented Jul 17, 2021 • edited

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

Choose a reason for hiding this comment

peterrum commented Jul 24, 2021

peterrum commented Jul 24, 2021

kronbichler commented Jul 25, 2021 • edited

kronbichler left a comment

Choose a reason for hiding this comment

peterrum commented Jul 25, 2021

peterrum commented Jul 25, 2021

peterrum commented Jul 26, 2021

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

peterrum Jul 27, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kronbichler left a comment

Choose a reason for hiding this comment

peterrum commented Jul 28, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

peterrum commented Jul 10, 2021 •

edited

peterrum commented Jul 17, 2021 •

edited

kronbichler commented Jul 25, 2021 •

edited

peterrum Jul 27, 2021 •

edited