Perfect restart for granular pair style

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

I have questions concerning the restart for granular pair style. It is written in the manual that granular pair styles cannot restart exactly because they uses half step velocities for force computation. Or Lammps do not store forces at the end of the run but recomputed them at the beginning of the run, of course with wrong full step velocities.

The problem is that for my simulations the newly recomputed force has a relative error (compared to the old one) of about 1e-3, which can seem small but in my case is sufficient to screw up my previously barostated state. Therefore, I want to change the code to obtain a perfect restart for granular pair styles.

The simplest solution I see is just to add a fix (that will be directly called from the pair gran style like the "fix shear history") which will store the forces at the end of the run (at the end_of_step() stage of the last timestep of the run) and restore them at the setup() stage substituting the forces computed in the initial phase.

I have two simple questions before starting the implementation:

-Did I correctly understand the problem and the solution to obtain exact restarts?
-Is there a more elegant way to do it than the one I proposed?

Thanks,
Anton Peshkov

2

To be honest, I doubt the procedure as is incurs any more error in the system trajectory than changing the communication procedure (box decomposition). If you’re force is changing significantly between half and full steps, that’s probably more an indication that something else is going awry (esp. for granular particles with relatively smooth potentials on the scale of the amount a particle should move in a timestep). One would probably need more info to truly understand your problem/difficulties and to comment if the work would be necessary(what barostat, etc) or beneficial.

Nonetheless, if its easy to code, than there’s certainly no harm in it. Your reasoning for the start-up seems consistent with velocity-verlet.

3

Hi Eric,

I study the rheology of a granular suspension with viscous force acting on granular particles under shear stress with Lees-Edwards boundary conditions and an imposed pressure using the Berendsen barostat (in 2d with shear in x direction). For some parameters the system volume and the stress in the xy direction takes a long time of several 1e8 or 1e9 steps to equilibrate (while the pressure equilibrates very quickly). After a restart, the volume of the system quickly increase by less than one percent but it is enough to sensibly change the stress in xy direction by about 5 percent (the pressure stay constant during this time). After that, it takes once again a long time for the volume and stress in xy direction to converge to a stationary value. I attached an example plot of the pressure, stress in xy direction and volume for two successive restarts (in red and green respectively).

Given that the pressure stays constant, I do not think that this is a barostat problem (I will however try Nose-Hoover). I could probably solve the problem by just reducing the timestep, but it is already very low (5e-5 for a k_n of 1e4). For evident time reasons I don’t want to put it even lower, and I don’t see any option in Lammps that would allow me to change the timestep during the run (to drastically reduce it just at the beginning and at the end of the run). Apart from that restart problem, the system behaves correctly.

As of my solution, I already found a problem in it. If between runs they will be some changes to the system that will change the forces (like new atoms added or some fixes that change the force), they will not be taken into account. Therefore, a better solution will be to:
-save the half step velocity at post_force() stage of the last timestep of the run
-substitute this velocity at setup_pre_exchange() stage
-put back the full step velocity at setup() stage

Anton

pressure.png

stress_xy.png

volume.png

4

Hi Eric,

I study the rheology of a granular suspension with viscous force acting on granular particles under shear stress with Lees-Edwards boundary conditions and an imposed pressure using the Berendsen barostat (in 2d with shear in x direction).

By Lees-Edwards do you mean fix deform or some other implementation? I assume by viscous forces you mean, the dashpot at particle contact. I hope you are not referring to ‘fix viscous’ which would seem wholly unphysical for a granular shearing problem. If you are interested in hydrodynamics, fast lubrication dynamics are available. However, under this implementation, strictly speaking, granular collisions are not allowed because it is Stokes flow. Consequently, rigid particles may never come into contact. There are reasons for this to physically break down, but the no overlap constraint is an assumption in Stokesian Dynamics(which fast lubrication dynamics is based on).

For some parameters the system volume and the stress in the xy direction takes a long time of several 1e8 or 1e9 steps to equilibrate (while the pressure equilibrates very quickly). After a restart, the volume of the system quickly increase by less than one percent but it is enough to sensibly change the stress in xy direction by about 5 percent (the pressure stay constant during this time). After that, it takes once again a long time for the volume and stress in xy direction to converge to a stationary value. I attached an example plot of the pressure, stress in xy direction and volume for two successive restarts (in red and green respectively).

These plots are troubling… I would expect for granular systems, if the system dilates(and positions are rescaled?) you would see a significant decrease in the magnitude of the shear stress. Of course if you’re using fix viscous this might happen due to increased velocities at the edge of the box in the shear direction. Anyways, to me, it doesn’t seem rational just yet to place blame on the barostat. Your pressure fluctuations seem to behave consistently with dilation. But I am curious why the fluctuations seem to build up, someone more familiar with these systems could probably comment

Given that the pressure stays constant, I do not think that this is a barostat problem (I will however try Nose-Hoover). I could probably solve the problem by just reducing the timestep, but it is already very low (5e-5 for a k_n of 1e4). For evident time reasons I don’t want to put it even lower, and I don’t see any option in Lammps that would allow me to change the timestep during the run (to drastically reduce it just at the beginning and at the end of the run). Apart from that restart problem, the system behaves correctly.

