14. Barnes-Hut treecode
赤木 ：前回、もうちょっと色々プロット、みたいな話をしたけど、
その前に時間積分できないと面白くないので、時間積分のプログラムつくりま
しょう。
学生C：前に作った particle1.cr でよくないですか？
赤木 ：あれはもちろんあれで動くんだけど、粒子数が多いと粒子数の自乗に
比例して時間かかるからあんまり大きな粒子数扱えないじゃない？なので、
今回は 11.4 で話がでてきた Barnes-Hut ツリー法を書いて
ね、と。
学生C：書いてね、いわれても、、、
赤木 ：まあこれも、作者が昔 Ruby で書いたのがあって、それは Josh
Barnes が C で書いたのが元になってるの。なので、それを Crystal にして
みて。で、ファイル入出力は PSDF で。
学生C：はあ、、、
赤木 ： Ruby のはこれ:
1:require "../command_line/clop.rb"
2:
3:def get_self_other_acc(myself, other, eps)
4: return @pos*0 if myself == other
5: rji = myself.pos - other.pos
6: r2 = eps * eps + rji * rji
7: myself.mass * rji / (r2 * sqrt(r2))
8:end
9:
10:class Body
11:
12: attr_accessor :mass, :pos, :vel, :acc
13:
14: def pairwise_acc(other, eps)
15: rji = other.pos - @pos
16: r2 = eps * eps + rji * rji
17: r = sqrt(r2)
18: r3 = r * r2
19: da = rji / r3
20: self.acc += other.mass * da
21: other.acc -= self.mass * da
22: end
23:
24: def ekin # kinetic energy
25: 0.5*@mass*(@vel*@vel)
26: end
27:
28: def epot(body_array, eps) # potential energy
29: p = 0
30: body_array.each do |b|
31: unless b == self
32: r = b.pos - @pos
33: p += -@mass*b.mass/sqrt(r*r + eps*eps)
34: end
35: end
36: p
37: end
38:
39: def get_node_acc(other, tol, eps)
40: get_self_other_acc(self, other, eps)
41: end
42:
43: def to_s
44: "mass = " + @mass.to_s +
45: " pos = " + @pos.join(", ") +
46: " vel = " + @vel.join(", ")
47: end
48:
49: def pp(indent = 0) # pretty print
50: print " "*indent + to_s + "\n"
51: end
52:
53: def ppx(body_array, eps) # pretty print, with extra information (acc)
54: STDERR.print to_s + " acc = " + @acc.join(", ") + "\n"
55: end
56:
57: def write
58: printf("%22.15e", @mass)
59: @pos.each do |x| printf("%23.15e", x) end
60: @vel.each do |x| printf("%23.15e", x) end
61: print "\n"
62: end
63:
64: def read
65: a = gets.split.collect{|x| x.to_f}
66: @mass = a[0]
67: @pos = a[1..3].to_v
68: @vel = a[4..6].to_v
69: end
70:
71:end
72:
73:class NBody
74:
75: attr_accessor :time, :body, :rootnode
76:
77: def initialize(n=0, time = 0.0)
78: @body = [Body.new]
79: for i in 0...n
80: @body[i] = Body.new
81: end
82: @time = time
83: end
84:
85: def evolve(tol, eps, dt, dt_dia, dt_out, dt_end, init_out, x_flag)
86: @dt = dt
87: @tol = tol
88: @eps = eps
89: @nsteps = 0
90: get_acc
91: e_init
92: write_diagnostics(x_flag)
93: t_dia = dt_dia - 0.5*dt
94: t_out = dt_out - 0.5*dt
95: t_end = dt_end - 0.5*dt
96:
97: write if init_out
98:
99: while @time < t_end
100: leapfrog
101: @time += @dt
102: @nsteps += 1
103: if @time >= t_dia
104: write_diagnostics(x_flag)
105: t_dia += dt_dia
106: end
107: if @time >= t_out
108: write
109: t_out += dt_out
110: end
111: end
112: end
113:
114: def leapfrog
115: @body.each do |b|
116: b.vel += b.acc*0.5*@dt
117: b.pos += b.vel*@dt
118: end
119:# get_acc
120: get_tree_acc
121: @body.each do |b|
122: b.vel += b.acc*0.5*@dt
123: end
124: end
125:
126: def get_acc
127: @body.each{|b|b.acc = b.pos*0}
128: i = 0
129: while (i < @body.size)
130: j = i+1
131: while (j < @body.size)
132:@body[i].pairwise_acc(@body[j], @eps)
133:j += 1
134: end
135: i += 1
136: end
137: end
138:
139: def get_tree_acc
140: maketree
141: @rootnode.center_of_mass
142:# @rootnode.pp(0)
143: @body.each{|b| b.acc = @rootnode.get_node_acc(b, @tol, @eps)}
144: end
145:
146: def ekin # kinetic energy
147: e = 0
148: @body.each{|b| e += b.ekin}
149: e
150: end
151:
152: def epot # potential energy
153: e = 0
154: @body.each{|b| e += b.epot(@body, @eps)}
155: e/2 # pairwise potentials were counted twice
156: end
157:
158: def e_init # initial total energy
159: @e0 = ekin + epot
160: end
161:
162: def write_diagnostics(x_flag)
163: etot = ekin + epot
164: STDERR.print <<END
165:at time t = #{sprintf("%g", time)}, after #{@nsteps} steps :
166: E_kin = #{sprintf("%.3g", ekin)} ,\
167: E_pot = #{sprintf("%.3g", epot)} ,\
168: E_tot = #{sprintf("%.3g", etot)}
169: E_tot - E_init = #{sprintf("%.3g", etot - @e0)}
170: (E_tot - E_init) / E_init = #{sprintf("%.3g", (etot - @e0)/@e0 )}
171:END
172: if x_flag
173: STDERR.print " for debugging purposes, here is the internal data ",
174: "representation:\n"
175: ppx
176: end
177: end
178:
179: def pp # pretty print
180: print " N = ", @body.size, "\n"
181: print " time = ", @time, "\n"
182: @body.each do |b| b.pp end
183: end
184:
185: def ppx # pretty print, with extra information (acc)
186: print " N = ", @body.size, "\n"
187: print " time = ", @time, "\n"
188: @body.each{|b| b.ppx(@body, @eps)}
189: end
190:
191: def write
192: print @body.size, "\n"
193: printf("%22.15e\n", @time)
194: @body.each do |b| b.write end
195: end
196:
197: def read
198: n = gets.to_i
199: @time = gets.to_f
200: for i in 0...n
201: @body[i] = Body.new
202: @body[i].read
203: end
204: end
205:
206: def maketree
207: @rootnode = makerootnode
208: @body.each do |b|
209: @rootnode.loadtree(b)
210: end
211: end
212:
213: def makerootnode
214: r = @body.inject(0){|oldmax, b| [oldmax, b.pos.map{|x| x.abs}.max].max}
215: s = 1
216: s *= 2 while r > s
217: Node.new([0.0, 0.0, 0.0], s)
218: end
219:
220:end
221:
222:class Node
223:
224: attr_accessor :mass, :pos
225:
226: def initialize(center, size)
227: @center, @size = center.to_v, size
228: @child = Array.new(8)
229: end
230:
231: def octant(pos)
232: result = 0
233: pos.each_index do |i|
234: result *= 2
235: result += 1 if pos[i] > @center[i]
236: end
237: result
238: end
239: def loadtree(b : Body)
240:# print "loadtree for center,size= #{@center}, #{@size}\n"
241: corner = octant(b.pos)
242:# print "new octant=#{corner}\n"
243: if @child[corner] == nil
244: @child[corner] = b
245:# p "inserted\n"
246: return
247: end
248: if @child[corner].class == Body
249: tmp_b = @child[corner]
250: child_size = @size / 2.0
251: @child[corner] = Node.new(@center + child_size*offset(corner),child_size)
252: @child[corner].loadtree(tmp_b)
253:# p "new cell made\n"
254: end
255:# p "recursive call"
256: @child[corner].loadtree(b)
257: end
258: def offset(corner)
259: r=[]
260: 3.times{ r.unshift( (corner & 1)*2 - 1 ) ; corner>>=1 }
261: r.to_v
262: end
263:
264: def pp(indent = 0)
265: print " "*indent+"node: center = #{@center.join(" ")} ; size = #{@size}\n"
266: if @mass
267: print " "*indent+" mass = #{@mass} pos = #{@pos.join(", ")}\n"
268: end
269: @child.each{|c| c.pp(indent + 2) if c}
270: end
271:
272: def check_body_in_cell
273: @child.each do |c|
274: if c.class == Body
275: (c.pos - @center).each do |x|
276: raise("\nbody out of cell:\n#{c.to_s}\n") if x.abs > @size
277: end
278: elsif c.class == Node
279: c.check_body_in_cell
280: end
281: end
282: end
283:
284: def center_of_mass
285: @mass = 0
286: @pos = [0, 0, 0].to_v
287: @child.each do |c|
288: c.center_of_mass if c.class == Node
289: if c
290: @mass += c.mass
291: @pos += c.mass * c.pos
292: end
293: end
294: @pos /= @mass
295: end
296:
297: def get_node_acc(b, tol, eps)
298: distance = b.pos - @pos
299: if 2 * @size > tol * sqrt(distance*distance)
300: acc = @pos*0
301: @child.each{|c| acc += c.get_node_acc(b, tol, eps) if c}
302: acc
303: else
304: get_self_other_acc(self, b, eps)
305: end
306: end
307:
308:end
309:
310:options_text = <<-END
311:
312: Description: First very simple version of Barnes-Hut tree code
313: Long description:
314: First very simple version of Barnes-Hut tree code
315:
316: (c) 2005, Piet Hut and Jun Makino; see ACS at www.artcompsi.org
317:
318: example:
319: ruby #{$0} -t 1 < cube1.in
320:
321:
322: Short name: -T
323: Long name:--opening_tolerance
324: Value type:float
325: Default value: 0.5
326: Global variable: tol
327: Description:Opening tolerance
328: Long description:
329: This option sets the tolerance value that governs the maximum size
330: of a tree cell that can remain closed; cells (nodes) with a size
331: large than the product of tolerance and distance to that cell will
332: be opened, and acceleration to its children will be computed.
333:
334:
335: Short name: -s
336: Long name:--softening_length
337: Value type:float
338: Default value: 0.0
339: Global variable: eps
340: Description:Softening length
341: Long description:
342: This option sets the softening length used to calculate the force
343: between two particles. The calculation scheme comforms to standard
344: Plummer softening, where rs2=r**2+eps**2 is used in place of r**2.
345:
346:
347: Short name: -c
348: Long name:--step_size
349: Value type:float
350: Default value:0.001
351: Global variable:dt
352: Description:Time step size
353: Long description:
354: This option sets the size of the time step, which is constant and
355: shared by all particles. It is wise to use option -s to specify a
356: softening length that is significantly larger than the time step size.
357:
358:
359: Short name: -d
360: Long name:--diagnostics_interval
361: Value type:float
362: Default value:1
363: Global variable:dt_dia
364: Description:Interval between diagnostics output
365: Long description:
366: The time interval between successive diagnostics output.
367: The diagnostics include the kinetic and potential energy,
368: and the absolute and relative drift of total energy, since
369: the beginning of the integration.
370: These diagnostics appear on the standard error stream.
371: For more diagnostics, try option "-x" or "--extra_diagnostics".
372:
373:
374: Short name: -o
375: Long name:--output_interval
376: Value type:float
377: Default value:1
378: Global variable:dt_out
379: Description:Time interval between snapshot output
380: Long description:
381: The time interval between output of a complete snapshot
382: A snapshot of an N-body system contains the values of the
383: mass, position, and velocity for each of the N particles.
384:
385: This information appears on the standard output stream,
386: currently in the following simple format (only numbers):
387:
388: N: number of particles
389: time: time
390: mass: mass of particle #1
391: position: x y z : vector components of position of particle #1
392: velocity: vx vy vz : vector components of velocity of particle #1
393: mass: mass of particle #2
394: ...: ...
395:
396: Example:
397:
398: 2
399: 0
400: 0.5
401: 7.3406783488452532e-02 2.1167291484119417e+00 -1.4097856092768946e+00
402: 3.1815484836541341e-02 2.7360312082526089e-01 2.4960049959942499e-02
403: 0.5
404: -7.3406783488452421e-02 -2.1167291484119413e+00 1.4097856092768946e+00
405: -3.1815484836541369e-02 -2.7360312082526095e-01 -2.4960049959942499e-02
406:
407:
408: Short name: -t
409: Long name:--duration
410: Value type:float
411: Default value:10
412: Global variable:dt_end
413: Description:Duration of the integration
414: Long description:
415: This option sets the duration t of the integration, the time period
416: after which the integration will halt. If the initial snapshot is
417: marked to be at time t_init, the integration will halt at time
418: t_final = t_init + t.
419:
420:
421: Short name:-i
422: Long name: --init_out
423: Value type: bool
424: Global variable: init_out
425: Description:Output the initial snapshot
426: Long description:
427: If this flag is set to true, the initial snapshot will be output
428: on the standard output channel, before integration is started.
429:
430:
431: Short name:-x
432: Long name: --extra_diagnostics
433: Value type: bool
434: Global variable:x_flag
435: Description:Extra diagnostics
436: Long description:
437: If this flag is set to true, the following extra diagnostics
438: will be printed:
439:
440: acceleration (for all integrators)
441:
442:
443: END
444:
445:parse_command_line(options_text)
446:
447:include Math
448:
449:nb = NBody.new
450:nb.read
451:nb.evolve($tol, $eps, $dt, $dt_dia, $dt_out, $dt_end, $init_out, $x_flag)
学生C：まあ、、、やってみます。
とりあえずみてみるか。class Body は、、、
attr_accessor :mass, :pos, :vel, :acc
はまた
YAML.mapping(
id: {type: Int64, default: 0i64,},
time: {type: Float64, key: "t", default: 0.0,},
mass: {type: Float64, key: "m",default: 0.0,},
pos: {type: Vector3, key: "r",default: [0.0,0.0,0.0].to_v,},
vel: {type: Vector3, key: "v",default: [0.0,0.0,0.0].to_v,},
acc: {type: Vector3, key: "a",default: [0.0,0.0,0.0].to_v,},
pot: {type: Float64, default: 0.0,},
)
に置き換えて、入出力もそっちになるとして、28行目の epot、これ tree 使
いようになってなくないか？ツリー使うって加速度計算するとこでポテンシャ
ルも計算しないと駄目だこれ。
49行目の pp とかその次の ppx とかは Crystal なら pp! ですむからいらな
いか。 57行目からの read/write も、nacsio の関数使うから不要と。
class NBody は、、、 initialize はいるのか？粒子配列は型だけ与えて長さ
0でいいんじゃないかなあ？どうせファイルから読むし。
そのあとの evolve, leapfrog はコンパイル通れば大丈夫そうかな。
get_acc これは使って、、、あ、ちゃんと get_tree_acc でポテンシャル計算
してこっち使わないと駄目か。なのでこれはいらないと。
get_tree_acc は、get_tree_acc_and_pot かなんかにして加速度とポテンシャ
ルと両方計算するようにすればいいかな。
あと、179行目の pp, ppx はやっぱり pp! でいいと。191行目からの
write/read は nacsio 使うからだいぶ変更必要だな。あとはまあコンパイラ
が文句いわなければ大丈夫かなあ、、、
ではまあやってみますか。
(数日後)
学生C：なんかできたみたいですが、、、実行すると
gravity> crystal mkplummer.cr -- -n 100 -s 1 -Y | crystal hackcode1.cr
[2mShowing last frame. Use --error-trace for full trace.[0m
In [4mlib/nacsio/src/nacsio.cr:21:10[0m
[2m 21 | [0m[1mYAML.mapping([0m
[32;1m^------[0m
[33;1mError: undefined method 'mapping' for YAML:Module[0m
[2mShowing last frame. Use --error-trace for full trace.[0m
In [4mlib/clop/src/clop.cr:25:8[0m
[2m 25 | [0m[1m{{system("sh ./lib/clop/src/clop_process.sh #{l} #{f} #{d} #{strname}")}}[0m
[32;1m^-----[0m
[33;1mError: error executing command: sh ./lib/clop/src/clop_process.sh 109 "/home/makino/papers/intro_crystal/hackcode1.cr" "/home/makino/papers/intro_crystal" "optionstr", got exit status 1[0m
で、エネルギーはまあまあ保存してます。
赤木 ：プログラムはどんな感じ？
学生C：以下ですが、
1:require "clop"
2:include Math
3:require "nacsio"
4:include Nacsio
5:
6:optionstr = <<-END
7: Description: First very simple version of Barnes-Hut tree code
8: Long description:
9: First very simple version of Barnes-Hut tree code
10: Crystal version - (c) 2020- Jun Makino
11: Original Ruby version -
12: (c) 2005, Piet Hut and Jun Makino. see ACS at www.artcompsi.org
13: example
14: hackcode1 < cube1.in
15:
16: Short name: -T
17: Long name:--opening_tolerance
18: Value type:float
19: Default value: 0.5
20: Variable name: tol
21: Description:Opening tolerance
22: Long description:
23: This option sets the tolerance value that governs the maximum size
24: of a tree cell that can remain closed; cells (nodes) with a size
25: large than the product of tolerance and distance to that cell will
26: be opened, and acceleration to its children will be computed.
27:
28: Short name: -s
29: Long name:--softening_length
30: Value type:float
31: Default value: 0.05
32: Variable name: eps
33: Description:Softening length
34: Long description:
35: This option sets the softening length used to calculate the force
36: between two particles. The calculation scheme comforms to standard
37: Plummer softening, where rs2=r**2+eps**2 is used in place of r**2.
38:
39: Short name: -c
40: Long name:--step_size
41: Value type:float
42: Default value:0.0078125
43: Variable name:dt
44: Description:Time step size
45: Long description:
46: This option sets the size of the time step, which is constant and
47: shared by all particles. It is wise to use option -s to specify a
48: softening length that is significantly larger than the time step size.
49:
50:
51: Short name: -d
52: Long name:--diagnostics_interval
53: Value type:float
54: Default value:0.25
55: Variable name:dt_dia
56: Description:Interval between diagnostics output
57: Long description:
58: The time interval between successive diagnostics output.
59: The diagnostics include the kinetic and potential energy,
60: and the absolute and relative drift of total energy, since
61: the beginning of the integration.
62: These diagnostics appear on the standard error stream.
63: For more diagnostics, try option "-x" or "--extra_diagnostics".
64:
65: Short name: -o
66: Long name:--output_interval
67: Value type:float
68: Default value:2
69: Variable name:dt_out
70: Description:Time interval between snapshot output
71: Long description:
72: The time interval between output of a complete snapshot
73: A snapshot of an N-body system contains the values of the
74: mass, position, and velocity for each of the N particles.
75:
76: Short name: -t
77: Long name:--duration
78: Value type:float
79: Default value:1
80: Variable name:dt_end
81: Description:Duration of the integration
82: Long description:
83: This option sets the duration t of the integration, the time period
84: after which the integration will halt. If the initial snapshot is
85: marked to be at time t_init, the integration will halt at time
86: t_final = t_init + t.
87:
88: Short name:-i
89: Long name: --init_out
90: Value type: bool
91: Variable name: init_out
92: Description:Output the initial snapshot
93: Long description:
94: If this flag is set to true, the initial snapshot will be output
95: on the standard output channel, before integration is started.
96:
97: Short name:-x
98: Long name: --extra_diagnostics
99: Value type: bool
100: Variable name:x_flag
101:
102: Description:Extra diagnostics
103: Long description:
104: If this flag is set to true, the following extra diagnostics
105: will be printed;
106: acceleration (for all integrators)
107:END
108:
109:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
110:options=CLOP.new(optionstr,ARGV)
111:
112:struct AccPot
113: property :a, :p
114: def initialize(a : Vector3 = Vector3.new, p : Float64 = 0.0)
115: @a=a; @p=p
116: end
117: def +(other : AccPot)
118: AccPot.new(@a+other.a, @p+other.p)
119: end
120:end
121:
122:class Body
123:
124: YAML.mapping(
125: id: {type: Int64, default: 0i64,},
126: time: {type: Float64, key: "t", default: 0.0,},
127: mass: {type: Float64, key: "m",default: 0.0,},
128: pos: {type: Vector3, key: "r",default: [0.0,0.0,0.0].to_v,},
129: vel: {type: Vector3, key: "v",default: [0.0,0.0,0.0].to_v,},
130: acc: {type: Vector3, key: "a",default: [0.0,0.0,0.0].to_v,},
131: pot: {type: Float64, default: 0.0,},
132: )
133:
134: def get_other_acc_and_pot( other, eps)
135: return AccPot.new if self == other
136: rji = @pos - other.pos
137: r2 = eps * eps + rji * rji
138: rinv = 1.0/sqrt(r2)
139: AccPot.new(@mass * rji*rinv*rinv*rinv, -@mass *rinv)
140: end
141:
142: def ekin # kinetic energy
143: 0.5*@mass*(@vel*@vel)
144: end
145:
146: def get_node_acc_and_pot(other, tol, eps)
147: get_other_acc_and_pot(other, eps)
148: end
149:
150: def loadtree(b) exit end
151:
152: def center_of_mass
153: @pos
154: end
155:
156:end
157:
158:class NBody
159: property :time, :body, :rootnode, :ybody
160:
161: def initialize(time = 0.0)
162: @body = Array(Body).new
163: @ybody= Array(YAML::Any).new
164: @time = time
165: @eps = 0.0
166: @e0=0.0
167: @dt=0.015625
168: @rootnode = Node.new([0.0,0.0,0.0].to_v, 1.0)
169: @tol = 0.0
170: @nsteps = 0
171: end
172:
173: def evolve(tol : Float64,
174: eps : Float64,
175: dt : Float64, dt_dia : Float64,
176: dt_out : Float64, dt_end : Float64,
177: init_out, x_flag)
178: @dt = dt
179: @tol = tol
180: @eps = eps
181: @nsteps = 0
182: get_tree_acc
183: e_init
184: write_diagnostics(x_flag)
185: t_dia = dt_dia - 0.5*dt
186: t_out = dt_out - 0.5*dt
187: t_end = dt_end - 0.5*dt
188: write if init_out
189: while @time < t_end
190: leapfrog
191: @time += @dt
192: @nsteps += 1
193: if @time >= t_dia
194: write_diagnostics(x_flag)
195: t_dia += dt_dia
196: end
197: if @time >= t_out
198: write
199: t_out += dt_out
200: end
201: end
202: end
203:
204: def leapfrog
205: @body.each do |b|
206: b.vel += b.acc*0.5*@dt
207: b.pos += b.vel*@dt
208: end
209: get_tree_acc
210: @body.each do |b|
211: b.vel += b.acc*0.5*@dt
212: end
213: end
214:
215: def get_tree_acc
216: maketree
217: @rootnode.center_of_mass
218: @body.each{|b|
219: # b.acc = @rootnode.get_node_acc(b, @tol, @eps)
220: ap= @rootnode.get_node_acc_and_pot(b, @tol, @eps)
221: b.acc=ap.a
222: b.pot = ap.p
223: }
224: end
225:
226: def ekin # kinetic energy
227: e = 0.0
228: @body.each{|b| e += b.ekin}
229: e
230: end
231:
232: def epot # potential energy
233: e = 0.0
234: @body.each{|b| e += b.pot*b.mass}
235: e/2 # pairwise potentials were counted twice
236: end
237:
238: def e_init # initial total energy
239: @e0 = ekin + epot
240: end
241:
242: def write_diagnostics(x_flag)
243: etot = ekin + epot
244: STDERR.print <<-EOF
245:at time t = #{sprintf("%g", time)}, after #{@nsteps} steps :
246: e_kin = #{sprintf("%.3g", ekin)}, \
247: e_pot = #{sprintf("%.3g", epot)}, \
248: e_tot = #{sprintf("%.3g", etot)}
249: e_tot - e_init = #{sprintf("%.3g", etot - @e0)}
250: (e_tot - e_init) / e_init = #{sprintf("%.3g", (etot - @e0)/@e0 )}\n
251:EOF
252:
253: if x_flag
254: STDERR.print " for debugging purposes, here is the internal data ",
255: "representation:\n"
256: pp! self
257: end
258: end
259:
260: def write
261: @body.size.times{|i|
262: @body[i].time = @time
263: CP.new(@body[i], @ybody[i]).print_particle
264: }
265: end
266: def read
267: update_commandlog
268: while (sp= CP(Body).read_particle).y != nil
269: @body.push sp.p
270: @ybody.push sp.y
271: end
272: @time = @body[0].time
273: end
274:
275: def makerootnode : Node
276: r = @body.reduce(0){|oldmax, b| [oldmax, b.pos.to_a.map{|x| x.abs}.max].max}
277: s = 1.0
278: while r > s
279: s *= 2
280: end
281: Node.new([0.0, 0.0, 0.0].to_v, s)
282: end
283: def maketree
284: @rootnode = self.makerootnode
285: i=0
286: @body.each do |b|
287:# print "loading body #{i} #{b.pos}\n"
288: @rootnode.loadtree(b)
289: i+=1
290: end
291: end
292:
293:end
294:
295:class Node
296: property :mass, :pos
297:
298: def initialize(center : Vector3, size : Float64)
299: @child = Array(Node|Body|Nil).new(8,nil)
300: @pos = Vector3.new
301: @mass=0.0
302: @center, @size = center, size
303: end
304:
305: def get_other_acc_and_pot( other, eps)
306: return AccPot.new if self == other
307: rji = @pos - other.pos
308: r2 = eps * eps + rji * rji
309: rinv = 1.0/sqrt(r2)
310: AccPot.new(@mass * rji*rinv*rinv*rinv, -@mass *rinv)
311: end
312:
313: def octant(pos)
314: result = 0
315: p=pos.to_a
316: c=@center.to_a
317: p.each_index do |i|
318: result *= 2
319: result += 1 if p[i] > c[i]
320: end
321: result
322: end
323:
324: def loadtree(b : Node)
325: end
326: def loadtree(b : Nil)
327: end
328: def loadtree(b : Body)
329: corner = octant(b.pos)
330: c=@child[corner]
331: unless c
332: @child[corner] = b
333: else
334: if @child[corner].is_a?(Body)
335: tmp_b = @child[corner]
336: child_size = @size / 2.0
337: c = Node.new(@center + child_size*offset(corner),child_size)
338: c.loadtree(tmp_b)
339: @child[corner]=c
340: end
341: c.loadtree(b)
342: end
343: end
344:
345: def offset(corner)
346: r=[] of Float64
347: 3.times{ r.unshift( ((corner & 1)*2 - 1 )+0.0) ; corner>>=1 }
348: r.to_v
349: end
350:
351: def check_body_in_cell
352: @child.each do |c|
353: if c.is_a?(Body)
354: (c.pos - @center).each do |x|
355: raise("\nbody out of cell:\n#{c.to_s}\n") if x.abs > @size
356: end
357: elsif c.is_a?(Node)
358: c.check_body_in_cell
359: end
360: end
361: end
362:
363: def center_of_mass
364: @mass = 0.0
365: @pos = [0.0, 0.0, 0.0].to_v
366: @child.each do |c|
367: if c
368: c.center_of_mass if c.is_a?(Node)
369: @mass += c.mass
370: @pos += c.mass * c.pos
371: end
372: end
373: @pos /= @mass
374: end
375:
376: def get_node_acc_and_pot(b, tol, eps)
377: distance = b.pos - @pos
378: if 2 * @size > tol * sqrt(distance*distance)
379: ap = AccPot.new
380: @child.each{|c| ap += c.get_node_acc_and_pot(b, tol, eps) if c }
381: ap
382: else
383: self.get_other_acc_and_pot(b, eps)
384: end
385: end
386:
387:end
388:
389:nb = NBody.new
390:nb.read
391:STDERR.print "after nb.read\n"
392:nb.evolve(options.tol, options.eps, options.dt,options.dt_dia,
393: options.dt_out, options.dt_end, options.init_out, options.x_flag)
基本的に Ruby のをコンパイルできるようにちょっと変えただけなので、あん
まり独自にどうというのはありません。
大きく違うのは、まず、 112行目あたりで、 AccPot という struct を作って、
加算も定義して、加速度とポテンシャルを1つの関数で返すことができてそれ
を加算もできるようにしました。この型があるとだとポテンシャル返すもので
も、ツリーアルゴリズムの基本である、あるノードからの相互作用は子ノード
からの相互作用の合計、というのがそのまま書けるので。
あと、Nacsio を使うので、 read が
def read
update_commandlog
while (sp= CP(Body).read_particle).y != nil
@body.push sp.p
@ybody.push sp.y
end
@time = @body[0].time
end
で、最初にコマンドログを読み書きしてから粒子を読んで、で、
このコードの中で使う粒子を @bodyに、 YAML のデータを @ybody にいれます。
出力は
def write
@body.size.times{|i|
@body[i].time = @time
CP.new(@body[i], @ybody[i]).print_particle
}
end
で、更新した @body と元の @ybody をあわせて出力です。まあ、
mkplummer.cr の出力だと余計なものはないので、これしてもしなくても同じ
ですが、、、
赤木 ：あと、これ、粒子動くところみたいわね。
学生C：あ、、、これ、今、出力は、最初に CommandLog がでますが、あとは
ずーっと粒子がでるだけで、どこで次の時刻になるかとかわからないんですが、、、
赤木 ：あ、それは、PSDF の考え方としてはそれでいいの。粒子が ID と時刻
もってるでしょ？だから、少なくとも概念としては、この出力は、各粒子につ
いて、4次元時空の中での軌跡を与えてるから。
学生C：うーん、概念としてなんか格好いいのはそうかもしれないですが、実
際に絵にするのはどうすればいいんですか？
赤木 ：そうね、、、メモリが一杯ある計算機なら、まず時間方向全部読み込
んで、それから分割、でもいいんだけど、粒子数ちょっと増えると破綻するわ
ね。出力が時間方向に1000枚あるとメモリが1000倍いるとかになるから。
だから、例えば、今読んでるのと違う(先の)時刻になったら中断、とかでいい
んじゃない？
学生C：中断、ってどうすればいいんですか？
赤木 ：まあ、単にプログラムの制御構造の中で処理なら、例えば今の
nacsplot.cr では
pp=Array(Particle).new
while (sp= CP(Particle).read_particle).y != nil
pp.push sp.p
end
としているけど、気分として
粒子読む
while まだ粒子がある (ここでは読まない)
粒子配列を粒子1つで新しく作る
whileまだ粒子がある && 読んだ粒子の時刻が他と同じ
粒子配列に追加
end
絵を書く
粒子配列を捨てる
end
みたいな感じかな。関数とかクロージャとか使って格好良く書いてもいいんだ
けど、まあここでは特にそんなことするほどでもないから。
学生C：あ、これならなんとか。
1:require "grlib"
2:require "clop"
3:require "./nacsio.cr"
4:include Math
5:include GR
6:include Nacsio
7:
8:optionstr= <<-END
9: Description: Plot program for multiple nacs snapshot
10: Long description: Plot program for multiple nacs snapshot
11:
12: Short name:-w
13: Long name: --window-size
14: Value type: float
15: Variable name:wsize
16: Default value:1.5
17: Description:Window size for plotting
18: Long description:
19: Window size for plotting orbit. Window is [-wsize, wsize] for both of
20: x and y coordinates
21:END
22:
23:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
24:options=CLOP.new(optionstr,ARGV)
25:update_commandlog
26:ENV["GKS_DOUBLE_BUF"]= "true"
27:
28:wsize=options.wsize
29:setwindow(-wsize, wsize,-wsize, wsize)
30:setcharheight(0.05)
31:setmarkertype(4)
32:setmarkersize(1)
33:sp= CP(Particle).read_particle
34:while sp.y != nil
35: pp=[sp.p]
36: time = pp[0].time
37: pp! time
38: while (sp= CP(Particle).read_particle).y != nil && sp.p.time == time
39: pp.push sp.p
40: end
41: clearws()
42: box
43: mathtex(0.5, 0.06, "x")
44: mathtex(0.06, 0.5, "y")
45: text(0.6,0.91,"t="+sprintf("%.3f",pp[0].time))
46: polymarker(pp.map{|p| p.pos[0]}, pp.map{|p| p.pos[1]})
47: updatews()
48:end
49:c=STDERR.gets
こんなのですかね？
赤木 ：まあ動いたらこれでいいんじゃない？
学生C：動くのは動きました。ただ、1万粒子くらいで結構遅いです。YAML っ
て格好いいですが、やっぱり遅くないですか？
赤木 ：そこはやっぱり問題なのよね、、、浮動小数点形式をバイナリフォー
マットじゃなくて人間が読める形で、というだけで結構遅いから。
まあ、実は、大きな計算の時には、どうせそんなにディスク容量がないから逆
に問題なかったりするんだけど。読み書きのとこを並列化すれば結構速くなる
し。
学生C：そんなものですか？
赤木 ：わりと。
14.1. ツリーコード解説
赤木 ：一応ここでツリーコードってどういうものかをコードみながら解説お
願い。
学生C：え、私よくわかってないですが、、、
赤木 ：まあ、わかってない人が読んでくほうが読者のためになるかもしれないし。
学生C：えー、そうですかね？まあやりますが。
上のプログラムで、4行目までは色々 requireとか include、110行目まではは
コマンドラインオプション関係ですね。そこまで飛ばします。
112行目からの struct AccPot は、加速度とポテンシャルをまとめて、加算を
定義したクラスです。これなしで、加速度やポテンシャルをタプルでまとめて
返り値にしてもいいのですが、加算があるほうがみやすいコードになるので
これ定義してます。使ってるのは 380行目です。
122行目からは Body クラスですね。これ Particle でないのは
元々の Joshua Barnes のプログラムから名前がきてるからですね。
124行目からは I/O 用の YAML.mapping でいつもの通りです。
134行目 からの get_other_acc_and_pot は、自分から相手への重力とそのポ
テンシャルです。実際の計算に使うものですね。
146行目の get_node_acc_and_pot は get_other_acc_and_pot を呼ぶだけです
が、引数が違います。これは、ツリーからの重力を計算する関数と同じ名前
にしているものです。その下にある loadtree は、実際には使わない
はずですが形式的にコンパイラが文句をいわないように作ってあるものです。
center_of_mass は粒子の重心なので位置そのものとなります。
赤木 ：まだこの辺ツリーアルゴリズム本体じゃない感じ？
学生C：そうですね。 get_node_acc_and_pot は最終的には使われますが。
158行目からの NBody にはいります。こっちもツリーアルゴリズムというより
は、N体時間積分のプログラム全体、というところです。
161行目からは initialize で、色々な内部変数に値をいれます。
173行目からの evolve が時間積分プログラム本体です。
184行目の write_diagnostics はエネルギー保存とか書く関数です。
189行目からが時間積分ループ本体で、基本的には leapfrog 関数を呼ぶだけ
で、あとは write_diagnostics とか 粒子を書く write を指定された
時間毎に呼んでます。
203行目からが leapfrog ですね。これは Ruby バージョンからそのままなの
で、今まで作った関数渡すとかのではなくて素朴に書いてあります。
get_tree_acc がツリー法で相互作用計算するアルゴリズムの実体です。
215行目からが get_tree_acc で、これは maketree でツリー構造を作り、
rootnode.center_of_mass でツリーの各ノードの重心を再帰的に計算します。
それから、218行目からの @body.each で全粒子へのツリーからの重力を計算
します。
あと細かい関数色々ありますが、ここで見ておく必要があるのは 283行目から
の maketree で
す。
まず、284行目の makerootnode ですが、これはもどって 275行目からです。
全粒子の座標絶対値の最大値をとって、それがはいるとようツリーのトップレ
ベルの箱の大きさを決めてます。
で、そのあとは 286行目からの @body.each で、ツリーに粒子を1個づつ
挿入するアルゴリズムを使ってます。
あとは、なので、その下の Node クラスの loadtree、center_of_mass、
get_node_acc_and_pot ですね。
loadtree は再帰的です。これは、8個の子セルのどこにいくべきか
を決めて、そこでが空なら単に子セルの代わりに粒子にするして、それでおし
まいです。
で、既に粒子だったら、新しいセルを作って、今までそこにあった
粒子を loadtree でいれて
そこがセルなら既に少なくとも1つ粒子が入ってるいるセルなので、
loadtree自身を呼び出します。
ここまでくると、必ず子セルが Node になっているので、
子セルに対して loadtree を呼びます。
赤木 ：これ短くて格好いいけどトリッキーよね、、、
学生C：まあこれ、元々の Joshua Barnes のプログラムは C で、ポインタが
粒子さしたりノードさしたりするんですよね、、、Crystal だとその辺
自分が何かというのは is_a? とかでわかるので、わりとまだ、、、
赤木 ：そうねえ、、、
学生C：あともうちょっとなので、 center_of_mass は再帰的に重心を計算す
るだけです。 sum とか reduce でもうちょっと綺麗に書ける気がしますが、、、
赤木 ：まあ動いてればね。
学生C：376行目からが get_node_acc_and_pot です。これは、
まず、距離と自分のサイズを比べて、サイズが大きいなら子セルからの力の合
計、そうでなければ自分の重心からの力、となります。後は389行目以降の
メイン部分で、粒子を nb.read de「「読んで、evolve を呼んで終わりです。
赤木 ：なんかかけあしだけど、ないよりはいいかしらね。
14.2. 課題
-
エネルギーの他、重心の位置、速度、系の全角運動量について、どの程度
保存しているかを確認し、それが適切な程度かどうかを検討せよ。
-
粒子数、時間刻み、ソフトニング、見込み角を系統的に変化させて、
エネルギー保存がどうなるかを検討せよ、エネルギーの指標には、
t=0 から 1 まで積分して、その間の最大値を用いよ。
14.3. まとめ
-
Barnes-Hut ツリーアルゴリズムの単純な実装を作った。
14.4. 参考資料
A hierarchical O(N log N) force-calculation algorithm