19. ちょっとしたツール
赤木 ：プラマーモデル作って時間積分だけではあんまり動きがないから、 もうちょっと動くのを作ってみない？
学生C：というと？
赤木 ：簡単なのはプラマーモデル2個の衝突ね。もっと沢山でもいいわ。
学生C：時間積分のプログラムに2つ読む機能とかつけますか？
赤木 ：そこは、処理速度としては問題あったりするけど、複数の単純な プログラムの組合せで実現する、っていう、まあ、Unix の精神でどうかしら？ つまり、
- スナップショットの位置、速度をずらすプログラム
-
複数のスナップショットをまとめて1つにするプログラム
学生C： とりあえずそれで2つとか3つぶつけるのができるわけですね。 やってみます。
ずらすのは、hackcode1 から読むとこと書くとこ残して、あとオプションで もってきたベクトルを位置とか速度に足すのか。えーと、、、こんなかな。
まずずらすほうだけ。こんな感じです。
1:require "clop"
2:require "./nacsio.cr"
3:include Math
4:include Nacsio
5:
6:optionstr= <<-END
7: Description: program to shift the position and velocity of a snapshot
8: Long description: program to shift the position and velocity of a snapshot
9:
10: Short name:-x
11: Long name: --shift-pos
12: Value type: float vector
13: Variable name:dx
14: Default value:0.0,0.0,0.0
15: Description:Shift value for position
16: Long description: Shift value for position
17:
18: Short name:-v
19: Long name: --shift-vel
20: Value type: float vector
21: Variable name:dv
22: Default value:0.0,0.0,0.0
23: Description:Shift value for velocity
24: Long description: Shift value for velocity
25:END
26:
27:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
28:options=CLOP.new(optionstr,ARGV)
29:
30:class Body
31: YAML.mapping(
32: pos: {type: Vector3, key: "r",default: [0.0,0.0,0.0].to_v,},
33: vel: {type: Vector3, key: "v",default: [0.0,0.0,0.0].to_v,},
34: )
35:end
36:
37:Nacsio.update_commandlog
38:Nacsio.repeat_on_snapshots(Body.from_yaml(""), true){|pp|
39: pp.each{|p|
40: p.p.pos += options.dx.to_v
41: p.p.vel += options.dv.to_v
42: p.print_particle
43: }
44:}
赤木 ：じゃあまずは説明ね。
学生C：はい。25行目までは基本的にオプションの文字列で、 -r, -v で位置、 速度のシフト量です。変数名は dx. dv です。
30行目からは粒子型の宣言で、今回位置と速度だけ使うので、それだけの 粒子にしてます。
37行目は update_commandlog で、標準入力からスナップショットファイルの ヘッダ部分読んで、自分のコマンド追加して書くものです。
38行目の repeat_on_snapshots は、nacsio.cr に追加した関数で、これ自体 は
1: def repeat_on_snapshots(p = Particle.new, read_all = false) 2: sp= CP(typeof(p)).read_particle 3: no_time = (sp.y["t"] == nil ) 4: while sp.y != nil 5: pp=[sp] 6: time = sp.y["t"].as_f unless read_all 7: while (sp= CP(typeof(p)).read_particle).y != nil && 8: (read_all || no_time || sp.y["t"].as_f == time) 9: pp.push sp 10: end 11: yield pp 12: end 13: endこういうものです。スナップショットを1個読んだら、というのは、時刻が変わった ら次とみなすとして、なんかする、という処理一般に使える関数があると便利、 ということで、 nacsplot2.cr とかであったループ構造を関数にしてます。
2行目でまずは粒子1つ読み込みます。それから、 一般的に使えるように、ということで、スナップショットの中に時刻データが なかったら、というのを3行目でチェックしてます。4行目の while は、粒子が読める限り読むものですが、ここでは読んでないことに 注意して下さい。ここでは、最初は sp には2行目で読んだものが入ってます。
while の下で配列 pp を作って、最初の要素をさっき読んだsp にして、6行目で時刻 も設定します。但し、引数の read_all が真なら時刻無視して全部読むので時 刻設定しません。7行目からの while で、粒子読んで、 時刻見る時は同じ時刻かどうかチェックします。で、チェックしないか、 同じ時刻なら pp に粒子追加します。 この while 抜けるのは、読んだ粒子の時刻が違うか、あるいは 全部読んだ時ですね。そこで、11行目の
11: yield ppが実行されます。この yield は、この関数呼び出し側の
Nacsio.repeat_on_snapshots(Body.from_yaml(""), true){|pp|
pp.each{|p|
p.p.pos += options.dx.to_v
p.p.vel += options.dv.to_v
p.print_particle
}
}
のブロックの {} の中に「置き換えられる」ということのようです。但し、
{} の中はこの関数の中というより呼び出し側の中で、呼び出し側の
変数とか使えるみたいです。 a が配列だとして、
s=1
a.each{|x| x +=s}
とか書くのと同じ感じですね。
repeat_on_snapshots は粒子と、 read_all の2つの引数ですが、
粒子は CP(T) 型を T に型いれて実体化したもので、これは実は
typeof(p) とするだけなので型だけあればよいので
Body.from_yaml("") を渡してます。 new でなくて from_yaml なのは、
Body クラスの中身が YAML.mapping だけなので、from_yaml しか
粒子作る方法がないからです。
で、 {|pp| と書いておくと、この pp に yield pp の中身が入って そのあとの } までが実行されて、それが終わると repeat_on_snapshots の次の行に戻るわけです。
で、処理本体は
pp.each{|p|
p.p.pos += options.dx.to_v
p.p.vel += options.dv.to_v
p.print_particle
}
で、粒子の、宣言した構造体になってる側で pos, vel に
それぞれコマンドラインオプションで指定された値を加えて、
そのあとで yaml 側に残っているデータと合わせて出力、をくり返します。
赤木 ：なるほど。プログラム本体は repeat_on_snapshots の中身が ほぼ全部ということね。わりと全体短くていい感じね。 実行結果は？
学生C：例えばこんな感じですね。
gravity> ./mkplummer -n 3 -Y -s 1 --- !CommandLog command: ./mkplummer -n 3 -Y -s 1 log: Plummer model created --- !Particle id: 0 t: 0.0 m: 0.3333333333333333 r: x: 0.6959348429044219 y: 0.5014594064614821 z: -0.05509228685899097 v: x: -0.07083562690974746 y: 0.369915247866214 z: -0.4293037912428686 --- !Particle id: 1 t: 0.0 m: 0.3333333333333333 r: x: -0.3818939282329637 y: -0.10402804359047872 z: 0.10433845010474857 v: x: 0.5625890492844403 y: 0.14162824345185604 z: 0.5622358681719519 --- !Particle id: 2 t: 0.0 m: 0.3333333333333333 r: x: -0.31404091467145817 y: -0.39743136287100317 z: -0.04924616324575757 v: x: -0.49175342237469283 y: -0.51154349131807 z: -0.13293207692908335 gravity> ./mkplummer -n 3 -Y -s 1|./nacsshift -x 2,3,4 --- !CommandLog command: './mkplummer -n 3 -Y -s 1 ./nacsshift -x 2,3,4' log: Plummer model created Unhandled exception: Expected sequence, not YAML::Nodes::Mapping at line 7, column 3 (YAML::ParseException) from /usr/share/crystal/src/yaml/nodes/nodes.cr:30:9 in 'raise' from /usr/share/crystal/src/yaml/from_yaml.cr:104:5 in 'new' from nacsio.cr:13:7 in 'initialize' from nacsio.cr:12:3 in 'new' from nacsshift.cr:45:3 in 'initialize' from nacsshift.cr:45:3 in 'new' from /usr/share/crystal/src/yaml/from_yaml.cr:12:3 in 'from_yaml' from nacsio.cr:108:30 in 'read_particle' from nacsio.cr:96:5 in 'read_particle' from nacsio.cr:149:5 in '__crystal_main' from /usr/share/crystal/src/crystal/main.cr:105:5 in 'main_user_code' from /usr/share/crystal/src/crystal/main.cr:91:7 in 'main' from /usr/share/crystal/src/crystal/main.cr:114:3 in 'main' from __libc_start_main from _start from ???座標の数字が 2,3,4 増えてるのがわかるかと。
赤木 ：なかなか素晴らしいわね。コマンドも何を実行したかが残るのね。
学生C：そうです。次は複数スナップショットですが、、、これ、 今の nacsio ではできないですよね？そもそもスナップショットを 標準入力からしか読まないじゃないですか？
赤木 ：そうねえ、まあ、だから、複数ファイルをつなげて標準入力から読ん で1つに、でもいいんだけど、一応複数ファイルを、でやってみて。 なので、 nacsio を色々な修正していいから。
学生C：わかりました。
といっても、えーと、どうするんだ？ read_particle はファイル引数にして、 今までのが動くようにするには STDIN がデフォルト値ならいいか。これは簡 単だな。
update_commandlog は、、、これは今、入力からログ読んで、書くまでやっ ちゃってるけど、そうじゃなくて読んでその値が入ってるインスタンス変数が 返るようにすればいいと。じゃあそういう関数を read_commandlog で作っ て、で、 read もファイルをもらって、今の配列に追加、、、するように もうなってるか。だから
c = CommandLog.new
fnames.each{|fn|
File.open(fn,"r"){|f|
c1 = read_commandlog(f)
read(body,ybody,f)
c.command += "\n"+c1.command
c.log += "\n"+c1.log
}
}
こんなので、粒子もコマンドログも追加するだけかな。
というわけで、、、なんかできたかな。
こんなふうです。
1:require "clop"
2:require "./nacsio.cr"
3:include Math
4:include Nacsio
5:
6:optionstr= <<-END
7: Description: program to add multiple snapshots to create one file
8: Long description: program to add multiple snapshots to create one file
9:
10: Short name:-i
11: Long name: --input-files
12: Value type: string
13: Variable name:fnames
14: Default value:none
15: Description:comma-separated list of input files
16: Long description: comma-separated list of input files
17:
18: Short name:-r
19: Long name: --reset-id
20: Value type: bool
21: Variable name:reset_id
22: Description:if true, reset id of particles
23: Long description: if true, reset id of particles
24:
25:END
26:
27:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
28:options=CLOP.new(optionstr,ARGV)
29:def write(body, ybody)
30: body.size.times{|i|
31: CP.new(body[i], ybody[i]).print_particle
32: }
33:end
34:def read(body,ybody, f)
35: while (sp= CP(typeof(body[0])).read_particle(f)).y != nil
36: body.push sp.p
37: ybody.push sp.y
38: end
39:end
40:
41:class Body
42: YAML.mapping(
43: id: {type: Int64, default: 0i64,},
44: )
45:end
46:
47:body = Array(Body).new
48:ybody= Array(YAML::Any).new
49:
50:c = CommandLog.new
51:options.fnames.split(",").each{|fn|
52: File.open(fn,"r"){|f|
53: c1 = read_commandlog(f)
54: c.command += "\n"+c1.command
55: c.log += "\n"+c1.log
56: read(body,ybody,f)
57: }
58:}
59:print c.add_command.to_nacs
60:body.each_with_index{|b,i| b.id=i.to_i64} if options.reset_id
61:write(body,ybody)
赤木 ：あら、割合いい感じ？では「解説お願いね。
学生C：まずオプションですが、 -i, -r の2つで、 -i でファイルネーム、 -r は、これいれると id を連続番号にします。
29行目からの write はさっきと同じですが、 34行目からの read は、 update_commandlog をなくして、入力ファイルをを引数にしてそこから読んで ます。
41行目からの Body は、 id だけですね。あとのデータさわらないので。
50行目で、空のコマンドログを作っておきます。
で、 body, ybody を初期化して、 fnames を配列にしてそれぞれから読むのが 51行目です、52行目の File.open(fn, "r") のあとに {|f| ...}がくるのは、
f= File.open(fn, "r") ...と変わらないですが、このブロック抜けたらファイル閉じるとかがあります。 で、53行目で コマンドログを読んで、54,55 行目で最初空のコマンドログに今読んだのを追加していきます。それから56行目で粒子を読んで配列に追加です。
ファイル読むのが終わったら、59 行目で、コマンドログに現在のコマンド名 を追加して出力です。で、オプションで指定があれば60行目でidふりなおして、 最後に write で粒子を出力でおしまいです。
赤木 ：2048粒子のプラマーモデルちょっと離したところからぶつけるとかしてみて。
学生C：コマンドはこんなのですかね。
./mkplummer -n 2048 -Y -s 1 > p2k.yml nacsshift < p2k.yml > tmp0 -x -3,0,0 -v 0.3,-0.2,0 nacsshift < p2k.yml > tmp1 -x 3,0,0 -v -0.3,0.2,0 nacsadd -r -i tmp0,tmp1 > 4k.yml env LD_LIBRARY_PATH=../../src/fdps-crystal/ ~/src/fdps-crystal/fdpscr -O 4k-out.yml -t 20 -d 1 -T 0.5 < 4k.yml env GKS_WSTYPE=jpg nacsplot3 < 4k-out.yml -w 10結果の絵はこんな感じです。
Figure 27: 初期条件
Figure 28: t=6
Figure 29: t=12
Figure 30: t=18
赤木 ：まあそれっぽいわよね。
この章は息抜きなので課題はありません。
特になし。この文書で作っているプログラムの
github 19.1. 課題
19.2. まとめ
19.3. 参考資料