espressomd / espresso

The ESPResSo package
https://espressomd.org
GNU General Public License v3.0
222 stars 183 forks source link

Inertialess tracers are broken #4902

Closed jngrad closed 2 months ago

jngrad commented 2 months ago

In the LB particle coupling code, real particles are folded onto the local domain to eliminate edge cases, such as when a particle is located at pos=[-1e-30,-1e-30,-1e-30] in unfolded coordinates. Here is the corresponding code:

https://github.com/espressomd/espresso/blob/4285bb685c4c5e5b1fca8a3c4abfada8c5f0319e/src/core/lb/particle_coupling.cpp#L179-L180

The LB inertialess tracers do not take this precaution:

https://github.com/espressomd/espresso/blob/4285bb685c4c5e5b1fca8a3c4abfada8c5f0319e/src/core/virtual_sites/lb_tracers.cpp#L63-L65

In addition, the inertialess tracers code converts positions from MD units to LB units, even though add_md_force() expects positions in MD units since it calls LB::Solver::add_force_density(), which converts positions from MD units to LB units:

https://github.com/espressomd/espresso/blob/4285bb685c4c5e5b1fca8a3c4abfada8c5f0319e/src/core/lb/Solver.cpp#L198-L200

When using agrid values other than unity, inertialess tracers are converted twice. The confusion is probably due to the fact that LB::Solver methods don't convert all quantities. For example, forces aren't converted, but positions and velocities are.

Here is a MWE:

import espressomd
import espressomd.lb
import espressomd.propagation

system = espressomd.System(box_l=3 * [6.0])
system.time_step = 0.01
system.cell_system.skin = 0.4
lbf = espressomd.lb.LBFluidWalberla(
    tau=0.01, agrid=0.5, density=1., kinematic_viscosity=1.,
    ext_force_density=[0.02, 0.04, 0.06], single_precision=False)
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=1.)
p = system.part.add(pos=3 * [-1e-30], v=[-1, 0, 0],
                    propagation=espressomd.propagation.Propagation.TRANS_LB_TRACER)
system.integrator.run(1)
print(p.v)

Output:

RuntimeError: Cannot apply force to LB

Expected output:

[0.0003 0.0006 0.0009]

The virtual_sites_tracers.py test fails when agrid=0.5 in virtual_sites_tracers_common.py.

This bug most likely affects all waLBerla versions of ESPResSo. Here is the same bug in the original commit that introduced waLBerla on the python branch:

https://github.com/espressomd/espresso/blob/3fd170980bed30c430a9b0264e9504632b4b7326/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp#L74-L76