Skip to content
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.

Sign up for GitHub

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

Jump to bottom

NonMatching::MappingInfo: precompute mapping data for an IteratorRange of cells/faces #14809

Merged
merged 1 commit into from
Mar 14, 2023
Merged

NonMatching::MappingInfo: precompute mapping data for an IteratorRange of cells/faces #14809

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
5 changes: 4 additions & 1 deletion include/deal.II/fe/mapping.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,9 @@ namespace NonMatching
{
template <int dim>
class FEImmersedSurfaceValues;
}
template <int dim, int spacedim>
class MappingInfo;
} // namespace NonMatching


/**
Expand Down Expand Up @@ -1305,6 +1307,7 @@ class Mapping : public Subscriptor
friend class FEFaceValues<dim, spacedim>;
friend class FESubfaceValues<dim, spacedim>;
friend class NonMatching::FEImmersedSurfaceValues<dim>;
friend class NonMatching::MappingInfo<dim, spacedim>;
};


Expand Down
172 changes: 144 additions & 28 deletions include/deal.II/matrix_free/fe_point_evaluation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -470,6 +470,20 @@ class FEPointEvaluation
reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const ArrayView<const Point<dim>> &unit_points);

/**
* Reinitialize the evaluator to point to the correct precomputed mapping of
* the cell in the MappingInfo object.
*/
void
reinit(const unsigned int cell_index);

/**
* Reinitialize the evaluator to point to the correct precomputed mapping of
* the face in the MappingInfo object.
*/
void
reinit(const unsigned int cell_index, const unsigned int face_number);

/**
* This function interpolates the finite element solution, represented by
* `solution_values`, on the cell and `unit_points` passed to reinit().
Expand Down Expand Up @@ -585,6 +599,22 @@ class FEPointEvaluation
DerivativeForm<1, spacedim, dim>
inverse_jacobian(const unsigned int point_index) const;

/**
* Return the Jacobian determinant multiplied by the quadrature weight. This
* class or the MappingInfo object passed to this function needs to be
* constructed with UpdateFlags containing `update_JxW_values`.
*/
Number
JxW(const unsigned int point_index) const;

/**
* Return the normal vector. This class or the MappingInfo object passed to
* this function needs to be constructed with UpdateFlags containing
* `update_normal_vectors`.
*/
Tensor<1, spacedim>
normal_vector(const unsigned int point_index) const;

/**
* Return the position in real coordinates of the given point index among
* the points passed to reinit().
Expand All @@ -599,6 +629,11 @@ class FEPointEvaluation
Point<dim>
unit_point(const unsigned int point_index) const;

/**
* Number of quadrature points of the current cell/face.
*/
const unsigned int n_q_points;

private:
/**
* Common setup function for both constructors. Does the setup for both fast
Expand Down Expand Up @@ -713,10 +748,20 @@ class FEPointEvaluation
*/
SmartPointer<NonMatching::MappingInfo<dim, spacedim>> mapping_info;

/**
* The current cell index to access mapping data from mapping info.
*/
unsigned int current_cell_index;

/**
* The current face number to access mapping data from mapping info.
*/
unsigned int current_face_number;

/**
* The reference points specified at reinit().
*/
std::vector<Point<dim>> unit_points;
ArrayView<const Point<dim>> unit_points;
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed this here to ArrayView to avoid the copy of unit points. I can also address this in a follow-up, when we get rid of this class member as proposed in https://github.com/dealii/dealii/pull/14809/files#r1125725149 .


