Implement new LB boundary infrastructure

[jngrad]

Old infrastructure

The LBBoundaries framework offers a LBBoundary wrapper class for Shape objects and a container to store and serialize LBBoundary objects.

Advantages:

• removal of individual boundaries in an arbitrary order (the geometry is recalculated after each addition or removal)
• the sequence of Shapes that were used to build that geometry can be recovered, including their slip velocities
• automated checkpointing

Disadvantages:

• adding/removing a boundary causes the entire boundary geometry to be reconstructed, which silently erases per-node slip velocity information, see https://github.com/espressomd/espresso/issues/4252#issuecomment-948893947
• the framework depends on the ObjectList manager, which has a high maintenance cost

New infrastructure

It is now possible to pass Shape objects directly to the LBFluidWalberla instance to build the desired geometry without the above-mentioned disadvantages. It's also possible to pass a list of node indices to the LBFluidWalberla instance, or to update nodes individually via VelocityBounceBack.

It is however not possible to recover the sequence of shapes that were used to build that geometry, so if this information is important, one has to keep the Shape objects in a Python list and manually register this list in the checkpoint file. Likewise, the visualizer cannot know the sequence of shapes that was used to build the boundary geometry, so it needs to be rasterized with LB_draw_node_boundaries.

Outline

The new framework could be completed with the following steps:

