Add a test for moving particles in a complex geometry and

[blaisb]

This PR aimed at addressing #13246, but a better solution was found by @kronbichler in #13278 . However, the tolerance of the particle_handler locating algorithm is still a bit too harsh and, consequently, particles can still be lost when crossing cells. This PR fixes the tolerance and also a unit test that I created to reproduce the bug I encountered in #13426. I think it would be worthwhile to keep the test since it helped me debug the issue quite nicely.

The test is very simple. 48 particles are inserted in an hyper_ball. They are moved slightly using Euler's method. Before #13278 a particle would be lost on the second step. I think we need more tests for particles ( It has been a few times that I found issues in particles using Lethe :) ).

[kronbichler]

This must happen in

dealii/source/fe/mapping_q.cc

Lines 793 to 818 in 5478749

	Point<dim, VectorizedArray<double>> unit_point =
	internal::MappingQImplementation::
	do_transform_real_to_unit_cell_internal<dim, spacedim>(
	p_vec,
	inverse_approximation.compute(p_vec),
	support_points,
	polynomials_1d,
	renumber_lexicographic_to_hierarchic);
	// If the vectorized computation failed, it could be that only some of
	// the lanes failed but others would have succeeded if we had let them
	// compute alone without interference (like negative Jacobian
	// determinants) from other SIMD lanes. Repeat the computation in this
	// unlikely case with scalar arguments.
	for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
	if (unit_point[0][j] == std::numeric_limits<double>::infinity())
	unit_points[i + j] = internal::MappingQImplementation::
	do_transform_real_to_unit_cell_internal<dim, spacedim>(
	real_points[i + j],
	inverse_approximation.compute(real_points[i + j]),
	support_points,
	polynomials_1d,
	renumber_lexicographic_to_hierarchic);
	else
	for (unsigned int d = 0; d < dim; ++d)
	unit_points[i + j][d] = unit_point[d][j];

We already repeat things if we happen to fail, so it surprises me that we do not catch that case (as you described in more detail on slack). I will take a look.

I think the current solution is a bit unfortunate because it forces a re-computation of all points that left the cell, which can be a substantial cost. So we should first try to see if there is a more fundamental thing we can fix.

[blaisb]

This must happen in

dealii/source/fe/mapping_q.cc

Lines 793 to 818 in 5478749

	Point<dim, VectorizedArray<double>> unit_point =
	internal::MappingQImplementation::
	do_transform_real_to_unit_cell_internal<dim, spacedim>(
	p_vec,
	inverse_approximation.compute(p_vec),
	support_points,
	polynomials_1d,
	renumber_lexicographic_to_hierarchic);
	// If the vectorized computation failed, it could be that only some of
	// the lanes failed but others would have succeeded if we had let them
	// compute alone without interference (like negative Jacobian
	// determinants) from other SIMD lanes. Repeat the computation in this
	// unlikely case with scalar arguments.
	for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
	if (unit_point[0][j] == std::numeric_limits<double>::infinity())
	unit_points[i + j] = internal::MappingQImplementation::
	do_transform_real_to_unit_cell_internal<dim, spacedim>(
	real_points[i + j],
	inverse_approximation.compute(real_points[i + j]),
	support_points,
	polynomials_1d,
	renumber_lexicographic_to_hierarchic);
	else
	for (unsigned int d = 0; d < dim; ++d)
	unit_points[i + j][d] = unit_point[d][j];

We already repeat things if we happen to fail, so it surprises me that we do not catch that case (as you described in more detail on slack). I will take a look.
I think the current solution is a bit unfortunate because it forces a re-computation of all points that left the cell, which can be a substantial cost. So we should first try to see if there is a more fundamental thing we can fix.

@kronbichler It would be really nice to find a better solution than the one I propose. I do not like my solution at all, because, like you said, it duplicate the costs.
Maybe this identified some sort of remaining issue in transform_points_real_to_unit_cell that we can fix and then we can only keep the test of this PR and not the rest :). That would make me very happy :).

[kronbichler]

@blais I now made a correction in #13278 - I believe this obviates the current patch, requiring only the test case (which is indeed valuable). I added another test in that PR that picks only the point that did not succeed: The algorithm would find another inverse point far away from the unit cell that also happened to have zero residual. We came there by bad luck, coming from a very bad initial guess.

[blaisb]

@kronbichler you're the best :).
Do you want me to remove the changes to particle_handler in this PR and to just keep the test in this one so this PR just becomes the addition of the test?

[kronbichler]

Yes, @blaisb , that would be great to have the test.

[blaisb]

@gassmoeller It seems the tolerance issue was also there on top of the bug addressed by @kronbichler . I think it's good if you review this one.

[tjhei]

particle_handler.cc:1379:44: error: implicit conversion from 'double' to 'bool' changes non-zero value from 1.0E-12 to true [-Werror,-Wfloat-zero-conversion]

[bangerth]

You've got the parentheses wrong here. You write

  if (foo(), 1e-12)

where you mean

  if (foo(1e-12))

Nice catch by our warnings infrastructure :-)

[blaisb]

nice catch! My bad, I should have seen it from the automatic indentation.
I'm not sure if 1e-12 is the best threshold here. It seems to work in most cases I tried.

[gassmoeller]

If you fix the parenthesis this looks good to me. I would suggest a threshold of 1e-10 as this is the default for GridTools::find_active_cell_around_point, that way #13202 combined with this PR should change as little as possible.

[gassmoeller]

The PR looks good to me, thanks for narrowing down the problem (and thanks to @kronbichler for fixing the issue inside the mapping!). See my comment, 1e-10 should be the best threshold if only for compatibility.

[gassmoeller]

If you fix the parenthesis this looks good to me. I would suggest a threshold of 1e-10 as this is the default for GridTools::find_active_cell_around_point, that way #13202 combined with this PR should change as little as possible.

[blaisb]

Dear all,
I reverted the changes on the particle_handler. Changing the tolerance did not fix the bug I am encountering in some very narrow edge case (we have one falling test remaining in Lethe where a particle disappears under very small motion).
I suggest we move on and merge this PR with just the test. I will try to fix the remaining issue with the disappearing particles and make a subsequent PR + another test. @gassmoeller maybe we can try to have a chat about this this week?

My conclusions is that I will need to write additional test for particles :).

[bangerth]

On 1/24/22 6:19 PM, Bruno Blais wrote: My conclusions is that I will need to write additional test for particles :).

More tests = always more better.

[bangerth]

/rebuild

[blaisb]

@bangerth I think this might be ready for merge since both @kronbichler and @gassmoeller approved :).