active matter

[alibh1996]

Hi. I am currently exploring the fascinating world of active particles using the lattice Boltzmann method. For my research, I have installed both Espresso 4.2.2 and Espresso Python, and have utilized Walberla to simulate Poiseuille flow.
In my pursuit of applying torque from Poiseuille flow to raspberry particles, In order to apply the torque from Poiseuille flow to the raspberry particles, I have used virtual sites according to instructions in "Rigid arrangements of particles". But the problem is that Walbrella works only in Espresso Python and does not work in Espresso 4.2.2 and virtual sites work in Espresso 4.2.2 but does not work in Espresso Python. Consequently, the torque does not seem to be applied to the active particles as intended (my code is available in the attachment).
I would be deeply grateful if you could assist me with the following inquiries:

1. Could you kindly elucidate why torque is not being applied to my particles?
2. Have I conducted the simulation of active particles in a Poiseuille flow correctly? If there are any discrepancies, may I request your guidance on where I might have erred?
3. Is it possible to have walberla in Espresso 4.2.2?
4. Is the use of Walberla imperative for creating a Poiseuille flow, or might there be a simpler alternative method available?

Your expertise and assistance in this matter would be immensely appreciated.

my code:

import espressomd.lb
import espressomd.math
# import matplotlib.pyplot as plt
# plt.rcParams.update({'font.size': 18})
import espressomd.swimmer_helpers
import tqdm
import numpy as np
# import espressomd.visualization
import espressomd.shapes
import espressomd
import espressomd.observables
import espressomd.accumulators
# required_features = ["VIRTUAL_SITES_RELATIVE", "TRANS_LB_TRACER"]

# # import espressomd.virtual_sites
# espressomd.assert_features(required_features)


HYDRO_PARAMS = {'box_l': [50,50,200],
                'time_step': 0.01,
                'skin': 1,
                'agrid': 1,
                'dens': 0.1,
                'visc': 5,
                'kT': 1,
                'gamma_rotation': 1,
                'gamma': 1,
                'mass': 1,
                'dipole_length': 4,
                'f_swim': 5,
                'dipole_particle_type': 1,
                'corr_tmax': 100
                }

HYDRO_N_STEPS = 10000

system = espressomd.System(box_l=HYDRO_PARAMS['box_l'])
system.cell_system.skin = HYDRO_PARAMS['skin']
system.time_step = HYDRO_PARAMS['time_step']
# system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()

system.periodicity = [True, True, True]

RECT_PARAMS = {'wca_sigma': 0.5,
               'wca_epsilon': 0.1,}

RECT_STEPS_PER_SAMPLE = 100
RECT_N_SAMPLES = 500

TYPES = {'particles': 0,
         'boundaries': 1}
BOX_L = HYDRO_PARAMS['box_l'][0]

wall_right= espressomd.shapes.Wall(dist=0., normal=[0, 1, 0])
system.constraints.add(shape=wall_right , particle_type=TYPES['boundaries'])
wall_left= espressomd.shapes.Wall(dist=-50., normal=[0, -1, 0])
system.constraints.add(shape=wall_left , particle_type=TYPES['boundaries'])
wall_up= espressomd.shapes.Wall(dist=-50., normal=[-1, 0, 0 ])
system.constraints.add(shape=wall_up , particle_type=TYPES['boundaries'])
wall_down= espressomd.shapes.Wall(dist=0., normal=[1, 0, 0])
system.constraints.add(shape=wall_down , particle_type=TYPES['boundaries'])

# LB parameters
AGRID = 0.5
VISCOSITY = 2.0
FORCE_DENSITY = [0, 0.0, 0.1]
DENSITY = 1.5

# LB boundary parameters
WALL_OFFSET = AGRID


system.min_global_cut = HYDRO_PARAMS['dipole_length']


lbf = espressomd.lb.LBFluidWalberla(agrid=HYDRO_PARAMS['agrid'],
                                    density=HYDRO_PARAMS['dens'],
                                    kinematic_viscosity=HYDRO_PARAMS['visc'],
                                    tau=HYDRO_PARAMS['time_step'],
                                    ext_force_density=FORCE_DENSITY,
                                    kT=HYDRO_PARAMS['kT']
                                    )
system.lb = lbf

top_wall = espressomd.shapes.Wall(normal=[1, 0, 0], dist=WALL_OFFSET)
bottom_wall = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET))
left_wall = espressomd.shapes.Wall(normal=[0, 1, 0], dist=WALL_OFFSET)
right_wall = espressomd.shapes.Wall(normal=[0, -1, 0], dist=-(BOX_L - WALL_OFFSET))

lbf.add_boundary_from_shape(top_wall)
lbf.add_boundary_from_shape(bottom_wall)
lbf.add_boundary_from_shape(left_wall)
lbf.add_boundary_from_shape(right_wall)

system.thermostat.set_lb(LB_fluid=lbf, gamma=HYDRO_PARAMS['gamma'], seed=42)