1. convert all uses of system.lbboundaries.add(LBBoundary(shape=my_shape), velocity=my_vel) to lb_fluid.add_boundary_from_shape(shape=my_shape, velocity=my_vel) in the testsuite, samples, tutorials and user guide using 61b2e0549d083ef9 as a template
2. remove the LB_draw_boundaries feature from the OpenGL visualizer: the boundary geometry is now rasterized with the LB_draw_node_boundaries feature
3. move the VelocityBounceBack class from lbboundaries.py to lb.pyx
4. remove the LBBoundary and LBBoundaries classes in the core, script interface and python interface
5. remove the LB_BOUNDARIES macro from the core, testsuite, samples, docs and maintainer/configs/*.hpp
6. decide where to move the edge_detection() and calc_cylinder_tangential_vectors() utility functions from lbboundaries.py; these are used to set up a cylindrical slip velocity in the testsuite and samples

[jngrad]

2. remove the LB_draw_boundaries feature from the OpenGL visualizer: the boundary geometry is now rasterized with the LB_draw_node_boundaries feature

Actually, the LB_draw_node_boundaries is quite inefficient and slows down the visualizer considerably. Maybe keep LB_draw_boundaries and make it instead draw black dots at the center of every boundary nodes, i.e. offset = 0.5. The code is already available in the shape constraint visualization function for non-trivial shapes (e.g. the Union shape).

[RudolfWeeber]

lb_fluid.add_boundary_from_shape(shape=my_shape, velocity=my_vel): Some thoughtss on names: add_boudnary could suggest that the boundary is a cohesive object that one can talk to.

The most clear syntax would probably be

sphere = Sphere(...)
boundary_nodes = lb_fluid.nodes_within_shape(sphere)
n.boundary = VelocityBounceBack(v=...) for n in boudary_nodes

If we want a convenience funciton

• mark_as_boundary is more clear than add_boudnary, I guess
• the boundary type should probably be made explicit (VelocityBounceBack or VelocityBC)

[jngrad]

lb_fluid.add_boundary_from_shape(shape=my_shape, velocity=my_vel): Some thoughtss on names: add_boudnary could suggest that the boundary is a cohesive object that one can talk to.

I'm open to name changes.

The most clear syntax would probably be

sphere = Sphere(...)
boundary_nodes = lb_fluid.nodes_within_shape(sphere)
n.boundary = VelocityBounceBack(v=...) for n in boudary_nodes

I'm afraid people won't do this for shape-based constraints. Setting the slip velocity of the 4-roller mill with this syntax on a 128x128x4 LB grid took me 30 seconds on a Release build, and probably several minutes on a Debug build. Before the convenience functions were introduced, the boundary velocity assignment was wrapped in a tqdm progress bar to show users the script was still alive. For comparison, the LB fluid reaches a steady state in that geometry in less than 10 seconds. This is because, every time we update the velocity of one node, the entire velocity field needs to be recalculated from scratch.

If we want a convenience funciton

* mark_as_boundary is more clear than add_boudnary, I guess

* the boundary type should probably be made explicit (VelocityBounceBack or VelocityBC)

Yes, that's a good idea.

What we would also need to figure out, is where to store the boundary field. Currently this is contained in the LB class, but the EK code also needs it, and EK doesn't necessarily have a blockforest. We could duplicate the lb.mark_as_boundary method as ek.mark_as_boundary, but the method would need to raise a runtime error if no blockforest exists, or if it has been invalidated (e.g. when the lb actor is removed from the list of active actors).

[RiccardoFrenner]

Replacing the calls to system.lbboundaries.add by lbf.add_boundary_from_shape in testsuite/python/lb_shear.py causes the test to fail. More specifically, one element of the pressure tensor does not match the expected value by an absolute difference of 0.0296 and relative difference of 4.29e-5. The failing element correlates with shear_plane_normal, e.g. with shear_plane_normal = [1,0,0] node_pressure_tensor[1][2] does not match, and with shear_plane_normal = [0,1,0] it's node_pressure_tensor[0][2].

Maybe you @jngrad or @RudolfWeeber have a clue what the problem is. Otherwise my next steps would be to debug it and see what exactly changes the values of the pressure tensor. I already looked at the code of both the old and the new method quite deeply but could not see what the issue could be.

Here is the diff causing the issue:

@@ -99,6 +99,7 @@ class LBShearCommon:

         self.lbf = self.lb_class(**LB_PARAMS)
         self.system.actors.add(self.lbf)
+        self.lbf.clear_boundaries()

         wall_shape1 = espressomd.shapes.Wall(
             normal=shear_plane_normal, dist=AGRID)
@@ -109,8 +110,10 @@ class LBShearCommon:
         wall2 = espressomd.lbboundaries.LBBoundary(
             shape=wall_shape2, velocity=.5 * SHEAR_VELOCITY * shear_direction)

-        self.system.lbboundaries.add(wall1)
-        self.system.lbboundaries.add(wall2)
+        self.lbf.add_boundary_from_shape(
+            wall_shape1, velocity=-.5 * SHEAR_VELOCITY * shear_direction)
+        self.lbf.add_boundary_from_shape(
+            wall_shape2, velocity=.5 * SHEAR_VELOCITY * shear_direction)

         t0 = self.system.time
         sample_points = int(H / AGRID - 1)

[jngrad]

I had a quick look into it, and I can't tell what the issue is. The lbf[:,:,:].velocity and lbf[:,:,:].is_boundary match perfectly for both ways of setting boundaries. The old LBBoundaries framework has additional steps because it reconstructs the entire sequence of shapes every time a new shape is added, but other than that the slip velocities are correct. This also can't be an issue with leaking boundaries because that was fixed on October 6.

[jngrad]

I printed out the content of the core hash map and can now report that the old LBBoundaries implementation was not correctly ported to walberla: it rasterizes the shape on all grid points including the ghost layer grid points without wrapping around the periodic boundary. Therefore wall1 will populate the slip velocities of slices with at x=-1 (ghost layer that copies the real layer at x=N) and x=0 (first real layer) when in reality it should populate x=0 (first real layer) and x=N+1 (ghost layer that copies the real layer at x=0).

I and @reinaual still don't understand why this is an issue, because fluid cells in contact with the boundaries physically cannot be in contact with the ghost layer boundaries (in this particular setup), therefore it shouldn't matter that the slip velocities are wrong in the old LBBoundaries implementation. Yet it does:

    def check_profile(self, shear_plane_normal, shear_direction):
        """
        Integrate the LB fluid and regularly compare with
        the exact solution.

        """
        self.tearDown()
        self.setUp()
        self.system.box_l = np.max(
            ((W, W, W), shear_plane_normal * (H + 2 * AGRID)), 0)
        self.system.actors.add(self.lbf)

        wall_shape1 = espressomd.shapes.Wall(
            normal=shear_plane_normal, dist=AGRID)
        wall_shape2 = espressomd.shapes.Wall(
            normal=-1.0 * shear_plane_normal, dist=-(H + AGRID))
        wall1 = espressomd.lbboundaries.LBBoundary(
            shape=wall_shape1, velocity=-.5 * SHEAR_VELOCITY * shear_direction)
        wall2 = espressomd.lbboundaries.LBBoundary(
            shape=wall_shape2, velocity=.5 * SHEAR_VELOCITY * shear_direction)
        wall_shape3 = espressomd.shapes.Wall(
            normal=shear_plane_normal, dist=0)
        wall_shape4 = espressomd.shapes.Wall(
            normal=-1.0 * shear_plane_normal, dist=-(H + 2*AGRID))
        wall3 = espressomd.lbboundaries.LBBoundary(
            shape=wall_shape3, velocity=-20. * SHEAR_VELOCITY * shear_direction)
        wall4 = espressomd.lbboundaries.LBBoundary(
            shape=wall_shape4, velocity=20. * SHEAR_VELOCITY * shear_direction)

        implementation = 'new'
        if implementation == 'old':
            self.system.lbboundaries.add(wall1)
            self.system.lbboundaries.add(wall2)
            self.system.lbboundaries.add(wall3)
            self.system.lbboundaries.add(wall4)
        else:
            self.lbf.add_boundary_from_shape(
                wall_shape1, velocity=-.5 * SHEAR_VELOCITY * shear_direction)
            self.lbf.add_boundary_from_shape(
                wall_shape2, velocity=.5 * SHEAR_VELOCITY * shear_direction)
            self.lbf.add_boundary_from_shape(
                wall_shape3, velocity=-20. * SHEAR_VELOCITY * shear_direction)
            self.lbf.add_boundary_from_shape(
                wall_shape4, velocity=20. * SHEAR_VELOCITY * shear_direction)

With implementation = 'new' the test will fail but this is most likely just a numerical precision issue (the two extra walls don't change any of the tensor values, because their slip velocities are not stored in the hash map). With implementation = 'old' the test will completely break, even though the extra 2 walls (with extreme slip velocities) only exist in the inaccessible ghost layers.

@reinaual There must be something wrong in the streaming step.

[jngrad]

@RiccardoFrenner for the time being, you can comment out the failing test in testsuite/python/CMakeLists.txt and move forward with the changes.

[RiccardoFrenner]

Okay, sure.

[RiccardoFrenner]

@jngrad How should the LBBoundary::get_force() method be replaced? I tried with something like np.sum([n.boundary_force for n in lbf.get_nodes_in_shape(shape)], axis=0) but didn't get the same results, as can be seen in e.g. testsuite/python/lb_boundary_volume_force.py. The reason for that is, that LBBoundary::get_force() also sums up boundary forces in the ghost layers.

For my understanding: If each direction is periodic, as in lb_boundary_volume_force.py, boundary nodes at the edge of the simulation box do not contain the complete boundary force, since the other part is stored in the ghost layer of the force field?

[jngrad]

You can comment out the lines in python tests that query boundary forces. These forces are not reliable at the moment, because the dynamic UBB doesn't natively support a macroscopic value getter to sum up the forces. That might change in the near future.

[RiccardoFrenner]

The last thing I need to do is to move the edge_detection and calc_cylinder_tangential_vectors functions from lbboundaries.py. Since they are both LB specific and lb.pyx is the only python file with LB code, I think they should be put in there. Otherwise a new file e.g. lb_utils.py would also be an option.

[jngrad]

You can move it to lb.pyx for now. If we later need it in EK, @reinaual can decide for a better place.

[RiccardoFrenner]

Where/how should I do my pull request?

[jngrad]

Simply post your branch name here, I'll review and merge it from the command line.

[RiccardoFrenner]

https://github.com/RiccardoFrenner/espresso/tree/fix-4381-rebase

[jngrad]

Thanks! The merge was successful.