Hello everyone,
Based on the LAMMPS manual(fix wall/gran/region command — LAMMPS documentation), (region command — LAMMPS documentation)
and a related discussion on the forum(Understanding fix controller),
I wrote a script that places four spheres (each with a 3 mm diameter) on a bottom plate at z = -5 mm.
After they settle, a “topcap” is created at z = 0.0 and moves downward at a velocity of -0.001 (m/s).
The simulation continues until the topcap reaches the bottom plate.
I use dump files to monitor the behavior of both the bottom plate and the topcap.
While the data for the bottom plate seems correct, the data for the topcap is not.
ITEM: TIMESTEP
9000000
ITEM: NUMBER OF ATOMS
4
ITEM: BOX BOUNDS ff ff ff
-7.4999999999999997e-03 7.4999999999999997e-03
-7.4999999999999997e-03 7.4999999999999997e-03
-7.4999999999999997e-03 1.0000000000000000e-02
ITEM: ATOMS f_bottom_plate[1] f_bottom_plate[2] f_bottom_plate[3] f_bottom_plate[4] f_bottom_plate[5] f_bottom_plate[6] f_bottom_plate[7] f_bottom_plate[8]
1 -3.51995e-194 -5.82521e-194 0.000277371 -0.00281015 -0.00309985 -0.005 0.0015
1 4.28305e-193 -3.54802e-193 0.000277371 0.00176391 -0.00158697 -0.005 0.0015
1 3.00496e-194 2.56651e-194 0.000277371 -0.0010733 -0.000612156 -0.005 0.0015
1 0 0 0.000277371 -0.00257804 0.00285807 -0.005 0.0015
ITEM: TIMESTEP
9000000
ITEM: NUMBER OF ATOMS
4
ITEM: BOX BOUNDS ff ff ff
-7.4999999999999997e-03 7.4999999999999997e-03
-7.4999999999999997e-03 7.4999999999999997e-03
-7.4999999999999997e-03 1.0000000000000000e-02
ITEM: ATOMS f_topcap[1] f_topcap[2] f_topcap[3] f_topcap[4] f_topcap[5] f_topcap[6] f_topcap[7] f_topcap[8]
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
Below is the minimal working example (MWE):
# 1. Setup simulation environments
units si
dimension 3
atom_style hybrid sphere
boundary f f f
newton off
processors * * *
# 2. Setup Communication and Neighboring list
comm_modify mode single vel yes
neigh_modify delay 0 every 1 check yes
# 3. Setup variables
variable critDt equal 1e-6
region domain_3D block -0.0075 0.0075 -0.0075 0.0075 -0.0075 0.01
create_box 2 domain_3D
# 4. Setup Walls
fix bottom_plate all wall/gran hertz/history 4.98E10 6.4E10 50.0 30.0 0.5 1 &
zplane -0.005 NULL contacts
fix xside_plate all wall/gran hertz/history 4.98E10 6.4E10 50.0 30.0 0.5 1 &
xplane -0.005 0.005 contacts
fix yside_plate all wall/gran hertz/history 4.98E10 6.4E10 50.0 30.0 0.5 1 &
yplane -0.005 0.005 contacts
# 5. Setup pair
pair_style gran/hertz/history 4.98E10 6.4E10 50.0 30.0 0.5 1
pair_coeff * *
# 6. Apply gravity and Eq of motion
fix grav_acc all gravity 9.81 vector 0.0 0.0 -1.0
fix viscous_damping all viscous 0.01
fix eqMotion all nve/sphere
# 7. Setup pouring
region gen_area block -0.003 0.003 -0.003 0.003 0.005 0.009
fix ins all pour 4 1 4767548 diam one 0.003 dens 2000 2000 region gen_area
# 8. Preparation for run
timestep ${critDt}
# 9. Setup dump
shell if [ -d "post_1" ]; then rm -rf post_1; fi # make directory for post
shell mkdir post_1
dump dump_atoms all vtk 10000 post_1/atoms*.vtk fx fy fz xu yu id type radius diameter x y z vx vy vz
# 10. Calculate forces acting on the base
compute computeBaseForceZ all reduce sum f_bottom_plate[4]
variable varBaseForceZ equal c_computeBaseForceZ
variable nStep equal step
fix printOutBaseZforce all print 1000000 "Step:${nStep}, BaseFoeceZ:${varBaseForceZ}"
shell if [ -d "dump_1" ]; then rm -rf dump_1; fi # make directory for post
shell mkdir dump_1
dump bottomPlateInfo all custom 1000000 dump_1/dump1*.contact f_bottom_plate[1] f_bottom_plate[2] f_bottom_plate[3] f_bottom_plate[4] &
f_bottom_plate[5] f_bottom_plate[6] f_bottom_plate[7] f_bottom_plate[8]
# 11. Pouring
run 4000000
# 12. Gen topcap and move
variable topcapVz equal -0.001
#region topcap_region plane 0.0 0.0 0.0 0.0 0.0 -1.0 move NULL NULL v_topcapVz
region topcap_block block -0.005 0.005 -0.005 0.005 0.0 0.001 move NULL NULL v_topcapVz side out
fix topcap all wall/gran/region hertz/history 4.98E10 6.4E10 50.0 30.0 0.5 1 region topcap_block contacts
shell if [ -d "dump_2" ]; then rm -rf dump_2; fi # make directory for post
shell mkdir dump_2
dump topcapInfo all custom 1000000 dump_2/dump2*.contact f_topcap[1] f_topcap[2] f_topcap[3] f_topcap[4] &
f_topcap[5] f_topcap[6] f_topcap[7] f_topcap[8]
unfix ins
#reset_timestep 0
run 9000000 upto # 0.001 m/s * 5s = -0.005 cm (Reach to the bottom wall), 5s /(1e-6) = 5000000 step
I have tried defining the wall as both a block region and a plane region, but neither approach appears to work correctly.
If you spot any mistakes or if I misunderstood anything in the manual, please let me know.
Thank you for your help!
Best regards,
- Information of LAMMPS build
$ ./lmp -help
Large-scale Atomic/Molecular Massively Parallel Simulator - 19 Nov 2024 - Development
Git info (develop / patch_19Nov2024-916-ga0fcbc9b71)
OS: Linux "Ubuntu 24.04.1 LTS" 6.8.0-52-generic x86_64
Compiler: Intel LLVM C++ 202500.0 / Intel(R) oneAPI DPC++/C++ Compiler 2025.0.4 (2025.0.4.20241205) with OpenMP 5.1
C++ standard: C++17
Embedded fmt library version: 10.2.0
MPI v3.1: Intel(R) MPI Library 2021.14 for Linux* OS
Accelerator configuration:
OPENMP package API: OpenMP
OPENMP package precision: double
OpenMP standard: OpenMP 5.1
INTEL package API: OpenMP
INTEL package precision: single mixed double
INTEL package SIMD: enabled
FFT information:
FFT precision = double
FFT engine = mpiFFT
FFT library = KISS
Installed packages:
GRANULAR INTEL MOLECULE OPENMP PYTHON RIGID VTK