# add particles

pos=[20,25,50] 
for i in range(1):     
        particle = system.part.add(pos=pos, 
                                #    director=director,
                                   mass=HYDRO_PARAMS['mass'],
                                   rotation=3*[True],
                                #    rinertia=3*[1], 
                                   swimming={'f_swim': HYDRO_PARAMS['f_swim']})
        rel_offset1= [0, 2, 0.]
        p1 = system.part.add(pos=particle.pos + rel_offset1)
        p1.vs_auto_relate_to(particle)
        rel_offset2 = [0, -2, 0.]
        p2 = system.part.add(pos=particle.pos + rel_offset2)
        p2.vs_auto_relate_to(particle)
        rel_offset3 = [0, 4, 0]
        p3 = system.part.add(pos=particle.pos + rel_offset3)
        p3.vs_auto_relate_to(particle)
        rel_offset4 = [0, -4, 0]
        p4 = system.part.add(pos=particle.pos + rel_offset4)
        p4.vs_auto_relate_to(particle)
system.integrator.run(0, recalc_forces=True)


system.non_bonded_inter[TYPES['particles'], TYPES['boundaries']].wca.set_params(
    epsilon=RECT_PARAMS['wca_epsilon'], sigma=RECT_PARAMS['wca_sigma'])


import pathlib
vtk_base_dir = pathlib.Path('./vtk_out').resolve()
vtk_identifier = 'RESULTS_FLOW_FIELD'
vtk_outdir = vtk_base_dir / vtk_identifier
vtk_outdir.mkdir(parents=True, exist_ok=True)
lb_vtk = espressomd.lb.VTKOutput(lb_fluid=lbf, 
                                 identifier=vtk_identifier,
                                 observables=["velocity_vector"],
                                 base_folder=vtk_base_dir.as_posix(),
                                 prefix='lb_velocity')
lbf.add_vtk_writer(vtk=lb_vtk)

for i in tqdm.tqdm(range(1000)):
    system.integrator.run(10)
    system.part.writevtk(vtk_outdir / f'position_{i}.vtk', types=[0])
    lb_vtk.write()

[jngrad]

1. The 4.2.2 release and the python branch are quite different. I would recommend you pick one, because writing a script that would work for both versions is hard, even for seasoned ESPResSo developers. In addition, the way of setting the thermostat for virtual particles is completely different; see LB/Virtual sites: incomplete API change #4935 for a related issue.
2. The way the walls are currently set up, particles will bounce off, but the LB fluid won't bounce back. For that to happen, the walls should be offset by about one agrid away from the box side, inward. See 13.9.2. Shape-based boundary conditions for more details, and samples/visualization_poiseuille.py for a minimal simulation example to set up a Poiseuille flow. The sample exists in both the 4.2.2 and python versions of ESPResSo, although the 3D visualizer it uses is known to be glitchy.
3. waLBerla cannot be used in 4.2.2.
4. The 4.2.2 release has a different LB engine, but is just as capable of generating a Poiseuille flow.

[alibh1996]

Hello. Thank you very much for your tips. I made the necessary changes according to your instructions and got a good result. As I said before, I am looking to simulate active particles in a channel containing Poiseuille flow with the LB method. But the problem that exists now is that the velocity of the fluid increases with each step when it moves forward, and after a while the velocity of the fluid becomes very high. I want to have a fluid with constant velocity . Am I missing something in my code?
thanks for your help

my code:

import espressomd.lb
import espressomd.math
import tqdm
import numpy as np
from espressomd.shapes import Wall
import espressomd.shapes
import espressomd
import espressomd.observables
import espressomd.accumulators
required_features = ["VIRTUAL_SITES_RELATIVE", "TRANS_LB_TRACER"]
from espressomd.virtual_sites import VirtualSitesRelative
import espressomd.virtual_sites


HYDRO_PARAMS = {'box_l': [50,50,200],
                'time_step': 0.01,
                'tau': 0.01,
                'skin': 1,
                'agrid': 2,
                'dens': 1,
                'visc': 1,
                'kT': 0,
                'gamma_rotation': 1,
                'gamma': 1,
                'mass': 1,
                'dipole_length': 4,
                'f_swim': 0,
                'force_density':[0.0, 0.0, 1],
                'dipole_particle_type': 1,
                'corr_tmax': 100
                }


HYDRO_N_STEPS = 10000

system = espressomd.System(box_l=HYDRO_PARAMS['box_l'])
system.cell_system.skin = HYDRO_PARAMS['skin']
system.time_step = HYDRO_PARAMS['time_step']

system.periodicity = [True, True, True]

RECT_PARAMS = {'wca_sigma': 0.5,
               'wca_epsilon': 0.1,}

RECT_STEPS_PER_SAMPLE = 100
RECT_N_SAMPLES = 500

TYPES = {'particles': 0,
         'boundaries': 1}
BOX_L = HYDRO_PARAMS['box_l'][0]


