Bisection bugs in dipolar P3M

[jngrad]

The dipolar P3M algorithm fails for specific particle configurations. The bug depends on the box size and operating system. MWE:

import espressomd
import espressomd.magnetostatics
BOX_L = 1
system = espressomd.System(box_l=3*[BOX_L])
system.time_step = 0.01
system.part.add(pos=[[0, 0, 0], [0.5, 0.5, 0.5]], rotation=2 * [(1, 1, 1)], dip=2 * [(1, 0, 0)])
solver = espressomd.magnetostatics.DipolarP3M(prefactor=1, accuracy=1e-2)
system.actors.add(solver)

• with box size 1: ./pypresso mwe.py triggers an assertion Root must be bracketed for bisection in dp3m_rtbisection in function double dp3m_rtbisection(...), initially reported in https://github.com/espressomd/espresso/pull/3869#issuecomment-676686275
  • this assertion was introduced to prevent an infinite loop
• with box size of 10: mpiexec -n 4 ./pypresso mwe.py triggers an assertion and generates a stack trace (mwe.log) randomly, depends on the operating system, reported in https://github.com/espressomd/docker/pull/189#issuecomment-699689428
  • with dip = 2 * [(0.1, 0, 0)], it instead triggers the bisection assertion

This issue makes the test_09_no_errors_dp3m_cpu() check in testsuite/python/p3m_tuning_exceptions.py fragile in CI, which prevents us from updating the ROCm image to Ubuntu 20.04.

[RudolfWeeber]

Assigned #3915 https://github.com/espressomd/espresso/issues/3915 to @RudolfWeeber https://github.com/RudolfWeeber .

I’m not sure, what to do about this. I do not understand the dipolar p3m code, unfortunately.

[jngrad]

It now fails on ICP computer sheep at random. This could be a problem during the summer school.

[RudolfWeeber]

My understanding is the following:

• The error occurs in the bisection which finds the best alpha for a given real-space cutoff.
• The error message means "You can't bisect if the evaluated function is not <0 on the one end and >0 on the other end of the bisection interval."
• In the MWE, this is the case, because the real space error for the cutoff which comes out of the previous tuning steps is never small enough.
• The maximum real space cutoff is limited by the size of the (mpi-local) box + skin Otherwise pairs would be missed in the short range loop for the real space part.

I assume that the previous tuning steps, which determins the number of mesh points, assignment order and r_cut doesn't handle the case where r_cut becomes bigger than half the local box properly. Unfortunately, this code is really hard to understand. Rather than digging further, we should probably rewrite this in a more functional style and clearly write down the pre- and postconditions for each step. While we are at it, we can introduce the possibility to inject artificial timings, so the entire thing can be tested. And include skin tuning for that matter.

@pkreissl, @pstaerk, where/under what condition did the issue occur in the tutorial?