/**
* Bool indicating if fast path is chosen.
Expand All @@ -733,13 +778,16 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
const FiniteElement<dim> &fe,
const UpdateFlags update_flags,
const unsigned int first_selected_component)
: mapping(&mapping)
: n_q_points(numbers::invalid_unsigned_int)
, mapping(&mapping)
, fe(&fe)
, update_flags(update_flags)
, mapping_info_on_the_fly(
std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
update_flags))
, mapping_info(mapping_info_on_the_fly.get())
, current_cell_index(numbers::invalid_unsigned_int)
, current_face_number(numbers::invalid_unsigned_int)
{
setup(first_selected_component);
}
Expand All @@ -751,10 +799,13 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
NonMatching::MappingInfo<dim, spacedim> &mapping_info,
const FiniteElement<dim> & fe,
const unsigned int first_selected_component)
: mapping(&mapping_info.get_mapping())
: n_q_points(numbers::invalid_unsigned_int)
, mapping(&mapping_info.get_mapping())
, fe(&fe)
, update_flags(mapping_info.get_update_flags())
, mapping_info(&mapping_info)
, current_cell_index(numbers::invalid_unsigned_int)
, current_face_number(numbers::invalid_unsigned_int)
{
setup(first_selected_component);
}
Expand Down Expand Up @@ -853,14 +904,45 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
fe_values->reinit(cell);
}

this->unit_points =
std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
this->unit_points = unit_points;

const_cast<unsigned int &>(n_q_points) = unit_points.size();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see why you do this. I do the same in FEFaceValues to support mixed meshes. The reason for me to do that there was because n_q_points was already a const and public member variable. I was not happy with the solution there, since it is very error prone on user side if one should decide to store n_q_points in a variable. I would much more prefer if we could either introduce a getter function or, if we just need this to iterate over quadrature points, a function like quadrature_point_indices() in FEEvaluation (https://www.dealii.org/developer/doxygen/deal.II/classFEEvaluationData.html#a0eca800437c07c9f23ec6815e58f1432).


if (update_flags & update_values)
values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
values.resize(n_q_points, numbers::signaling_nan<value_type>());
if (update_flags & update_gradients)
gradients.resize(unit_points.size(),
numbers::signaling_nan<gradient_type>());
gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
}



template <int n_components, int dim, int spacedim, typename Number>
void
FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
const unsigned int cell_index)
{
current_cell_index = cell_index;
current_face_number = numbers::invalid_unsigned_int;

const_cast<unsigned int &>(n_q_points) =
mapping_info->get_unit_points(current_cell_index, current_face_number)
.size();
}



template <int n_components, int dim, int spacedim, typename Number>
void
FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
const unsigned int cell_index,
const unsigned int face_number)
{
current_cell_index = cell_index;
current_face_number = face_number;

const_cast<unsigned int &>(n_q_points) =
mapping_info->get_unit_points(current_cell_index, current_face_number)
.size();
}


Expand All @@ -874,7 +956,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
if (precomputed_mapping)
{
unit_points = mapping_info->get_unit_points();
unit_points =
mapping_info->get_unit_points(current_cell_index, current_face_number);
Comment on lines +959 to +960
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this allocate memory in the operator= in case the old size of unit_points is already correct? Probably not, but it would be good to check.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do I check this? We could get rid of unit_points in FEPointEvaluation entirely and always ask MappingInfo (maybe in a follow-up).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we could move all code into MappingInfo (also the case where the user never sets up such an object explicitly, e.g. in particle codes that want to evaluate information once), that would be great. Regarding how to check it: I use a compiler/trace tool and monitor the allocations, first using toy codes, to see how the compiler reacts. Here I use this one:

#include <vector>
#include <iostream>

int main()
{
  std::vector<unsigned int> data(5, 3);
  std::vector<unsigned int> second(3, 6);
  for (unsigned int i=0; i<10; ++i)
    data = second;
  std::cout << data.size() << " " << data.capacity() << " " << second.size() << " " << second.capacity() << std::endl;
}

Using valgrind as valgrind ./a.out to check for memory errors, I now see

==269622==   total heap usage: 4 allocs, 4 frees, 73,760 bytes allocated

This is less than 10 repetitions, so it can't allocate memory inside the loop (one might wonder if compilers might optimize it away, as it does not change the result, but at least gcc does not here). The reason for the loop is that I see more than 2 allocs, which I would expect for the two vectors - gcc counts the full program here. The cause is that std::cout << does another one to put the data to the output stream, and the last one is in the global program init outside main (it is done inside call_init.part.0 and part of the program startup, outside the control of my program).
The output is also revealing, which is why I included it:

3 5 3 3

You see from the first two numbers that the vector data kept its original capacity of 5, even though the length was changed to 3 after assignment from second. This is a strong indication that it kept its memory region. You could output the address of &data[0] to ensure it kept the same. You can also look at the source code of stl_vector.h to see this (for me, the file is /usr/include/c++/12/bits/stl_vector.h, but the implementation is in )/usr/include/c++/12/bits/vector.tcc:

  template<typename _Tp, typename _Alloc>
    _GLIBCXX20_CONSTEXPR
    vector<_Tp, _Alloc>&
    vector<_Tp, _Alloc>::
    operator=(const vector<_Tp, _Alloc>& __x)
    {
       ...
          const size_type __xlen = __x.size();
       ...
          else if (size() >= __xlen)
            {
              std::_Destroy(std::copy(__x.begin(), __x.end(), begin()),
                            end(), _M_get_Tp_allocator());
            }
       ...

From this, you see that we only call std::copy (and destroy the elements that now got invalid due to reducing the size), but we did not deallocate/reallocate.


if (update_flags & update_values)
values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
Expand All @@ -883,7 +966,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
numbers::signaling_nan<gradient_type>());
}

if (unit_points.empty())
if (unit_points.size() == 0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's keep this show it was.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with @peterrum - just an interesting remark, compiler people have been working on detecting size() == 0 and turn it into something like empty() for performance reasons. To give insights on why there might be a difference (apart from code reading arguments, which are also there): The latter checks where the addresses of begin() and end() are the same (one comparison, plus possibly loads on begin() and end()), whereas size() first computes temp = end() - begin(), and then we compare temp == 0, so two instructions instead of one.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed this because ArrayView does not provide an empty() function. If we decide to keep the local vector copy I will revert this change.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about adding empty() to ArrayView. We can do that in a follow-up PR!

return;

AssertDimension(solution_values.size(), fe->dofs_per_cell);
Expand Down Expand Up @@ -946,11 +1029,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
Number>::set_gradient(val_and_grad.second,
j,
unit_gradients[i + j]);
gradients[i + j] =
apply_transformation(mapping_info->get_mapping_data()
.inverse_jacobians[i + j]
.transpose(),
unit_gradients[i + j]);
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index,
current_face_number);
gradients[i + j] = apply_transformation(
mapping_data.inverse_jacobians[i + j].transpose(),
unit_gradients[i + j]);
}
}
}
Expand Down Expand Up @@ -1022,7 +1106,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
if (precomputed_mapping)
{
unit_points = mapping_info->get_unit_points();
unit_points =
mapping_info->get_unit_points(current_cell_index, current_face_number);

if (update_flags & update_values)
values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
Expand Down Expand Up @@ -1086,9 +1171,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
if (integration_flags & EvaluationFlags::gradients)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
{
gradients[i + j] = apply_transformation(
mapping_info->get_mapping_data().inverse_jacobians[i + j],
gradients[i + j]);
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index,
current_face_number);
gradients[i + j] =
apply_transformation(mapping_data.inverse_jacobians[i + j],
gradients[i + j]);
internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
gradient, j, gradients[i + j]);
Expand Down Expand Up @@ -1250,9 +1338,10 @@ inline DerivativeForm<1, dim, spacedim>
FEPointEvaluation<n_components, dim, spacedim, Number>::jacobian(
const unsigned int point_index) const
{
AssertIndexRange(point_index,
mapping_info->get_mapping_data().jacobians.size());
return mapping_info->get_mapping_data().jacobians[point_index];
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index, current_face_number);
AssertIndexRange(point_index, mapping_data.jacobians.size());
return mapping_data.jacobians[point_index];
}


Expand All @@ -1262,9 +1351,35 @@ inline DerivativeForm<1, spacedim, dim>
FEPointEvaluation<n_components, dim, spacedim, Number>::inverse_jacobian(
const unsigned int point_index) const
{
AssertIndexRange(point_index,
mapping_info->get_mapping_data().inverse_jacobians.size());
return mapping_info->get_mapping_data().inverse_jacobians[point_index];
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index, current_face_number);
AssertIndexRange(point_index, mapping_data.inverse_jacobians.size());
return mapping_data.inverse_jacobians[point_index];
}



template <int n_components, int dim, int spacedim, typename Number>
inline Number
FEPointEvaluation<n_components, dim, spacedim, Number>::JxW(
const unsigned int point_index) const
{
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index, current_face_number);
AssertIndexRange(point_index, mapping_data.JxW_values.size());
return mapping_data.JxW_values[point_index];
}


template <int n_components, int dim, int spacedim, typename Number>
inline Tensor<1, spacedim>
FEPointEvaluation<n_components, dim, spacedim, Number>::normal_vector(
const unsigned int point_index) const
{
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index, current_face_number);
AssertIndexRange(point_index, mapping_data.normal_vectors.size());
return mapping_data.normal_vectors[point_index];
}


Expand All @@ -1274,9 +1389,10 @@ inline Point<spacedim>
FEPointEvaluation<n_components, dim, spacedim, Number>::real_point(
const unsigned int point_index) const
{
AssertIndexRange(point_index,
mapping_info->get_mapping_data().quadrature_points.size());
return mapping_info->get_mapping_data().quadrature_points[point_index];
const auto &mapping_data =
mapping_info->get_mapping_data(current_cell_index, current_face_number);
AssertIndexRange(point_index, mapping_data.quadrature_points.size());
return mapping_data.quadrature_points[point_index];
}


Expand Down