取り上げた問題

.PSfrag replacements

2. 取り上げた問題

本研究では，以下のような２次元平面内の正方領域

Ω = {(x, y)|0 < x < 1, 0 < y < 1}

におけるポアソン方程式のディレクレ型境界値問題を考える．

 



∂ ² u

∂x ² + ∂ ² u

∂y ² = 2x(x − 1) + 2y(y − 1) (x, y) ∈ Ω

u(x, y) = 0 (x, y) ∈ ∂Ω.

(1)

この問題の厳密解は，

u(x, y) = x(x − 1) × y(y − 1) (2)

であるが，正方領域

Ω

を

x

方向，

y

方向それぞれに間隔

h

の等間隔に分割し，

( x i , y j ) = ( ih, jh )， u ( x i , y j ) = u i,j

として，式

(1)

のポアソン方程式を以下のように離散化する．

u _i+1,j − 2u _i,j + u _i−1,j

h ² + u _i,j+1 − 2u _i,j + u _i,j−1

h ² = −f _i,j . (3)

ここで，−f

i,j = 2x i (x i − 1) + 2y j (y j − 1)

である．この離散化された式（連立一次方程式）を反復解法で解く場合，例えば

SOR

法を適用すると，

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k+1) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (4)

のようになる．ここで

ω

は加速係数，kはステップ数である．本研究では，この計算領域を図

1

のように分割した場合について，いくつかの反復解法について，その解法の性質や加速係数の依存性等について調べた．

Area2 Area3 Area1

j

i

図

1:

領域分割

Area1

Area2 A

B i

j

i j

図

2:

データ更新

3. 反復解法

3.1 JSOR

法，JSORω法

JSOR

法

¹⁾

は，分割された領域ごとに

SOR

法を用いて計算を行い，その後にデータを更新する手法である．図

1

のように分割された領域は各々の計算機上にメモリ確保されるため，他のタスクのデータを得るためにはデータを受け渡す必要がある．このとき図

2

の

Area1

の点

A

での計算には，Area2の点

B

のデータが必要となる．また，Area2の点

B

での計算には，Area1の点

A

のデータが必要となる．このため，

分割で割り当てられた配列とは余分に配列を確保する必要がある．また，余分に確保された配列は他のタスクにより計算が行われるため，反復の繰り返し毎にタスク間で境界部分のデータ通信を行い，更新する必要がある．JSOR法が計算を行ってからデータを更新するので，この順序では，k

+ 1

ステップ目の計算を行う際に，他のタスクの情報は

k

ステップ目のデータしか得ることができない．したがって，Area1 の点

A

のような

i = imax

での計算では，Area2の点

B

の計算値を用いることになるが，

k

ステップ目のデータしか必要としていないため，領域を分割しない

SOR

法の式

(4)

と同じ式で表すことができる．一方，Area2の点

B

のような

i = 1

での計算では，

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (5)

であり，u

i−1,j

に

k + 1

ステップ目ではなく

k

ステップ目のデータを用いていることがわかる．以上のように，領域を分割して

SOR

法を解く場合，タスク間でデータを受け渡さなければならないため，領域を分割した境界

i = 1

での計算は，古いデータ

(k

ステップ目のデータ)を三つ用いた計算となる．本論文では簡単のため，二分割の場合について計算を行ったが，それぞれの境界でのデータ更新と代入の違いによる収束性の影響を調べるため，境界の左側で用いる加速係数を

ω 1

右側を

ω 2

，内部領域では

ω

として収束回数が最も少なくなるような

(ω ₁ , ω ₂ , ω)

の組をもとめた．以下，この手法を

JSORω

法とする．

3.2 PSOR

法，PSORω法

PSOR

法は，新しい

k + 1

ステップ目のデータをより多く用いられるようにするため考案されたものであり，分割された計算領域の開始行列での計算が終了した直後に，データを転送するように手順を変更した手法である．図

3

に

PSOR

法の計算過程とデータ転送についての順序について概念的に示したものを示す．図中の

Area2

の点

B

や

Area3

の点

B’

のように，それぞれのタスクで

i = 1

の計算が行われた後，すぐにデータを

Area1

の点

C

あるいは

Area2

の点

C’

へ転送する．このため，Area1の点

A

や

Area2

の点

A’

のように

i = imax

の計算する時には，他のタスクの

i = 1

の

k + 1

ステップ目のデータを得ることができることになる．以下，JSOR法の時と同様に

PSOR

法の手法について示す．

A’

B

Area1 Area2 Area3

B’

A

i=1 i=imax i=1

i=imax

i

C’

C

図

3: PSOR

法

Area1 Area2

i=imax i=1

j j

i i

i,j i-1,j

i,j+1

i+1,j

i,j-1

A

B

図

4: i = imax

での計算

Area1 Area2

i=imax i=1

j j

i i

i,j i-1,j

i,j+1

i+1,j

i,j-1

B A

図

5: i = 1

での計算

図

4

に分割された領域の境界部分におけるデータ通信の様子を示した．Area1が着目している計算領域であり，Area2が他のタスクの計算領域である．点

A

のように計算領域の

i = imax

