Hi all.
I have written a pair potential which does local damping (Kelvin damping a la the fracture simulations by Holian et. al.). It is implemented in an analogous way to pair_dpd.cpp, but damps the relative velocity vector rather than just the radial component. It seems to work fine for the most part, but there is some bizarre behavior which I’m hoping someone has some thoughts on. It works well enough otherwise that I’ve been using it for a couple of months without noticing this bug.
In order to see the strange behavior, it seems like I need all of the following 3 things:
- a multi-processor simulation
- a volume/rescale fix on one direction and a zero-pressure nph on the other (simulations are 2D) (volume/rescale in both directions won’t produce the bug… maybe non-zero confining pressures will produce the bug too, but I haven’t tried)
- my home-brewed Pair class which implements the Kelvin dissipation
The behavior is as follows: eventually, one of the pairs of particles which lies near a CPU boundary gets locked into a very tight orbit (few tenths of a particle diameter) at large velocities. The interaction potential is lj_smooth, and max velocities are order 10 in LJ units… period is only a few timesteps, and timesteps are 0.0025 in LJ units. The displacements seem to stay along a single ray with some slight drift, and the velocity just keeps changing sign but not direction. The particles never really get that far apart, and one would only notice something pathological going on by looking at, e.g., kinetic energy, since I don’t add the viscous contribution in to the virial.
I call it a “bug” rather than a numerical instability, since I can only get it to happen when I do a multi-processor run.
Here’s a relevant snippet from my pair class which comes after computing the usual pairwise forces:
-------------------- pair_lj_smooth_kelvin.cpp -------------
… snip …
// Now put on the drag forces:
if(visc[itype][jtype]!=0){
double delvx = v[i][0] - v[j][0];
double delvy = v[i][1] - v[j][1];
double delvz = v[i][2] - v[j][2];
f[i][0] -= delvxvisc[itype][jtype];
f[i][1] -= delvyvisc[itype][jtype];
f[i][2] -= delvz*visc[itype][jtype];
… snip …