# LB boundary parameters
WALL_OFFSET = HYDRO_PARAMS['agrid']
system.min_global_cut = HYDRO_PARAMS['dipole_length']


n=10
Z = np.linspace(5, 45, n)
Y = np.linspace(5, 45, n)
a = []
for z in Z:
    for y in Y:
        director = np.array([1,0,0])
        pos = [y, z, 5]
        particle= system.part.add(
            pos=pos,
            director = np.array([0,0,1]),
            swimming={'f_swim': HYDRO_PARAMS ['f_swim']},
            mass= HYDRO_PARAMS ['mass'],
            rotation=3*[True],
            rinertia=3*[1]
            )


lbf = espressomd.lb.LBFluid(agrid=HYDRO_PARAMS['agrid'],
                                    dens=HYDRO_PARAMS['dens'],
                                    visc=HYDRO_PARAMS['visc'],
                                    tau=HYDRO_PARAMS['tau'],
                                    ext_force_density=HYDRO_PARAMS['force_density'],
                                    # ext_force_density = [0, 0, fluid_force],
                                    kT=HYDRO_PARAMS['kT'],
                                    seed=123
                                    )
system.actors.add(lbf)
system.thermostat.set_lb(LB_fluid=lbf, gamma=HYDRO_PARAMS['gamma'], seed=123 , act_on_virtual=False)


bound_walls = [None] * 4
bound_walls[0] = Wall(normal=[1, 0, 0], dist=WALL_OFFSET)
bound_walls[1] = Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET))
bound_walls[2] = Wall(normal=[0, 1, 0], dist=WALL_OFFSET)
bound_walls[3] = Wall(normal=[0, -1, 0], dist=-(BOX_L - WALL_OFFSET))
for wall in bound_walls:
    system.lbboundaries.add(espressomd.lbboundaries.LBBoundary(shape=wall))
    system.constraints.add(shape=wall, particle_type=TYPES['boundaries'], penetrable= True)


system.non_bonded_inter[TYPES['particles'], TYPES['boundaries']].wca.set_params(
    epsilon=RECT_PARAMS['wca_epsilon'], sigma=RECT_PARAMS['wca_sigma'])


import pathlib
vtk_base_dir = pathlib.Path('./vtk_out').resolve()
vtk_identifier = 'RESULTS_FLOW_FIELD'
vtk_outdir = vtk_base_dir / vtk_identifier
vtk_outdir.mkdir(parents=True, exist_ok=True)

for i in tqdm.tqdm(range(1000)):
    system.integrator.run(10)
    system.part.writevtk(vtk_outdir / f'position_{i}.vtk', types=[0])
    lbf.write_vtk_velocity(f'{vtk_outdir}/velocity_{i}.vtk')

[jngrad]

But the problem that exists now is that the velocity of the fluid increases with each step when it moves forward, and after a while the velocity of the fluid becomes very high. I want to have a fluid with constant velocity .

You are currently simulating in the supersonic regime (Mach 2), which is not something ESPResSo was designed for. The results of your simulation will not match with experimental results. There is no guarantee the flow will ever reach terminal velocity. Consider decreasing the external force by two orders of magnitude to stay below the Mach 0.35 limit. You can plot the Mach number as a function of the position along the x-axis with this snippet:

import matplotlib.pyplot as plt
xdata = (np.arange(lbf.shape[0]) + 0.5)[1:-1] * HYDRO_PARAMS['agrid'] - 0.5 * system.box_l[0]
H = system.box_l[0] - 2. * HYDRO_PARAMS["agrid"]
DYN_VISC = HYDRO_PARAMS['dens'] * HYDRO_PARAMS['visc']
ydata = HYDRO_PARAMS['force_density'][2] * 1. / (2. * DYN_VISC) * (H**2 / 4. - xdata**2)
lb_speed_of_sound = (HYDRO_PARAMS['agrid'] / HYDRO_PARAMS['tau']) / np.sqrt(3.)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xdata, ydata / lb_speed_of_sound)
ax.legend()
ax.set_xlabel("x")
ax.set_ylabel("Mach number")
plt.show()

Once the external force is within a valid range for ESPResSo, you will also need to make the simulation run for a longer time. It is currently not long enough to see the fluid velocity converge to terminal velocity. If you replace the last for loop in your script by the following loop, you will see the fluid velocity reach a plateau:

xdata = []
lb_node = lbf[lbf.shape[0]//2, lbf.shape[1]//2, lbf.shape[2]//2]
for i in tqdm.trange(20):
    system.integrator.run(steps=200)
    xdata.append(lb_node.velocity)
xdata = np.array(xdata)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xdata[:, 0], label=r"$v_x$")
ax.plot(xdata[:, 1], label=r"$v_y$")
ax.plot(xdata[:, 2], label=r"$v_z$")
ax.legend()
ax.set_xlabel("simulation time")
ax.set_ylabel("velocity")
plt.show()