での計算は

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k+1) _i−1,j + u ^(k+1) _i,j−1 + u ^(k+1) _i+1,j + u ^(k) _i,j+1 ) (6)

である．Area2の点

B

は既に計算され，すぐにデータを転送することから，u

_i+1,j

に

k + 1

ステップ目のデータを用いることができる．よって，新しいデータ

(k + 1

ステップ目のデータ)を三つ用いることで，より精度の高い計算が行えると考えられる．図

5

は，Area2が着目している計算領域であり，Area1が他のタスクの計算領域の場合である．点

B

のように計算領域の

i = 1

での計算は

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (7)

であり，式

(5)

と同じ式となる．実際には

PSOR

法はその構成から，SOR法と同じ収束性を持つことが

Xie

等

²⁾

によって示されている．これを実際に数値計算により検証した．計算実験は二分割の場合について計算を行ったが，それぞれの境界でのデータ更新と代入の違いによる収束性の影響を調べるため，また

JSORω

法の場合と同様にタスク境界での加速係数を変化させたものについても

PSORω

法とよび，計算

を行った．

4. 計算結果と考察

表

1

に，配列数をいくつか変えた場合それぞれについて，JSOR法，JSORω法，PSOR法，PSORω法で計算した最適な加速係数

(の組)

の収束回数を示す．表中の数値は収束回数であり，括弧内は収束回数が最小になった場合の

ω

の値である．比較のため，ヤコビ法，領域分割していない

SOR

法についても計算結果を示した．JSOR法，JSORω法，PSOR法，PSORω法は，二分割された計算領域をそれぞれ２台の計算機で計算することになるので，計算に要する時間は半分になる．（ただし，これにデータ転送の時間が加わる．）

表

1:

収束回数比較

配列数ヤコビ法

SOR

法

JSOR

法

JSORω

法

PSOR

法

PSORω

法

200 × 200 204030 886(1.971) 1971(1.968) 1508(1.981) 885(1.971) 842(1.970) 500 × 500 1267928 2244(1.988) 5273(1.986) 3713(1.992) 2243(1.988) 2120(1.988) 1000 × 1000 5072240 4526(1.994) 11767(1.992) 7356(1.996) 4525(1.994) 4274(1.994)

この表から，ヤコビ法と領域分割をしないで

SOR

法を用いた場合を比較すると，すべての配列数の場合について，大幅に収束回数が少なくなることがわかる．ヤコビ法は完全に並列に計算を行うことができるが，収束性を考えると

SOR

法を用いた反復解法を用いることが有効となるであろうことが容易に推測される．また，JSOR法では，それらの収束回数の半分より，１台で行った

SOR

法の収束回数が少なくなっており，分割の影響で，実質的な計算時間も多くかかることがわかる．一方，JSORω法は

JSOR

法に比べ，反復回数が少なくなっている．これはデータ受け渡し境界部分での加速係数を小さく取り，古いデータの過大評価を押さえた効果によるものであり，計算に要する時間は，１台で行った

SOR

法よりも若干少なくなることが予想される．しかし，SOR法と比較すると反復回数そのものは依然大きく，古いデータを用いていること自体の収束性への影響が非常に大きいことは明らかである．

次に，PSOR法の結果についてみてみると，SOR法と

PSOR

法がほとんど同じ最適加速係数を持っており，またほとんど差のない収束回数であることがわかる．これは，Xie等

²⁾

の報告を指示するものである．ただし，数回の差が出ているのは，初期値の影響によるものと考えられる．PSOR

ω

法は

PSOR

法と比べ，多少ではあるが収束回数が少ない．このことは，境界の加速係数を変えることで，収束性を上げる可能性があることを示している．以上から，PSOR法，もしくは

PSORω

法を用いることにより，並列計算においても十分大きな収束速度を得ることが期待できることがわかった．

次に加速係数の依存性について，表

2

に，JSORω法と

PSORω

法で収束回数が最小になった時の加速係数の組を示す．加速係数

ω

については

0 . 001

刻み，ω

1

と

ω 2

については

0 . 01

刻みで変化させている．収束回数が同じになる加速係数の組がある場合には，その時の誤差の小さい方を最小の組として選んだ．

まず

JSOR ω

法の場合には，

ω

に比べて

ω 1

が小さい時に収束回数が最小になっている．これはデータ受け渡し境界の古いデータの過大評価を押さえることで，収束性を改善できることを裏付けるものである．

また，JSORω法の場合には，いずれの場合も

ω ₂

が小さい．このことから，ω

₁

を用いた計算結果の影響を受け，ω

2

も小さく取る必要があったことを示している．

一方，PSORω法の場合には，タスク境界で用いられる加速係数

ω ₁

，ω

₂

は，ωに近い値となっている．

これは，JSORω法の境界の最適加速係数よりも大きい．新しいデータを多く用いることにより，境界の計算に，より大きな加速係数を用いることができたと考えられる．また

200 × 200，500 × 500

の場合では，

ω，ω 1

，ω

2

がほとんど同じ値であるが，1000

× 1000

の場合には

ω

と

ω 1

，ω

2

