Immersed boundary method: external forces and constraints may be missed

[RudolfWeeber]

@sgekle, @fweik

The CPU implementation likely misses part of the forces contributed by constraints and the ext_force property of particls close to node boundaries, both the actual particles and its ghosts are coupled to the fluid, evenn on a single MPI rank. IBM applies the content of p.f.f to the fluid. On ghosts, this, however, does not contain single particle forces such as those from constraints and the particle's ext_force property. This has to be so, otherwise, those forces would be applied several times in the integrator after summing up the forces contributed by the ghosts.

With the Walberla branch, it is possible to test the individual coupling schemes because the last applied forces are accessible in the tests. Above issue was discovered when writing down the test for inertialess tracers.

@sgekle, I'll sort out the Walberla branch. If you'd like to fix the issue in 4.1, you would need to manually add external force and constraint contributions before coupling ghosts to the fluid.

The GPU version is probably not affected, since only non-ghost particles are considred when collecting particles for transfer to the gpu.

[RudolfWeeber]

Particles within one agrid of a node- or domain boundary are affected. Between 1/2 and 7/8th of the forces are neglected, depending on whether the particle is next to a wall/edge/corner.

[sgekle]

Dear Rudolf,

could you point me to the exact position in the source code where you think that the problem is?

Thanks,

  Stephan

Am 29.05.20 um 20:29 schrieb RudolfWeeber:

Particles within one agrid of a node- or domain boundary are affected. Between 1/2 and 7/8th of the forces are neglected, depending on whether the particle is next to a wall/edge/corner.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/espressomd/espresso/issues/3744#issuecomment-636121618, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABKTVKBI43EE7I3NHPIVF4DRT75INANCNFSM4NOG6EIQ.

[RudolfWeeber]

Stephan, I don't know, as I don't fully understand the code in 4.1.2. We'll try to sort it out in the Walberla branch first, because there, the code is simpler and the individual coupling step can be tested directly. https://github.com/RudolfWeeber/espresso/blob/walberla/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp Once we have a olution there, we can look into backporting it.

[RudolfWeeber]

@ffweik, I copied the ghost force code from vs relative as suggested, but it didn't do the trick.

The contents of p.p.ext_force still is not present as part as p.f.f on ghost particles, so if a particle close to a boudnary is coupled, all lb cells within range from the ghosts don't receive that part of the external force.

It was also my understanding that the ghost force reduction only works one way, i.e., it collects forces from all ghosts and stores the sum on the physical particles, but doesn't store the sum on the ghosts. Is that incorrect?

[fweik]

Upon further inspection it might be the case that the code in the virtual sites is in the wrong direction, sorry for that. Just to clarify: You want that the forces on the ghost particles are the same as on the corresponding real particles? If this is the case then you can set the forces on the ghosts to zero, and then do a regular ghost update with the force as data part. (You have to set them to zero first because the forces are added, not assigned).

[RudolfWeeber]

@fweik, using ghost update with forces as data parts worked. Thanks.

Offline discussion yielded that the VS relative version is likely correct, because there, forces are transferred from physical non-ghost particles to virtual ghost particles. In that case, the standard ghost force reduction is the correct thing to do.

Now looking into 4.1.3rc.

[RudolfWeeber]

Stephan (@sgekle), could you please review the following code and let me know if that's what the method is supposed to do: https://github.com/RudolfWeeber/espresso/blob/walberla/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp and the single step test https://github.com/RudolfWeeber/espresso/blob/walberla/testsuite/python/virtual_sites_tracers_common.py#L73

Should the update of the particle velocity happen before or after the the point force is applied to the fluid? It matters, because of the f/2 correction.

If you like, we can organize a video call to discuss further.

@jngrad, I don't think, this will make it into 4.1.3. Unfortunately, the test cannot be backported to older versions of Espresso, because the existing LB does not consider point forces in lb velocity interpolation and has no Python access to the full force double-buffering which only happens in the inertialess tracers method, in any case.

[sgekle]

Dear Rudolf,

sorry for the late reply. I tried to figure out the exact order that we are currently using. It appears to be the following: 1) force calculation 2) forces into fluid 3) LB streaming & collision 4) velocity interpolation & position update

I looked at VirtualSitesInertialessTracers.cpp https://github.com/RudolfWeeber/espresso/blob/walberla/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp where indeed the order seems to be different. To fix this, I would suggest to move the velocity interpolation lb_lbinterpolation_get_interpolated_velocity() to the beginning of after_lb_propagation(). This should fix the order and in addition it ensures that IBM particles are advected with exactly the local fluid velocity because interpolation & advection directly follow each other and no intermediate modifications to the particle velocity can occur. Furthermore, when calling the add_md_force there is a minus in front of the paticle force. Why is this? In IBM the positive force should be added, i.e. if the particle experiences a force to the right, the fluid should also feel a force to the right.

Finally, two small remarks:

• is it still possible to change the name? The term "tracers" suggests that such particles only follow the velocity field passively but do not modify the velocity field. Therefore, it might be clearer to call them "inertialessParticles" instead.
• I never fully understood what the "lattice_switch" variable does, but it does not appear to be connected to invidual particles. Therefore, I wonder whether the comparison "if (lattice_switch == ActiveLB::NONE)" really needs to be inside the particle loop? Wouldn't it be more logical (and faster) to check this condition once at the beginning of the function and not at every iteration?

Best regards,

  Stephan

Am 23.06.20 um 09:35 schrieb RudolfWeeber:

Stephan (@sgekle https://github.com/sgekle), could you please review the following code and let me know if that's what the method is supposed to do: https://github.com/RudolfWeeber/espresso/blob/walberla/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp and the single step test https://github.com/RudolfWeeber/espresso/blob/walberla/testsuite/python/virtual_sites_tracers_common.py#L73

Should the update of the particle velocity happen before or after the the point force is applied to the fluid? It matters, because of the f/2 correction.

If you like, we can organize a video call to discuss further.

@jngrad https://github.com/jngrad, I don't think, this will make it into 4.1.3. Unfortunately, the test cannot be backported to older versions of Espresso, because the existing LB does not consider point forces in lb velocity interpolation and has no Python access to the full force double-buffering which only happens in the inertialess tracers method, in any case.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/espressomd/espresso/issues/3744#issuecomment-647966248, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABKTVKAZGTPTDAJMUXQZAELRYBLMRANCNFSM4NOG6EIQ.

[RudolfWeeber]

Stephan, I moved the particle velocity update as suggested. I also adapted the test. Can you please take an other look whether this is what you meant?

lattice_switch contains the info which LB implementaiton (if any) is active. It probably should be renamed, but that's a different matter.

There is a per particle check, because switching the virtual sites to inertialess tracers + no LB + virtual particle is illegal.

The minus in add_md_force() is confusing but correct as far as I can see. There is an other minus in there, somewhere. This needs to be changed, still.

The name is, unfortunately, fixed for now, since it is already there in 4.0 and 4.1. We maintain interface compatibility, where possible.

More fundamentally, being a tracer is not a property of the particle in the first place. Actually, a different integration/propagation scheme is used for those particles. Once/when this is represented in the core, the interface will change once more.