ATOM LOOP

여러 DEM프로그램들 (PFC, Yade 등)은 일반적인 스크립트 언

어 또는 파이썬을 지원하므로 Loop를 돌리는 것이 매우 간단한 작업이다. 특히, 시뮬레이션에 사용되는 Particle의 수와 특성을 인덱싱하여 뽑아내는 것이 매우 쉽다.

 

그러나, LIGGGHTS의 Script언어는 꽤나 독특한 형태이고, 사용되는 variable또한 equal style, atom style로 구분되며 각각 scalar, vector자료형을 가지므로 자료호환을 하려면 꽤나 고민을 많이 해야 된다.

 

아래는 오랜동안 고민하여 작성된 LIGGGHTS의 script로, 시뮬레이션에 사용된 atom의 개수를 찾아 loop를 돌려 id와 radius를 출력하는 script이다.

 

중요 "compute는 running 중 실행이 되어야만, run 종료 후 compute를 소환할 수 있다."

그래서, run 이전에 정의된 compute를 실행하기 위한 fix print과정이 필요하다.

compute                    cal_radius all property/atom update_on_run_end yes radius

variable                   imsi equal c_cal_radius[3]
fix                        print_test  all print 100 "${imsi}" screen yes  # This work for later use of compute

###
run                        1000 up to

variable                   N_atoms equal count(all)
print                      "N_atoms is = ${N_atoms}"

label                      loop_atoms
variable                   index_i loop ${N_atoms}
variable                   atom_rad equal c_cal_radius[${index_i}]
print                      "index and radius is ${index_i}, ${atom_rad}"
next                       index_i
jump                       v_test.in loop_atoms

 

아래는 Full script

 

### This is test file

### 1. Initializaion Phase

# Preliminaries
atom_style                 granular
atom_modify                map array                                  # necessary?
boundary                   f f f                                      # default is p p p
newton                     off                                        # need more information and examples
communicate                single vel yes
units                      si                                         # This command cannot be used after the simulation box
hard_particles             yes
dimension                  3                                          # default value is 3

# Neighbor listing
neighbor                   0.003 bin
neigh_modify               delay 0

# Declare domain
#region                     domain_1 block -0.04 0.04 -0.04 0.04 -0.03 0.08 units box   # xmin-max ymin-max zmin-max
region                     domain_1 block -0.1 0.1 -0.1 0.1 -0.03 0.08 units box   # xmin-max ymin-max zmin-max
create_box                 2 domain_1  # Declare number of atom types used in the simulation, This affect "fix property/global"

variable                   dt_ equal 0.00001

timestep                   ${dt_}
#timestep                   0.00001
fix                        grav_acc all gravity 9.81 vector 0.0 0.0 -1.0

### 2. Setup properties
# Material and interaction properties required
fix                         m_1 all property/global youngsModulus peratomtype 2.e7 2.e7
fix                         m_2 all property/global poissonsRatio peratomtype 0.25 0.25
fix                         m_3 all property/global coefficientRestitution peratomtypepair 2 0.1 0.1 0.1 0.1
fix                         m_4 all property/global coefficientFriction peratomtypepair 2 0.01 0.01 0.01 0.01

pair_style                 gran model hertz tangential history #Hertzian without cohesion
pair_coeff                 * *

### 3. Create atoms
create_atoms               1  single 0.0 0.0  -0.01   # atom #1
create_atoms               1  single 0.0 0.0   0.00   # atom #2
create_atoms               1  single 0.0 0.0   0.01   # atom #3
create_atoms               1  single 0.0 0.0   0.02   # atom #4
create_atoms               1  single 0.0 0.0   0.03   # atom #5
create_atoms               1  single 0.0 0.01 -0.01   # atom #6
create_atoms               1  single 0.0 0.01  0.00   # atom #6
create_atoms               1  single 0.0 0.01  0.01   # atom #6
create_atoms               1  single 0.0 0.01  0.02   # atom #6
create_atoms               1  single 0.0 0.01  0.03   # atom #6

set                        atom 1  type 1 density 2000 diameter 0.0005
set                        atom 2  type 1 density 2200 diameter 0.0006
set                        atom 3  type 1 density 2000 diameter 0.0007
set                        atom 4  type 1 density 2200 diameter 0.0008
set                        atom 5  type 1 density 2000 diameter 0.0009
set                        atom 6  type 1 density 2200 diameter 0.0010
set                        atom 7  type 1 density 2000 diameter 0.0011
set                        atom 8  type 1 density 2200 diameter 0.0012
set                        atom 9  type 1 density 2000 diameter 0.0013
set                        atom 10 type 1 density 2200 diameter 0.0014

### 5. Set compute
compute                    count_id all property/atom update_on_run_end yes id
compute                    cal_radius all property/atom update_on_run_end yes radius


variable                   ak equal c_count_id[3]
variable                   akk equal c_cal_radius[3]

### 4. Setup variables
variable                   step_now equal step
variable                   time_now equal step*v_dt_

### 4. Import Mesh from STL
fix                        box_base   all mesh/surface file base.stl        type 2 scale 1.0
fix                        box_side   all mesh/surface file side_wall.stl   type 2 scale 1.0
fix                        box_moving_L all mesh/surface file moving_wall_L.stl type 2 scale 1.0
fix                        box_moving_R all mesh/surface file moving_wall_R.stl type 2 scale 1.0

fix                        gran_walls all wall/gran model hertz tangential history mesh n_meshes 4 meshes box_base box_side box_moving_L box_moving_R

### 4. Set equation of motion
fix                        integral_eq_motion all nve/sphere
compute                    rke all erotate/sphere update_on_run_end yes

### 5. Create Imaging informations
run 1
dump                       dump_mesh_base     all mesh/vtk 100 post/base_mesh_*.vtk   id vel box_base
dump                       dump_mesh_side     all mesh/vtk 100 post/side_mesh_*.vtk   id vel box_side
dump                       dump_stl_moving_L  all mesh/stl 100 post/moving_L_stl_*.stl  box_moving_L
dump                       dump_stl_moving_R  all mesh/stl 100 post/moving_R_stl_*.stl  box_moving_R
dump                       dump_atoms all custom/vtk 100 post/atoms*.vtk id type x y z vx vy vz omegax omegay omegaz density diameter

### Setup print
fix                        print_time  all print 100 "current step, dt, time = ${step_now}, ${dt_}, ${time_now}, ${ak}, ${akk}" screen yes
#fix                        print_time  all print 100 "ak = ${time_now}, ${amp_1}, ${ped_1}, ${omega_1}, ${pos_1}" screen yes

dump                       dump_imsi          all custom 1000 test230727*.txt c_count_id
### 7. Run
run                        10000 upto

print                      "     "
print                      "===="
print                      "simulation end"
print                      "===="

### Variable test

variable                   N_atoms equal count(all)
print                      "N_atoms is = ${N_atoms}"

label                      loop_atoms
variable                   index_i loop ${N_atoms}
variable                   atom_rad equal c_cal_radius[${index_i}]
print                      "index and radius is ${index_i}, ${atom_rad}"
next                       index_i
jump                       v_test.in loop_atoms

====
simulation end
====
N_atoms is = 10
index and radius is 1, 0.00025
index and radius is 2, 0.0003
index and radius is 3, 0.00035
index and radius is 4, 0.0004
index and radius is 5, 0.00045
index and radius is 6, 0.0005
index and radius is 7, 0.00055
index and radius is 8, 0.0006
index and radius is 9, 0.00065
index and radius is 10, 0.0007