での違いが生じる．これについては今後の課題とし，他の問題や格子点数の異なる計算を行うことで，違いが生じる要因を特定していきたい．

表

2:

最適加速係数

配列数

ω ω 1 ω 2

収束回数

200 × 200(JSOR ω

法)

1.981 1.81 1.79 1508 500 × 500(JSORω

法)

1.992 1.84 1.83 3713 1000 × 1000(JSOR ω

法)

1.996 1.85 1.85 7356 200 × 200(PSOR ω

法)

1.970 1.97 1.98 842 500 × 500(PSORω

法)

1.988 1.97 1.98 2120 1000 × 1000(PSORω

法)

1.994 1.95 1.94 4274

5. ３次元計算の場合の計算結果

並列計算が十分その威力を発揮できると思われる大規模計算について，以下のような領域

Ω = {(x, y, z)|0 <

x < 1, 0 < y < 1, 0 < z < 1}

における三次元のポアソン方程式の境界値問題

 



∂ ² u

∂x ² + ∂ ² u

∂y ² + ∂ ² u

∂z ² = −2 (x, y, z) ∈ Ω u(x, y, z) = 0 (x, y, z) ∈ ∂Ω

(8)

について同様の計算を行った．計算格子数は

200 ×200 ×200

であり，初期条件はすべての格子点で

u i,j,k = 0

とした．それぞれの方法で得られた結果について，最適の収束回数を表

3

に示す．JSORω法については，

２次元の場合の結果から，ω

= 1 . 980，ω 1 = 1 . 81，ω 2 = 1 . 81

として計算を行い，PSOR

ω

法については，

ω = 1.970，ω 1 = 1.98，ω 2 = 1.98

とした．

表

3:

三次元での計算結果

SOR

法

JSOR

法

JSORω

法

PSOR

法

PSORω

法

425 606 576 424 418

この表から，三次元計算で，収束判定法を変えた場合にも，SOR法と

PSOR

法がほぼ同じ収束をすることがわかる．またそれに比べ，JSOR法は収束回数が多い．また，JSORω法と

JSOR

法と比較すると，タスク境界での加速係数を小さく取ることにより，収束性が改善していることがわかる．さらに，PSORω法と

PSOR

法との比較から，タスク境界での加速係数を変えることで，収束性に多少の改善がみられた．以上のように，三次元計算の場合にも２次元の場合で得られたのと同様の結果が得られた．

ドキュメント内第7回環瀬戸内応用数理研究部会シンポジウム2003 (ページ 34-38)

.PSfrag replacements

2. 取り上げた問題

Ω = {(x, y)|0 < x < 1, 0 < y < 1}

 



∂ 2 u

∂x 2 + ∂ 2 u

∂y 2 = 2x(x − 1) + 2y(y − 1) (x, y) ∈ Ω

u(x, y) = 0 (x, y) ∈ ∂Ω.

(1)

u(x, y) = x(x − 1) × y(y − 1) (2)

Ω

x

y

h

( x i , y j ) = ( ih, jh )， u ( x i , y j ) = u i,j

(1)

u i+1,j − 2u i,j + u i−1,j

h 2 + u i,j+1 − 2u i,j + u i,j−1

h 2 = −f i,j . (3)

i,j = 2x i (x i − 1) + 2y j (y j − 1)

SOR

u (k+1) i,j = (1 − ω ) u (k) ij + ω

4 ( h 2 f ij + u (k+1) i−1,j + u (k+1) i,j−1 + u (k) i+1,j + u (k) i,j+1 ) (4)

ω

1

Area2 Area3 Area1

j

i

1:

Area1

Area2 A

B i

j

i j

2:

3. 反復解法

3.1 JSOR

JSOR

1)

SOR

1

2

Area1

A

B

B

A

+ 1

k

A

i = imax

B

k

SOR

(4)

B

i = 1

u (k+1) i,j = (1 − ω ) u (k) ij + ω

4 ( h 2 f ij + u (k) i−1,j + u (k+1) i,j−1 + u (k) i+1,j + u (k) i,j+1 ) (5)

i−1,j

k + 1

k

SOR

i = 1

(k

ω 1

ω 2

ω

(ω 1 , ω 2 , ω)

JSORω

3.2 PSOR

PSOR

k + 1

3

PSOR

Area2

B

Area3

∂ ² u

∂x ² + ∂ ² u

∂y ² = 2x(x − 1) + 2y(y − 1) (x, y) ∈ Ω

u _i+1,j − 2u _i,j + u _i−1,j

h ² + u _i,j+1 − 2u _i,j + u _i,j−1

h ² = −f _i,j . (3)

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k+1) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (4)

¹⁾

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (5)

(ω ₁ , ω ₂ , ω)

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k+1) _i−1,j + u ^(k+1) _i,j−1 + u ^(k+1) _i+1,j + u ^(k) _i,j+1 ) (6)

_i+1,j

u ^(k+1) _i,j = (1 − ω ) u ^(k) _ij + ω

4 ( h ² f ij + u ^(k) _i−1,j + u ^(k+1) _i,j−1 + u ^(k) _i+1,j + u ^(k) _i,j+1 ) (7)

²⁾