Problem with periodic condition on small physical domains

[Alchem334]

Greetings!

I found a bug in the periodic boundary condition that occurs on meshes with a small element size. Here is the test case https://gist.github.com/Alchem334/0f3d0c9b7805d4a1e5438c41c3dfee3e

In the test case, the Laplace equation is solved on a square mesh of size 2e-5 x 2e-5 with a periodic condition on the faces x = 0 and x = 2e-5. After constructing the MultiPointConstraint object, you can see the number of degrees of freedom for which the periodic boundary condition is specified. In case of (20,20) elements on the mesh, number of dofs mpc.num_local_slaves is equal 20, which is not correct, there should be 21 degrees of freedom. With mesh of (10,10) elements, the periodic condition is constructed correctly with mpc.num_local_slaves=11.

The loss of one degree of freedom in the periodic condition affects the results. In the case of the (10,10) grid, the values in the upper right and upper left corners are the same.

In the case of the grid (20,20) they are different.

The problem occurs in the compute_distance_gjk function, which calculates the distance vector from a point to a cell.

https://github.com/FEniCS/dolfinx/blob/74facfc3a587be94a8936bec0aece8b79d92f62e/cpp/dolfinx/geometry/gjk.cpp#L238-L253

In the case of a point lying on a cell, this function should return a zero vector, but on a grid (20,20) cells for one of the degrees of freedom, it returns a non-zero vector, apparently with a length of 1e-6. This is because the exit criterion specified by eps is met too early. If you reduce the value of eps for example to 1e-20 then no error occurs.

The compute_distance_gjk function is found in dolfinx and is called in dolfinx_mpc as follows. This function is used in the shortest_vector function to find the shortest vector, which is used in the squared_distance function to find the squared distance. Based on the result of the squared_distance function, a decision is made whether to look for the degree of freedom in a given cell. This happens in the compute_colliding_cells function.

dolfinx_mpc/cpp/utils.cpp

Lines 1011 to 1020 in f7809a6

	const std::vector<double> distances_sq
	= dolfinx::geometry::squared_distance(mesh, tdim, cells, _point);
	// Only push back closest cell
	if (auto cell_idx
	= std::min_element(distances_sq.cbegin(), distances_sq.cend());
	*cell_idx < eps2)
	{
	auto pos = std::distance(distances_sq.cbegin(), cell_idx);
	colliding_cells.push_back(cells[pos]);
	}

The criterion is the squared distance less than eps2 = 1e-20. Since the square of the distance is equal to 1e-12, this condition is not met and for one degree of freedom no cells will be found at all in which its pair could be located.

dolfinx_mpc/cpp/utils.cpp

Lines 1040 to 1052 in f7809a6

	// Compute exact collision
	auto cell_collisions = dolfinx_mpc::compute_colliding_cells(
	mesh, bbox_collisions, points, eps2);
	// Extract first collision
	std::vector<std::int32_t> collisions(points.shape(0), -1);
	for (int i = 0; i < cell_collisions.num_nodes(); i++)
	{
	auto local_cells = cell_collisions.links(i);
	if (!local_cells.empty())
	collisions[i] = local_cells[0];
	}
	return collisions;

Finally because of this the find_local_collisions function returns -1 for this degree of freedom instead of the cell id.

I'm not sure how to fix this error once and for all. For my purposes, I simply choose eps in the compute_distance_gjk function such that there is a corresponding pair for all periodic degrees of freedom. I'm not sure you can set eps = 1e-200 and forget about this problem. In this case, some more bugs may come out.

[jorgensd]

I would strongly recommend scaling your problem such that the cell size is not 2e-5/N this can cause all kinds of other issues, as dolfinx uses double precision, so setting eps2 smaller than about 4e-16 does not make sense.

[Alchem334]

Of course, it's always better to scale the equations, but the cell size 1e-6 is not that small. Specifically, in the compute_distance_gjk function, the value of eps is compared with the scalar product of vectors inside the cell, then even in the case of dimensionless equations and double precision, a reasonable choice of eps is not less than ~1e-32. Well, or compare not with eps, but with eps * eps, as done a dozen lines below.

https://github.com/FEniCS/dolfinx/blob/74facfc3a587be94a8936bec0aece8b79d92f62e/cpp/dolfinx/geometry/gjk.cpp#L290-L292

By the way, the eps2 variable in dolfinx_mpc itself is set to 1e-20

dolfinx_mpc/cpp/PeriodicConstraint.cpp

Lines 139 to 140 in b43e139

	std::vector<std::int32_t> local_cell_collisions
	= dolfinx_mpc::find_local_collisions(*mesh, tree, mapped_T_b, 1e-20);

dolfinx_mpc/cpp/utils.cpp

Lines 883 to 886 in b43e139

	std::vector<std::int32_t> dolfinx_mpc::find_local_collisions(
	const dolfinx::mesh::Mesh& mesh,
	const dolfinx::geometry::BoundingBoxTree& tree,
	std::span<const double> points, const double eps2)