Hi list,

I’ve been using the GRANULAR package in LAMMPS to simulate granular flow just fine (both Windows and Linux). Recently I try to simulate annular shear flow, whose boundary control involves converting from Cartesian to Polar coordinates and defining/calculating many intermediate variables.

Now unlike previous jobs, the issue with this job is that it finishes alright on my Windows machine, but on Linux it runs for some time before exhausting all allocated memory and generating error: “slurmstepd: Exceeded job memory limit at some point. srun: error: task 2: Killed”. This looks like a memory leak to me. (LAMMPS 15May2015 version)

Naturally, I then tried it on newer versions (11-Aug17). It seems the memory leak problem is fixed, but weirdly, it only works under 10 processes, and when using more than 10 cores it would get segmentation fault before even starting the execution.

The script included. It would need the file “circular_adapt.data” to run though, which I can attach along if needed. Any input/idea would be appreciated!

# Annular Shear

units si

dimension 2

atom_style sphere

boundary p f p

newton off

comm_modify vel yes

variable dl equal 0.003

variable ds equal 0.002

variable err equal 0.1

variable dllower equal {dl}*(1-{err})

variable dlupper equal {dl}*(1+{err})

variable dslower equal {ds}*(1-{err})

variable dsupper equal {ds}*(1+{err})

variable cl equal 0.5

variable cs equal 0.5

variable d equal {cl}*{dl}+{cs}*{ds} # 0.0025

variable width equal 300*$d

variable height equal 500*d
variable rhos equal 2450
variable rhol equal {rhos}*{ds}/{dl} # 1633

# Define origin point

variable x_orig equal {width}/2.0 # x coordinate of origin variable y_orig equal {x_orig} # y coordinate of origin

read_data circular_adapt.data

neighbor $d bin

neigh_modify delay 0

# Interaction model parameters

variable Pwall equal 100 # units N/m

variable Kappa equal 10000.0 # Keep this number above 10000

variable kn equal {Pwall}*{Kappa} # units N/m

variable kt equal {kn}*0.5
variable e equal 0.1 # Restitution coefficient
variable ms equal 3.14/6*({ds})^3*v_rhos # 1.0257*10^(-5)

variable ml equal 3.14/6*({dl})^3*v_rhol # 2.3079*10^(-5)
variable m equal {cl}*{ml}+{cs}*${ms} # Average grain mass: 1.6668*10^(-5)
variable gamman equal -2*ln($e)/(3.14^2+(ln($e))^2)^0.5*(2

*v_kn/$m)^0.5*

variable gammat equal 50(9.8/$d)^0.5*0

variable gammat equal 50

# Interaction model

pair_style hybrid/overlay gran/hooke/history {kn} {kt} {gamman} {gammat} 0.4 1 gran/hooke/history {kn} {kt} 0 0 0.4 1 gran/hooke/history 0 0 0 0 0 0

pair_coeff 1 1 gran/hooke/history 1

pair_coeff 1 2 gran/hooke/history 1

pair_coeff 2 2 gran/hooke/history 1

pair_coeff 1 3 gran/hooke/history 2

pair_coeff 2 3 gran/hooke/history 2

pair_coeff 1 4 gran/hooke/history 1

pair_coeff 2 4 gran/hooke/history 1

pair_coeff 3 3* gran/hooke/history 3

pair_coeff 4 4 gran/hooke/history 3

group inside type 1 2

group outbound type 3

group inbound type 4

fix 1 inside nve/sphere

compute stress all stress/atom NULL pair

fix 2 all enforce2d

# Convert to polar coordinates

variable r_coord atom sqrt((x-{x_orig})*(x-{x_orig})+(y-{y_orig})*(y-{y_orig}))

variable theta_coord atom atan2((y-{y_orig}),(x-{x_orig}))

variable f_r atom fx*cos(v_theta_coord)+fy*sin(v_theta_coord)

variable f_theta atom fy*cos(v_theta_coord)-fx*sin(v_theta_coord)

# Use velocity control method to apply pressure

compute 1 all erotate/sphere

compute R_outer outbound reduce ave v_r_coord # Radius of outer boundary center

compute R_inner inbound reduce ave v_r_coord # Radius of inner boundary center

compute F_rad_reac outbound reduce sum v_f_r # Radial reaction force summation on outer boundary grains (in units N)

variable P_rad_reac equal (c_F_rad_reac)/(2*PI*c_R_outer) # units N/m

compute F_theta_reac outbound reduce sum v_f_theta # Annular reaction force summation on outer boundary grains (in units N)

variable S equal (c_F_theta_reac*c_R_outer)/(2*PI*c_R_inner*c_R_inner) # units N/m

variable mu_wall equal v_S/v_P_rad_reac

variable gres equal sqrt(m*v_kn)*100 # gp as before
variable Vr equal (v_P_rad_reac-{Pwall})*(2*PI*c_R_outer)/${gres} # Radial velocity with which outer boundary grains should all move

# Convert radial velocity back to atom-wise x and y velocities

variable V_x atom v_Vr*cos(v_theta_coord)
variable V_y atom v_Vr*sin(v_theta_coord)

# Move the outer boundary

fix 3 outbound move variable NULL NULL NULL v_V_x v_V_y NULL

thermo_style custom step atoms ke vol c_1 c_R_inner c_R_outer c_F_rad_reac v_P_rad_reac v_Vr c_F_theta_reac v_S v_mu_wall # output equal/vector-style variables

thermo 100000

thermo_modify lost error norm no

compute_modify thermo_temp dynamic yes

# Output setting

dump 1 all custom 100000 output-shear-*.txt id type x y z vx vy vz fx fy fz v_r_coord v_theta_coord v_f_r v_f_theta v_V_x v_V_y diameter c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]

dump 2 all image 10000000 image.*.jpg type diameter &

zoom 1.8 size 600 1024

dump_modify 2 pad 5 backcolor white

dump 3 all movie 500000 movie.mpg type diameter &

zoom 1.8 size 600 1024

dump_modify 3 pad 5 backcolor white

# Set timestep

variable tauc equal (3.14^2+(ln($e))^2)^0.5*(m/4/v_kn)^0.5 # 7.9485^(-6) timestep (v_tauc*0.03)

# Let the system relax

run 0

# Start rotating the inner boundary

variable Vwall_tilde equal 0.01 # dimensionless wall velocity

variable Vwall equal {Vwall_tilde}*c_R_inner*sqrt({Pwall}/m) # dimensional wall velocity 367.4*0.01=3.674(m/s)
variable Twall equal (2*PI*c_R_inner)/{Vwall} # period of inner wall rotation (units s)

fix 4 inbound move rotate {x_orig} {y_orig} 0 0 0 1 ${Twall}

run 100000000

write_restart shearing-*.restart

Thanks,