I am a new user of LAMMPS in the context of my PhD project.
I recently learnt how to use it, but now I need to extend the code to add a feature.
This feature (which would be nice to have built-in the code) is for granular particles.
The idea is to use the following structure: fix ID group reinsert {dim} {dir} {value} {region} vel {vx} {vy} {vz}.
Basically, when a particle crosses the surface defined by {dim} {value} (example z 30) from the right direction (> if z>30 is the condition, < if z<30), the particle is reinserted with a random position in {region} with a given velocity [{vx};{vy};{vz}]. → fix reins all reinsert z > 30 reinsert_reg vel 0 -0.5 1.5
We can almost think of it as a periodic boundary condition, but with more flexibility since we can change the direction and velocities of the particles, the pbc cannot work for my purpose.
I made it work with a new fix function. It is:
a simple test of atom->x versus the value for each atom (x[i][2]>30 for i in nlocal here )
a sampling of the position in the reinsertion region with reins_region->generate_random(pos,true)
and a replacement of the old atom->v by [{vx};{vy};{vz}]
However, being a beginner, there are some points I am not completely familiar with. This function works with one core, but as soon as I use multiple processors (which work without the fix), it crashes printing “buffer overflow”.
This function, I think, should have the mask FINAL_INTEGRATE since I need the motion to happen first to know the next positions and compare them to the plane z=30 (first, is it right?).
Then, my guess is that the problem comes from the fact that the calculation of the new positions should move the atoms from one processor to another, and I did not implement that. Looking at the documentation and the forum, I saw that the potential candidates for doing that are: grow_arrays, copy_arrays, pack_exchange, unpack_exchange, pack_comm, unpack_comm, pack_reverse_comm, unpack_reverse_comm. However, I am very confused as of which ones I have to use for my case, as a lot of posts relate to store data over the timesteps, which is not my goal, just communicate the new positions/velocities to the right processor in some way. Then, what are the values I should include in the necessary function?
my first thought would be that the behavior you describe only makes sense if there is a preferential direction in which particles move (e.g. due to fix gravity). otherwise, how would you prevent a frequent removal/reinsertion of the same particles near the boundary?
Is there a specific reason to enforce a constant number of particles? Otherwise, it should be easy to set up the desired behavior with fix pour and either themo_modify lost ignore or fix evaporate.
another option to randomize could be using periodic boundaries and just defining two regions where randomized forces are applied
what you describe is essentially a combination of fix pour and fix evaporate in a single fix, but with the additional constraint of inserting exactly as many particles as you remove. This will require identifying connected subdomains and explicit MPI communication.
in general, please note that when you insert particles in a region you also need something that moves them away, or else you run a high risk of placing new particles on top of existing particles which will either have to result in an insertion retry or a diverging (= crashing) simulation.
Hello and thank you for your answer and suggestions. I will address your comments to give a bit more details.
It is the case, but also, since the reinsertion region is far from the removal plane, there is not way the same particle crosses the removal plane frequently.
My research will use the output of this calculation for something very specific which requires a constant amount of particles in the geometry. It also requires the id of the particles to be constant throughout their lifetime, and insertion/removal adds a high level of complexity here, which if I solve this parallelization problem is not worth exploring.
Randomized forces would prevent me from controlling the reinsertion velocity with accuracy, therefore I cannot use periodic boundary conditions. On top of that, the inlet and outlet do not have necessarily the same radius, or even the same orientation (we could think for instance of a conic geometry in which the particles are removed at one circular plane, but reinserted from an horizontal pipe somewhere in the cone).
Yes, I am expecting to need to identify connected subdomains and explicit MPI communication, but it is how to do it with the functions I mentioned above that is challenging to understand for me.
I have other sides of this function that I am not discussing here, but basically, I am working on an overlap check too, based on the fix pour function, but I did not want to complicate the discussion further.
If this gave you other thoughts, please do not hesitate to share. I am also looking forward to hearing about these MPI communication functions.
Neither of the functions you mention are MPI functions or applicable to your situation.
This whole process is far more complex than you currently seem to anticipate.
Basically, the only way I see this can work is to first remove particles according to the conditions on individual MPI processes and collect their atom IDs. Then do MPI communication (e.g. MPI_Allgather to collect that information on all MPI ranks). Then you can re-insert those particles with the same atom IDs (which is going to be very complex to avoid double insertion or missing one, if the insertion region is spanning multiple sub-domains which you must consider as a possibility). This is further complicated by the fact that many granular pair styles require to store particle history information (i.e. the positions of all neighbors at the previous step). that information must be updated/invalidated in suitable way or else the next force evaluation is likely going to be bogus.
This all must happen at the PRE_EXCHANGE stage, because after your reinsert event you must enforce rebuilding the neighbor lists, which includes rebuilding the local per-atom arrays.
Thank you, indeed it seems way more complex than anticipated.
Something I am confused about is that, before I read your last message, I found the “fix move” function, which kind of does what I was saying, since it moves the particles throughout the geometry and does not seem to fail with parallelization. Is there something I am missing here?
Fix move does not move atoms over a long distance (but you do). It is dependent on the regular neighbor list updates and the included redistribution of atoms into subdomains. Please see: 4.4.3. Neighbor lists — LAMMPS documentation
Fix move will cause “lost atoms” errors if atoms move too fast.
This reminds me to remind you that you need to reset the image flags for reinserted particles.
Thank you for all these answers. I think that I will try to use an intermediary solution, based on your suggestion:
First remove particles crossing the surface and count them
Then insert new particles with the same style as “pour”, but with a number of particles which is variable and depends on this counted number of crossing particles.
For now I will do with new IDs instead of the keeping the same and post-process to simplify the development work.
I think I will encounter many challenges on the way, so I might come back here to ask questions.
Thank you very much!
PS: out of curiosity, I was thinking of the periodic boundary conditions, how does it work then? Because it literally changes the position to a position far from the original one and it works across processors?
It is not far from the original position. Atoms can only move in small steps due to time integration.
Periodic boundaries in LAMMPS are effectively an extension of the domain decomposition. In DD you have a communication cutoff and all atoms up to that cutoff are recreated as “ghost” atoms in the neighboring subdomain. When you have PBC atoms can actually move a little bit outside their domain (for as long as there is a sufficient number of ghosts to cover all atoms up to the cutoff it doesn’t matter and that is why LAMMPS uses the “neighbor list skin” distance as additional buffer). When an atom wraps around the periodic boundary it actually doesn’t move but then its coordinate outside the box is wrapped back and its image flag is either incremented or decremented. So if you consider those unwrapped coordinates, the atoms just keep moving.