Hessians on general cells and faces #12744
Conversation
bergbauer
force-pushed
the
hessians_on_general_cells
branch
from
975efae
to
b660559
Compare
|
/rebuild |
bergbauer
force-pushed
the
hessians_on_general_cells
branch
from
b660559
to
d6b5fab
Compare
|
I applied your suggestions and fixed the failing tests. |
kronbichler
changed the title
There was a problem hiding this comment.
This is really a substantial improvement so I think it is well-placed there. Thanks a lot for the effort @bergbauer !
I agree with @kronbichler here! I am looking forward to try it out in my small matrix-free biharmonic solver 🥇
| hessians_quad = | ||
| gradients_from_hessians_quad = | ||
| this->scratch_data_array->begin() + | ||
| n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); | ||
| this->scratch_data = | ||
| this->scratch_data_array->begin() + n_components_ * dofs_per_component + | ||
| (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points); | ||
| hessians_quad = | ||
| this->scratch_data_array->begin() + | ||
| n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points); | ||
| this->scratch_data = this->scratch_data_array->begin() + | ||
| n_components_ * dofs_per_component + | ||
| (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) * | ||
| n_quadrature_points); |
There was a problem hiding this comment.
Not sure, but don't we only want the memory if MF was actually set up with Hessian enabled?
There was a problem hiding this comment.
This is what we did before. We maybe waste 4 or 5 instructions here to get the intermediate pointers, against an if and some bookkeeping to check if we had hessians filled in MappingInfo, so I think the code is cleaner this way. The local arrays are small enough to not really eat much into the memory of a global finite element computation, and in case we have no hessians we will never touch them (and hence not waste cache memory).
There was a problem hiding this comment.
We could move it to the end of scratch_data_array and disable it in case the flag UpdateFlags::update_jacobian_grads is not set in mapping_data->update_flags_inner_faces. This can also be done in a follow up if we think it is necessary.
| const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point]; | ||
| const VectorizedArrayType JxW = this->J_value[q_point]; |
There was a problem hiding this comment.
| const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point]; | |
| const VectorizedArrayType JxW = this->J_value[q_point]; | |
| const Tensor<2, dim, VectorizedArrayType>& jac = this->jacobian[q_point]; | |
| const VectorizedArrayType& JxW = this->J_value[q_point]; |
There was a problem hiding this comment.
We are a bit inconsistent here, as we do
dealii/include/deal.II/matrix_free/fe_evaluation.h
Lines 5805 to 5808 in b6bd942
dealii/include/deal.II/matrix_free/fe_evaluation.h
Lines 6279 to 6282 in b6bd942
In the end, it will not matter much because a sane compiler will likely treat both cases similarly and avoid temporaries as much as possible. I am not sure if aliasing could force the compiler to emit unnecessary loads in the case of a reference and if that was the reason to do a copy in the second appearance for
submit_gradient. I suggest we leave this as is (no reference) and just do something in case we have verified the reference does not lead to unnecessary loads in the multi-component case in the assembly code.
There was a problem hiding this comment.
You could show me how to look for this in assembly code.
| @@ -8281,12 +8348,19 @@ FEEvaluation<dim, | |||
| Number, | |||
| VectorizedArrayType>:: | |||
| evaluate(const VectorizedArrayType * values_array, | |||
| const EvaluationFlags::EvaluationFlags evaluation_flags) | |||
| const EvaluationFlags::EvaluationFlags evaluation_flag) | |||
There was a problem hiding this comment.
Don't we normally write flags (plural)?
There was a problem hiding this comment.
This was the only place where plural was used. I changed it to singular to be consistent.
bergbauer
force-pushed
the
hessians_on_general_cells
branch
from
d6b5fab
to
92774b2
Compare
This PR implements the evaluation and integration of hessians for general (non-affine) cells/faces in the matrix-free infrastructure.
Closes #10401