Why not just give separate small runs at the end and beginning? One thing you might do, if a barostat is the issue, would be to increase the pdamp at start up to prevent such huge dilations at start-up, and hopefully getting around any force at restart issues. Of course decreasing a time-step here would also seem to solve your problem, if force at restart is indeed the problem. These may of course just be skirting around any actual underlying problem.

As of my solution, I already found a problem in it. If between runs they will be some changes to the system that will change the forces (like new atoms added or some fixes that change the force), they will not be taken into account. Therefore, a better solution will be to:
-save the half step velocity at post_force() stage of the last timestep of the run
-substitute this velocity at setup_pre_exchange() stage
-put back the full step velocity at setup() stage

You might contact Dr. Sundaresan’s group at Princeton. I do know they’ve done some barostatted sheared granular systems using lammps, though for cohesive systems (circa 2008).

5

Anton,

Having thought on this a tick longer, I am somewhat convinced your problem lies not in your force calculation but in those giant pressure fluctuations. My two cents.

6

Hi Eric,

Yes, I use fix deform with velocity remapping for Lees-Edwards.
But, not I am not interestedin the dashpot model. I use a modified version of “fix viscous” (with the flow profile adapted to the shear rate). I study a modified version of a granular suspension studied in PRL 109, 118305 (arxiv:1204.2732). As you correctly pointed out in this case the fluid viscosity, increase the shear stress.

The increase in pressure fluctuations during the time is consistent with the decrease of the system volume. However, I am not sure what do you call giant pressure fluctuations? The relative fluctuations are less than one percent of the target pressure. Perhaps it is still not normal? I am new to granular systems.

I must admit that I did not read the manual carefully enough and missed the ‘pre’ and ‘post’ options of the run command, which can solve my problems without any need to code a perfect restart, by just making small runs with greatly reduced time step at the beginning and at the end.

Thanks for your suggestions,
Anton Peshkov

Hi Eric,

I study the rheology of a granular suspension with viscous force acting on granular particles under shear stress with Lees-Edwards boundary conditions and an imposed pressure using the Berendsen barostat (in 2d with shear in x direction).

By Lees-Edwards do you mean fix deform or some other implementation? I assume by viscous forces you mean, the dashpot at particle contact. I hope you are not referring to ‘fix viscous’ which would seem wholly unphysical for a granular shearing problem. If you are interested in hydrodynamics, fast lubrication dynamics are available. However, under this implementation, strictly speaking, granular collisions are not allowed because it is Stokes flow. Consequently, rigid particles may never come into contact. There are reasons for this to physically break down, but the no overlap constraint is an assumption in Stokesian Dynamics(which fast lubrication dynamics is based on).

For some parameters the system volume and the stress in the xy direction takes a long time of several 1e8 or 1e9 steps to equilibrate (while the pressure equilibrates very quickly). After a restart, the volume of the system quickly increase by less than one percent but it is enough to sensibly change the stress in xy direction by about 5 percent (the pressure stay constant during this time). After that, it takes once again a long time for the volume and stress in xy direction to converge to a stationary value. I attached an example plot of the pressure, stress in xy direction and volume for two successive restarts (in red and green respectively).

These plots are troubling… I would expect for granular systems, if the system dilates(and positions are rescaled?) you would see a significant decrease in the magnitude of the shear stress. Of course if you’re using fix viscous this might happen due to increased velocities at the edge of the box in the shear direction. Anyways, to me, it doesn’t seem rational just yet to place blame on the barostat. Your pressure fluctuations seem to behave consistently with dilation. But I am curious why the fluctuations seem to build up, someone more familiar with these systems could probably comment

Given that the pressure stays constant, I do not think that this is a barostat problem (I will however try Nose-Hoover). I could probably solve the problem by just reducing the timestep, but it is already very low (5e-5 for a k_n of 1e4). For evident time reasons I don’t want to put it even lower, and I don’t see any option in Lammps that would allow me to change the timestep during the run (to drastically reduce it just at the beginning and at the end of the run). Apart from that restart problem, the system behaves correctly.

Why not just give separate small runs at the end and beginning? One thing you might do, if a barostat is the issue, would be to increase the pdamp at start up to prevent such huge dilations at start-up, and hopefully getting around any force at restart issues. Of course decreasing a time-step here would also seem to solve your problem, if force at restart is indeed the problem. These may of course just be skirting around any actual underlying problem.

As of my solution, I already found a problem in it. If between runs they will be some changes to the system that will change the forces (like new atoms added or some fixes that change the force), they will not be taken into account. Therefore, a better solution will be to:
-save the half step velocity at post_force() stage of the last timestep of the run
-substitute this velocity at setup_pre_exchange() stage
-put back the full step velocity at setup() stage

You might contact Dr. Sundaresan’s group at Princeton. I do know they’ve done some barostatted sheared granular systems using lammps, though for cohesive systems (circa 2008).