-
Notifications
You must be signed in to change notification settings - Fork 35
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement additive Schwarz Preconditioner #589
Conversation
// store the cell matrix and renumber lexicographic | ||
auto & lapack_matrix = matrices[cell_batch_index * vectorization_length + v]; | ||
|
||
auto const & lex = matrix_free->get_shape_info().lexicographic_numbering; | ||
for(unsigned int i = 0; i < dofs_per_cell; i++) | ||
for(unsigned int j = 0; j < dofs_per_cell; j++) | ||
lapack_matrix.set(i, j, overlapped_cell_matrix[lex[i]][lex[j]]); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the only part of the code that is not fully clear in my opinion. In particular, I do not understand "renumber lexicographic".
It seems that the numbering needs to be consistent with read_dof_values()
/begin_dof_values()
used in the function apply_additive_schwarz_matrices()
. Is there a way to make this more clear? Why isn't there a function like read_dof_indices()
in analogy to read/begin_dof_values()
? Personally, I am against using auto
in such cases, but there might be different opinions.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe formulated differently: What is the reason that the code written here is / can not written as a matrix_free->celll_loop()
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
FE_Q
is numbered hierarchical but FEEvaluation
expects lexicographic numbering. Should we add the same reference like here
exadg/include/exadg/operators/operator_base.cpp
Lines 1894 to 1902 in ee073c7
if(not is_dg) | |
{ | |
// in the case of CG: shape functions are not ordered lexicographically | |
// see (https://www.dealii.org/8.5.1/doxygen/deal.II/classFE__Q.html) | |
// so we have to fix the order | |
auto temp = dof_indices; | |
for(unsigned int j = 0; j < dof_indices.size(); j++) | |
dof_indices[j] = temp[matrix_free.get_shape_info().lexicographic_numbering[j]]; | |
} |
FEEvaluation::get_dof_indices()
would indeed sometimes make sense, but we usually don't access the dof indices in performance-critical places.
We could also write this as a cell loop but I tend to write out the for loop over cell batches when we don't work on dof vectors.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I am somewhat confused. according to your explanations and https://dealii.org/current/doxygen/deal.II/structinternal_1_1MatrixFreeFunctions_1_1ShapeInfo.html#abad0471d12c4f07f3d83271377e24231, I would expect the code to read
lapack_matrix.set(lex[i], lex[j], overlapped_cell_matrix[i][j]);
because lapack_matrix
uses lexicographic numerbing and overlapped_cell_matrix
hierarchical numbering, and because shape_info::lexicographic_numbering is a "Renumbering from deal.II's numbering of cell degrees of freedom to lexicographic numbering used inside FEEvaluation"?
Generally, I think one should be able to get dof indices in the meaning of / consistent with integrator.read/begin_dof_values()
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I discussed this briefly with @peterrum. We think both the name shape_info.lexicograhpic_numbering
and the deal.II docu is confusing. It should be lexicographic_to_hierarchical_numbering
in my opinion. The current implementation of this PR would then be correct apart from renaming.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The documentation is really confusing, I agree. The renumbering in correct as it currently is.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I renamed the variable but the docu has to be fixed in deal.II. 13e5d64
The code looks overall good from my perspective. I raised some questions, mainly to improve my own understanding of the code (or |
src.update_ghost_values(); | ||
|
||
matrix_free->template cell_loop<VectorType, VectorType>( | ||
[&](auto const & matrix_free, auto & dst, auto const & src, auto const & cell_range) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Personally, I am a fan to keep code with lambdas as simple as possible, which is why I would probably define the lambda outside (e.g. cellwise_operation
). But I am fine with the present version as well.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
looks good!
This PR implements a matrix-based additive Schwarz preconditioner where a sparse matrix is assembled first, then cell blocks are cut out and factorized using LU factorization. It can then be used as a preconditioner or as a smoother in multigrid.