#------------------------------------#
# Todo(次の初期化から読んでください) #
#------------------------------------#
# 1. 専門用語の統一
#   - 粒子と分子
# 2. 複数の粒子間相互作用をどのように反映させるか
# 3. 領域境界に作用するマクロ力学変数(主に力)をどのように出力するのか
#   - region単位での出力が可能? あるいはregionに接書している粒子を同定して、平均値から作用する力や応力を求める?
# 4. 専門用語の理解
#   - Ghost Atom、Processor、force cutoff、skin distance、neighbor cutoff distance
# 

#------------------------#
#         初期化         #
#------------------------#

# 単位系の決定
units si

# どのような分子をシミュレーションで使用するかを決定します
#   多分地盤材料で使用するのはsphereかbpm/sphereぐらいかな
#   付けるかどうかは任意ですが、シミュレーション中に分子の大きさを変えるかの設定ができます(0:不変, 1:可変)
atom_style sphere 0

# ニュートンの第3法則(作用反作用の法則)を満たすかどうか
#   接触する2粒子間の力に関する情報共有をするかどうかのフラグです。通信量が増大する代わりに、計算量をわずかに減らせます
newton on

#------------------------#
#       変数の定義       #
#------------------------#

# 杭挿入のために必要な変数
#   dz という名前の変数(variable)を用意して、
#   それをrunが読みだされている間にxからyまで比例的に変化させてください、という意味です
variable dz equal ramp(0,-0.3)

#------------------------#
#       領域の定義       #
#------------------------#

# region: シミュレーションのために必要な「領域」を定義するためのコマンド
#   定義した領域は例えば、粒で埋めたりとか、境界領域として使ったりとか、
#   領域内に存在する粒子をグループとして名づけられたりとか、消したりとか色々なことに使えます
# 例えば下のようなコマンドがあった場合、それぞれの文字列は次のような意味を持ちます
#     1     2     3     4   5     6    7  8 9    10      11
#   region REG0 block -0.4 0.4 -0.15 0.15 0 2 side in units box
#   1: コマンド名(これは変えられません)
#   2: 領域名(自分で命名する必要あり, アルファベットと数字のみ?)
#   3: スタイル(領域の形, blockの他にもconeとかcylinderとかいろいろある)
#   4~9: 領域のサイズ(x方向の最小値, x方向の最大値, ...みたいな並び順)
#   10: 3で決めたスタイル(立体)において内側を領域とするのか(side in)、外側を領域とするのか(side out)
#   11: 距離単位をどのように決定するのか。boxであればunitsコマンドで指定した単位が採用される。
#       latticeの場合にはこの前に格子間隔を別途定義する必要があります。
# 今回はREG0が振動台実験用の箱、REG1が挿入される杭に相当する領域になっています。
region REG0 block -0.4 0.4 -0.15 0.15 0 2 side in units box
region REG1 block -0.4 0.4 -0.15 0.15 0.45 0.5 side in units box

# create_box: シミュレーション用のboxの定義(←よい日本語訳が思いつきません...)
#   一つ上の領域は、単に幾何学的な立体形状に名前を付けただけです
#   ここではシミュレーション用のboxの定義とboxの中で使用する粒子の種類数と結合の種類数を記述します。
#   コマンドの後の最初の引数が粒子の種類数で、次に領域名がきます。
create_box 1 REG0

# change_box: 一つ上のコマンドで定義したboxのパラメータを追加で変更するコマンド
# 例えばボックスを時間を経て拡大させるとか、せん断ひずみを加えるとかもできます。
#   ただここでは境界条件の変更のために使います。
# それぞれの引数の意味は以下の通りです。
#   all: ここには本来分子のグループIDというものがきます。
#        グループIDは本来これよりも前もって定義されるものですが、
#        allとすることでbox内の全ての粒子に次の境界条件が適用されるようになります。
#   boudary f f f: x, y, z軸方向の境界条件を決定します。fは非周期かつ固定境界を意味します。(他のはドキュメントを参照してください)
change_box all boundary f f f

#------------------------#
#  粒子間相互作用の定義  #
#------------------------#

# pair_style: 1組の分子間(粒子間)に作用する相互作用を定義します。
# ここは既往文献と合わせてしっかり考えたい場所です！
#   高速化の都合上、以下は塑性せん断変形を考慮したフックの接触理論(gran/hooke/history)を用いています。
#   ただ、他の文献では違う相互作用を使用している場合があるので、用いている理論と使用しているパラメータについてざっと調査をお願いします。
# gran/hooke/historyの後に続く引数は弾性係数、粘弾性減衰係数、静的降伏基準です。(詳しくはドキュメントを参考にしてください)
# ここから下は未検証事項です
#   複数の粒子間相互作用を与えたい場合にはおそらくpair_styleとpair_coeffを繰り返しやらないといけない?
#   おそらくcreate_boxで2種類の粒子があるとした場合は次のような書き方になる?
#   create_box 2 REG0
#   pair_style gran/hooke/history 1000.0 NULL 1000.0 NULL 0.5 0
#   pair_coeff * *                                              (一旦全ての粒子の組み合わせ(1-1, 1-2, 2-2)について一つ上のgran/hooke/historyのパラメータを適用)
#   pair_coeff 1 2 2.0 2.0 1.0 1.0 1.0 1.0                      (1-2の組み合わせの6個のパラメータをそれぞれ2, 2, 1, 1, 1, 1倍する)
#   異なる相互作用モデルを用いても大丈夫なのか??
pair_style gran/hooke/history 1000.0 NULL 1000.0 NULL 0.5 0

