∞YJ】〕一 - 九州大学学術情報リポジトリ

n (2-10)

次に右辺第2項目について考える。ただし引は定義によりZ/の関数であるが、

Maclaurin展開に寄与しない。これはレベルOにおけるセルの大きさとl番下の階層を決めると、 akは一憲に決まる定数であるからである。

エ αk r=エOK(z-Zl-l)-k

⁺

ヱ

k

⁽

^{Z -Z/}

) ^- k _-

I ^_ (

^{Z/ -Z/ー1}

)

::i(Z-勾) k=l k = !|

+

司ノ-7ι

7ι

7'』ア'~ 'κ一 + 一つ Lκ一。 ∞Y白山一 +

(

^k⁺^111-1

ì (

^Z/

_-Z

/_I

) ^川

= エエak I l ^，� '

¹⁷¹

��: ^Á ) I (

Z-Z/_I

)

^川7 ^(2-11)

ここで式(2-10)の右辺について111=0のとき

') __ ^u ，:..:._k

^--，-

= u

+ __ ^u_2 ^{門 +�刊+} 4 A +

白 (

^{Z -Z/}

r

^{Z -Z/_I}

⁽

^{Z -Z/ー1}

t (

^{Z -Z/_I}

)

^・^，

(

^Z-Z/_I

)

^斗

m = 1のとき

子初k _、 (

^{Z/ - Z/ー1}

) _{_} _ニ

^al

(

^Z/-Z/_I

) ₊

¹ ^2a2

(

^Z/-Z/_I

)

1 3a

3(

^{Z/ -Z/_I}

)

白 (

-Z/ーl)k+1 (Z-Z/ーl)2 (z-Zl-l)3 (z-Zl-i)4

111=2のとき

モ予 - a

^1."";

(k +

'，--_---，_ ---'---;- 内

1)k (Z/ -Z/ー 1)2

^_=

al (Z/ -Z/ー1)2

' '; +

2(Z/-Z/ 3a -

I? _

�

+

_， ^I

a 3 .， .4.3 (z/-zト1)2

^' ^" ^"

+...

^，.

Z A 2(z-Zl-l)山 (Z -Z/ーlY (Z-Z/_I)斗 2(Z-Z/-IY 111=3のとき

モ(k+2)(k+1)k (z/-z/_1)3

_〉 _αι _内 ^_₌

al(z/-z/-lf

^'^" '

+

_，

， 4a2(z/-Z/- -

lf

_' _，.'

+...

，

f:1

凡 6

(Z-Z/-I)"十九 (Z-Z/_I)斗 (2-Z/_I)

となる。これらを(Z-Z/_I)についてまとめると ( k+m-1 (2/ -Z/ー1)川 ì

エヱ

^αk

l _k m l j(z-Zl-l) k + 111 z-zl-l (z-Zlー1) _{= +} 21G 2Ml(Zl-Zl-l)l

� r + 2a a 3 2(Z/ -Z/ー1)+ al (Z/ -ZI_I )2 1

(z-Zl-l) - l

1 4 l h +3向( Z/-Z/ー1)+3a2(Z/-Z/_1)2+al(Z/-Z/ーlf]+

(z-z/_lf l

= 之 ^l _l 三 _σ 川 ^一1 1 ド(ヤZ/一勾ι叫一斗iJl J )

11川h凶12=凶I戸凶=斗l

バ(

_一

z勾ιZ/_Iト川一斗l

)y _H

¹^y ^h制;

" \ k 一 1 ) (2・12)

となる。よって式(2-9)、(2-10)及び(2-12)から、レベルlにおけるセルの中心Z/

に関する多重極展開は、親のセルの中心Z/- 1に関する多重極展開 U(z) =-ILln _ε。 (z-zlーl)+-L エ

2πεOH1= l(z-Zl-l)fl (2-13)

(Z/ -Z/-I)川H1 fm -1 1

=αo n1 + E G Kl k -l j(Zl-zlー 1)川 (2・14)

ho =α。 (2-15)

に変換できる。図2-7に示すように、lつの親に4つの子が存在する。それぞ

25

れ子に関して、セルの中心を2/.i、多重極展開係数をαk

，iと定義する。

。k

，1

^α^k^，

4 Z/ � /".4

ak/

z/ーl

�.3

2/.2 z/)

図2-7

アップワードパス概念

「子における多重極展開」を

「親における多重極展開」に変換する

ここでi

=

1，2，3，4及びk = 0， 1， ^{• • •}， kmaxである。それぞれの子内の粒子集団が外部に作る拡張した空間電荷ポテンシャルを式(2-1 3)、 (2-14)及び(2-15)に変換すると

叫(z)=ームln(z-zl-l)+-LZ

_2πεo_，�I

₍

-Z/ー1)' b

^(2-16)

(zli-ZlJI 川 (171- 1 ì

+ 2

⁰^k

^il _k- ^l |(Zli-zl-ly

^(2-17a)

bO，i

^=αO.i ^(2-17b)

となる。よって親のセル内の粒子集団が、外部に作る拡張した空間電荷ポテンシャルは、式(2-16)をi

=

1，2，3，4に関して重ね合わせた

U(z)=-2LIn(z-Zl-l)+-LZ 2m0111=l(z-Zlー1)'11

^(2-18)

blll =ヱbl/l，i

^(2-19)

となる。式(2-18)はレベル1-1における lつのセルが作る拡張した空間電荷ポテンシャルだが、レベルlにおけるポテンシャル式(2-9)について多重極展開係数引をんに、セルの中心をZIからZI-lに変換したものと同じ式である。よって式(2-7)、(2-8)及び(2-9)から式(2-14)、(2-15)、(2-18)及び(2-19)への変換公式を用いると、レベル1-2に関するポテンシャル式も同様にして求まる。

2-1-5 ダウンワードノてス

zにおける拡張した空間電荷ポテンシャルを計算する場合を考える。ダウンワードパスでは、「各レベルについて、( ₁₎ Zの親に最隣接の親の子で、かつ(2 )

zの隣でない子」の多重極展開のみ計算する。

具体例として、図 2-8及び図 2-9 におけるz での拡張した空間電荷ポテンシヤルを計算する場合を考える。図2-8は、レベル 3において多重極展開が適用されるセルを示している。図2・9において、まずレベル0及びレベルlについては該当する子がないので計算しない。レベル2 について、丸印が示すセルの

多重極展開のみ計算する。レベル 3についても同様に丸印が示すセルの多重極展開のみ計算する。図2・8 と図2-9 で、多重極展開を適用するセルが重なるこ

とカfわかる。

図2-8 レベル3で多重極展開適用可能領域(斜線) 黒丸はテスト粒子

• レベルO

レベル2

• レベルl

レベル3

刻2-9 ダウンワードパス概念

各レベルについて斜線が引かれた領域に着いてのみ計算する。

2-1-6 拡張したSCEFの計算

これまでは拡張した空間電荷ポテンシャルをツリー法で計算するアルゴリズムに関して説明してきた。 PICで SCEFを計算する場合は、空間電荷ポテンシヤルを数値的に計算し、さらに数値的に微分するという2段階に分けて計算する。今回採用したツリー法の場合、式(2-2)に(2-18)を代入した式

ん I 1 ぶ 1nb

i (z) = -L + 一一芝

2_nE

o

Z - Z，ー1

2 πεo h l (z - zl- l )111+

^(2-20)

について計算すればよいので、 1段階で計算可能である。尚、ツリー法が適用できない領域に関しては、式(2・2)に(2-1 )を代入した式

ムル刊E

(2-21 )

を用いて直接計算すれば良い。式(2-20)と(2-21)の和を取り、(2-3)及び(2・4)に

代入することで、 SCEFが求まる。

2-1 -7 近接効果

直接計算の短所の lつに、近接効果が挙げられる。テスト粒子が一様に分布されているとき、計算結果は理論値に近づく。しかし乱数を用いて分布を決定する際には、粒子間隔が非常に短くなる場合がある。このときSCEFの大きさが急激に増加し、理論値に対して非常に大きな誤差を生じる。これを近接効果という[22]。ツリー法でも、多重極展開を適用できないセルに関しては直接計算を行うため近接効果が生じる。開発したシミュレーションコードでは、テス

ト粒子の位置から、ある半径の外にある粒子についてのみ直接計算を行うことで近接効果を回避した。この半径を近接半径と呼ぶ。近接半径を変えて、 SCEF を計算し、理論値からの相対誤差に関してヒストグラムを計算した。粒子集団

の分布はKV分布を仮定し、粒子数10万個、レベル7及びp=1で計算を行った。結果を図2-1 0に示す。

12000

10000

.F‘ロ，コd -.4 『

8000

，

、C司c』qE司..1 6000

業者

認トト 4000 2000

。。

図2-10 粒子数

1mmm fu・u・uauドハUζJハUハUハUつ-ハυ戸、dti44「、J

2 ₃ 4 5

相対誤差[0/0]

計算誤差の近接半径依存性 10万個

レベル 7 多重極次数 l

半径lcmのKVビーム

近接半径=225μmが誤差が小さい

近接半径が225μmのとき、最もヒストグラムの裾が狭まることがわかる。よ

ってシミュレーションでは、最適な近接半径を225μmとして計算を行った。

ヌ1 2-11は近接効果を回避した場合と回避しなかった場合のSCEFを理論値と比較したものである。近接効果を回避することにより、明らかに理論値からのずれを小さくできる。

o 近接半径無し o 近接半径 225凶n

? ， 1"" 1 0理論値

S 1.5 .1

〉2

� H回

ハU-EEA

!言、

忠臣� 0.5 長

nu ハUnu

2 4 6 8 ^.，EA ^ハU

r [mm]

災12-11 近嬢半径適用の効果 2-ト8 高速化テクニック

この節では、開発したシミュレーションコードに関して、高速化のために行ったテクニックを下に整理した。 Oは最も効果的だったテクニック、 0はOに

次ぎ効果的だったテクニックである。

。複素計算について、一般的に公開されているライブラリcomplex.hは使用せず、シミュレーション内で個々に計算する。

。ツリー法では、上述した計算過程をとるため多数のループ計算を必要とする。

このとき粒子を基準とするループではなく、セルを基準とするループを用いる。

具体例

for(n=O; n<Ntest; n++) { } ではなく

for(cell_x=O; cell_x<cell_x_max; cell_x++) { for(cell_y=O; cell_y<cell_y_max; cell_y++){}

とする。ここでcell_x_max及びcell_y_max はレベルの関数で、そのレベルにおけるセルの数である。

。直接計算に関して、作用・反作用を利用する。自分自身及び最近接のセル8 つの計9個に関して直接計算をする必要がある。ここで作用・反作用を利用する。粒子iが粒子2に与える電場をElとすると、粒子2 は粒子 lに-Elの電場を与える。このことを用いると5つのセルについてのみ直接計算を行えば良いことになる。

。複数個の割り算について、分子は異なるが分母が同じ場合、先に分母の逆数を計算する。

02で割る演算y = ^x^/2は、右シフト演算子を用いる。 y = x>> 2

o if(x==A) { if(y==B) { } if(y==C) { } if(yヱコD) {}

とするより

if(x==A && y==B) { } if(x==A && y==C) { } if(x==A && y==D) { }

とした方が速い場合がある。

2-2 ステップサイズ最適化

シミュレーションパラメータである「輸送行列で計算する際に重要な縦方向のステップサイズ」と、「粒子集団の密度分布を計算する場合に必要な横方向のステップサイズ」について最適化を行った。

2- 2-1 横方向ステップサイズの最適化

ここでは、 x方向の密度分布を求める場合を考える。テスト粒子集団の実空間上での位置を

(

^Xi'Yi

)

で表す。ただしi= 1，2，…，N，ωで、N，ωはテスト粒子数であるo xiの最大値をXmax、最小値をXminとする。横方向のステップ数をNbin、ステップサイズをムxとし、

盲ー冒

ムx= 山max ^�^�m1n Nbin

で定義するo X方向のテスト粒子数に関するヒストグラムをhN'lill (X)とするo hN'lill (X)は、 k�:S;X < (k ⁺l)�に入るXiを累積していくことにより求まる。ただし k = 0，1，・..，Nbin -1である。図2・12にNbinを変えたときのhN/JlII (x)を示す。 Nbinを変

えると、 hN!Jill (x) の形状が変わることがわかる。これは、 hN"il/ (x)が k� :S; X < (k + 1)む区間における密度分布の積分値であるためである。密度分布を求めるためにはhN/川(X)をk�三X^<(k + l)�区間で微分する必要がある。しかしhN/'il/ (x)はk�:S; x < (k ⁺1) �区間で1点与えられるだけなので微分することは不可能である。このためシミュレーションでは k�:S;x < (k + 1)�区間での

j川(x)の平均値

川一b

を計算し密度分布とした。図2-13 に図2- 12と同様のNbinで計算したときの I1N"ill(X)を示す。図2-12と比較すると、 I1N/1Î1I (x)のおおよその形状がNbinに依らないことがわかる。しかしNbinの大きさで、密度分布に次のような違いが見られる。 Nbinが小さいとき、 �が大きくなるので密度分布は角張った形となる。

逆にNbinが大きいときは、密度分布に多くのノイズが入る。このノイズは、粒子分布を求めるときに用いた乱数のばらつきに起因する。 '{骨らかな密度分布を求めるために、最適なNlバnを求める必要があるo r骨らかな密度分布を求めるために、最適なNbinを求める必要がある。このためNbinを変えて11Nl t(x)を計算し、

Nbinを変えることによる誤差

QN!Jill (X) =

{n川';/1

_(X)一凡-J

(け}2

を計算し、 QNhill(X)の最大値QN!Jill，maxを求めた。密度分布が角張った形をしている場合又はノイズが入っている場合は、 QN!Jill，maxが増加する。よってQN!Jill，maxが最小となるとき、密度分布は最も滑らかになる。図2-14にテスト粒子数を変えたときのQNJ Jnaxを示す。 Nbinが小さい場合、 QNtt，maxはテスト粒子数に依存

しないことがわかる。一方N仰が大きい場合、テスト粒子数が大きくなると乱数のばらつきが少なくなるので、 QNhill，maxは小さくなる。シミュレーションでは各粒子数について、ほほQN，ドmaxが最小となるNbù1 = 100近傍を採用した。

( ) i× コJz D E

0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00

-120 -80 -40 _。

一一一- N.₀₁₁₁. =400

一一一 NhIl=80

一一- Nb140

40 80 120

x [mm]

図2-12 N，を変えたときのヒストグラム

。tn

テスト粒子数 10000個、ガウス分布

一一一-Nbin=400

一一一- Nbin=80

-一一一Nbin=20

80 120

o 40

x [mm]

ー80 -40 0.00

-120 16.00 14.00 12.00

ハUハUハU'EEA

8.00 6.00 4.00 2.00

( K)口

一心ZC

Nbinを変えたときの密度分布

テスト粒千数 10000個、ガウス分布災12-13

ハUハUハU'BEA

Q，，"TL� のN

Nbin，max ^-^�' bin依存性

ハunuhunUハリハUハUハUハリハUハUハUハUハUハUハUハU11ぺノ-ζJIl--一一一一一一Sパ、J51、eじCじ

N N

Nb 100

災12-14 20

ハU咽冒E・E・ハUハU .• ，A

5 話EED

Zα

2-2-2 縦方向ステップサイズの最適化

軌道計算を行う際に用いる輸送行列及び非線形キックは縦方向のステップサイズ仰の関数であり、最適化が必要である。 �sを変えながら軌道計算を行い、リング1周後のrms ビームサイズを比較した。このとき、bare tune は(7.203，

5.229)、 KEK-PSで周囲するパンチビーム強度 1.4_X1012 particle per p ulse (ppp) のピーク強度を持つコースティングビームについて計算を行った。結果を図 2-15に示す。 �sがO.lm以下になると結果が飽和することがわかる。 �sが小さくなると計算回数が多くなるため、シミュレーションでは最適なステップサイズとして計算結果が飽和する最大の却、すなわちω=O.lm を採用した。

horizont

5.4 10.0

5.2 9.5

包Eぞ〈B 、_ぞE

� ω

N A 『 5.0 9.0

8.5 ^ヨ

4.8 8.0

出^】⁰^尽^口tL口i ⁴³

4.6 7.5 ^Vトhr-u⁴³

-4.4 7.0

4.2 6.5

0.001 0.01 0.1 10

�s [m]

図2-15 RMSビームサイズの�s依存性

2-3 シミュレーションコードのベンチマークテスト

PATRASHについて正当性を評価するために、同様の目的で開発された

ACCSIM [22]と、同じビーム及び加速器パラメータを用いて計算を行い、結果を比較した。これは 2000年3月に ACCSIMの作者であるF. _W. Jones氏が KEK-PSに来所したときに依頼し、実現したものである。結果を図 2-16に不す。 PATRASHの計算結果はACCSIMの結果とほぼ一致した。図2-16 に見ら

ドキュメント内九州大学学術情報リポジトリ (ページ 35-47)

∞YJ】〕一

エ αk r=エOK(z-Zl-l)-k

ヱ

k

(

) - k -

I _ (

)

::i(Z-勾) k=l k = !|

+

司ノ-7ι

7ι

7'』 ア'~ 'κ一 + 一 つ­ Lκ一 。 ∞Y白山一 +

(

ì (

-Z

) 川

= エエak I l ，� '

��: Á ) I (

)

') __ u ，:..:._k

= u

+ __ u_2 門 +�刊+ 4 A +

白 (

r

(

t (

)

(

)

子初k 、 (

) _ ニ

(

) +

(

)

3(

)

白 (

-Z/ーl)k+1 (Z-Z/ーl)2 (z-Zl-l)3 (z-Zl-i)4

111=2のとき

モ 予 - a

(k +

1)k (Z/ -Z/ー 1)2

al (Z/ -Z/ー1)2

' '; +

2(Z/-Z/ 3a -

I? _

+

a 3 .， .4.3 (z/-zト1)2

+...

Z A 2(z-Zl-l)山 (Z -Z/ーlY (Z-Z/_I)斗 2(Z-Z/-IY 111=3のとき

モ(k+2)(k+1)k (z/-z/_1)3

al(z/-z/-lf

+

， 4a2(z/-Z/- -

lf

+...

f:1

(Z-Z/-I)"十九 (Z-Z/_I)斗 (2-Z/_I)

となる。これらを(Z-Z/_I)についてまとめると ( k+m-1 (2/ -Z/ー1)川 ì

エ ヱ

l k m l j(z-Zl-l) k + 111 z-zl-l (z-Zlー1) = + 21G 2Ml(Zl-Zl-l)l

� r + 2a a 3 2(Z/ -Z/ー1)+ al (Z/ -ZI_I )2 1

(z-Zl-l) - l

1 4 l h +3向( Z/-Z/ー1)+3a2(Z/-Z/_1)2+al(Z/-Z/ーlf]+

(z-z/_lf l

= 之 l l 三 σ 川 一1 1 ド(ヤZ/一勾ι叫一斗iJl J )

バ(

一

)y H

" \ k 一 1 ) (2・12)

となる。よって式(2-9)、(2-10)及び(2-12)から、レベルlにおけるセルの中心Z/

に関する多重極展開は、親 のセルの中心Z/- 1に関する多重極展開 U(z) =-ILln ε。 (z-zlーl)+-L エ

2πεOH1= l(z-Zl-l)fl (2-13)

(Z/ -Z/-I)川H1 fm -1 1

=αo n1 + E G Kl k -l j(Zl-zlー 1)川 (2・14)

ho =α。 (2-15)

に変換できる。図2-7に示すように、lつの親に4つの子が存在する。それぞ

25

⁽

) ^- k _-

I ^_ (

7'』ア'~ 'κ一 + 一つ Lκ一。 ∞Y白山一 +

_-Z

) ^川

= エエak I l ^，� '

��: ^Á ) I (

') __ ^u ，:..:._k

+ __ ^u_2 ^{門 +�刊+} 4 A +

⁽

子初k _、 (

) _{_} _ニ

) ₊

モ予 - a

エヱ

l _k m l j(z-Zl-l) k + 111 z-zl-l (z-Zlー1) _{= +} 21G 2Ml(Zl-Zl-l)l

= 之 ^l _l 三 _σ 川 ^一1 1 ド(ヤZ/一勾ι叫一斗iJl J )

_一

)y _H

に関する多重極展開は、親のセルの中心Z/- 1に関する多重極展開 U(z) =-ILln _ε。 (z-zlーl)+-L エ

₍

^il _k- ^l |(Zli-zl-ly

i (z) = -L + 一一芝

盲ー冒