# pair_coeff: 上で定義した相互作用を、異なる種類の粒子のどの組み合わせに適用させるかを指示します。
#   例えば2つの種類の粒子がある場合、2つの粒子の組み合わせは次の3通りです。(粒子1-粒子1, 粒子1-粒子2, 粒子2-粒子2)
#   本来であればこの組み合わせの一つ一つに対して、上で定義した相互作用モデルを当てはめていきます。
#       例えば次のようなコードは「粒子1-粒子1の組み合わせに対して、上で定義したgran/hooke/historyのパラメータを
#       2倍, 2倍, 1倍, 1倍, 1倍, 1倍したものを適用する」という意味です。
#       pair_coeff 1 1 2.0 2.0 1.0 1.0 1.0 1.0
#       (未検証事項) おそらくですが全てのパラメータが同じ場合にはパラメータの倍率の項は省略できます。
#   ただすべての組み合わせに対する相互作用を列挙するのは手間なので、*(アスタリスク)を使った一括適用ができます。
#       例えば"*"のみの場合には、全ての粒子を指します。
#       また"*3"などの場合には、1~3番目までの粒子を表すことができます。
pair_coeff	* *

# comm_style: 領域分割の手法を決定するコマンド
#   (ぶっちゃけよく分かってないです...)
#   LAMMPSは粒子の物理量のやり取りを行うために、Ghost Atom(お化け粒子)を用いています。
#   このお化け粒子を生成する際の領域分割手法について決定できます。
#   引数はbrickかtiledのどちらからです。
comm_style brick

# comm_modify: 上で定義した領域分割手法に関して、具体的に近傍の影響範囲やどのような情報をやり取りするかを決定するコマンドです
#   (こちらもぶっちゃけよく分かってないです...)
#   cutoffの引数を設定することで、設定された値以内に存在するお化け粒子の情報を各プロセッサは取得することが出来ます。
#   vel yesとすることで、粒子の速度情報も共有されることになります。
#       今回使用している相互作用モデルには粘弾性減衰の項の計算の際に、速度情報が必要なのでyesとなっています。
comm_modify vel yes

# neighbor: 近傍粒子リストを作成するためのcutoff距離と作成方法を定義するコマンド
#   (ぶっちゃけ分かってない...)
#   processorとは別に「近い」粒子同士の組み合わせをまとめたリストがあって、それをpairwise neighbor listと呼ぶ
#   
neighbor 0.1 bin

#------------------------#
#       解析の定義       #
#------------------------#

# fix: タイプステップあるいは最適化中に系に適用される操作を定義するコマンド
#   (簡単に言うとforループの中身的なもの)
#       引数の1つ目は自分でつけるfixの名前
#       引数の2つ目は粒子のグループID(allで全部の粒子グループを指す)
#       引数の3つ目は操作のコマンド名
#   fix fix-ID group-ID nve/sphereについて
#       nve/sphere: nveはconstant Number, Volume, Energyの略(小正準集団とも呼びます)
#       大きさがある粒子の一、速度、角速度の値を更新するために積分を実行します
#       (未検討事項)粒状体の位置更新の際にどのような手法があるのか, PFCは?
fix BOX1 all nve/sphere

# fix fix-ID group-ID gravity: 重力を与えるコマンド
#   (引数の説明は省略します...)
#   このコマンドによって指定する大きさ(入力例だと4つ目の9.81)はあくまで加速度で、力ではないことに注意してください。
#       質量に寄らない力を加えたい場合は fix addforceコマンドを使います
fix GRAV all gravity 9.81 vector 0 0 -1

# fix fix-ID group-ID wall/gran/region: 粒子と境界との条件を設定するコマンド
fix WALL0 all wall/gran/region hooke/history 1000.0 NULL 1000.0 NULL 0.5 0 region REG0
fix WALL1 all wall/gran/region hooke/history 1000.0 NULL 1000.0 NULL 0.5 0 region REG1

# fix INS all pour: 粒子を注ぐコマンド
fix INS all pour 2000 1 300719 vol 0.1 10 region REG1 diam one 0.035 dens 1 2


# compute: 一ステップごとに計算するパラメータの定義
compute     PRES0 all pressure TEMP0
compute     

# timestep: 解析時の時間間隔の決定
timestep    0.0001

# ログに出力する情報の設定
thermo_style	custom elapsed atoms
thermo		    1000
thermo_modify   lost ignore

# 何ステップ分解析を回すか
run     30000

# 粒子を注ぐのをやめる
unfix   INS

# 杭の挿入
region  REG2 cylinder z 0 0 0.1 0.35 0.8 side out units box move NULL NULL v_dz
fix     WALL2 all wall/gran/region hooke/history 100.0 100.0 1.0 1.0 0.5 0 region REG2

# 何ステップ分解析を回すか
run     20000