waterPMLコマンドのソースコードの概略

(2)位置の表現方法

(Outline of the source code of waterPML command; (2) Methods to represent a location)

◆概要 (Outline)

waterPMLコマンドでは速度場や応力場、格子セルの種類に関する情報、媒質の物性など 3次元空間内の様々な情報を配列を用いて表現する。配列要素番号が位置に対応する。食い違い格子を用いており、変数毎に定義される位置が異なるため、配列要素番号と位置の対応づけが複雑なものとなる。この対応付けを少しでも分かりやすく系統的なものにするため、 waterPMLコマンドでは位置を表す3種類の方法を併用している。
The waterPML command uses arrays to expresses various information in the 3-D space (e.g., the velocity and stress fields, the information on the classes of grid cells, and material properties). The index of each array corresponds to a location. The relation between the array index and location is complicated because of the staggered grid scheme in which different variables are defined at different locations. To make this complicated relation as complehensive and systematic as possible, the waterPML command uses the following three methods to express a location.

これらの方法の説明では 半格子点 という用語を用いる(著者の造語である)。格子セルの中心点と境界上の点をまとめたもの(図1)で、これらの点の間隔が格子セルサイズの半分になることから半格子点と呼んでいる。
The descriptions for these methods below uses a term half-grid nodes (an onomatopoiesis by the author). It points to the central and boundary points of grid cells (Fig. 1). They appear every half grid cell size; hence they are named as half-grid nodes.

図1. 半格子点の概念。黒枠で囲まれた茶色の正方形は個々の格子セルを表し、分かりやすさのため1つのセルを赤色で示した。青色の点が半格子点である。
Fig. 1. The concept of half-grid nodes. The brown squares bounded by black outlines are individual grid cells; one of them is shown by red for better illustration. The blue points are the half-grid nodes.

◆位置の3種類の表現方法 (Three methods to represent a location)

●方法1: 半格子点の方向別インデックス (Method 1: directional indices of half-grid nodes)

全ての半格子点に座標軸方向別に整数番号を割り当てたもので、それぞれの半格子点の位置は \(x\), \(y\), \(z\)座標を表す3つの整数の組で表現する(図2)。
Integer indices assigned to all half-grid nodes in each coordinate axis direction; the location of each half-grid node is expressed by a combination of three integers that represent the \(x\)-, \(y\)-, and \(z\)-coordinates (Fig. 2).

この方法の利点は番号と座標の間の換算が容易なことである。そのため設定の読み込みや計算結果の出力等、番号と座標の間の変換が必要になる場面で主に用いる。
An advantage of this method is the easiness of conversion between the indices and coordinates. Therefore, this method is mainly used in reading a configuration or outputting a result, where conversion between the indices and coordinates are needed.

プログラム全体を通して \(x\)座標に対応する番号を変数n0または配列要素n[0]、 \(y\)座標に対応する番号を変数n1または配列要素n[1]、 \(z\)座標に対応する番号を変数n2または配列要素n[2] で表す。
Throughout the program, the index for a \(x\)-coordinate is expressed by a variable n0 or an array component n[0], that for \(y\) by n1 or n[1], and \(z\) by n2 or n[2].

図2. 方向別インデックスの概念。
Fig. 2. The concept of directional indices.

●方法2: 半格子点の通し番号 (Method 2: a consecutive index of half-grid nodes)

全ての半格子点に一意な整数番号を割り当てたもので、 \((x,y,z)\)座標(3成分)を1つの整数で表現する(図3)。
Integer indices uniquely assigned to all half-grid nodes; the location of each half-grid node is expressed by a single integer that represents the three components of a \((x,y,z)\)-coordinate (Fig. 3).

この方法の利点は3次元空間内の全ての位置に関するループを 1つの制御変数で回せることである。そのため位置に関するループで主に用いる。
An advantage of this method is that it realizes a loop for all locations in the 3-D space with a single control variable. Therefore, this method is mainly used in loops for locations.

プログラム全体を通して半格子点の通し番号は変数igで表す。
Throughout the program, the consecutive index for half-grid nodes is expressed by a variable g0.

図3. 通し番号の概念。
Fig. 3. The concept of a consecutive index.

●方法3: 特定の量のための通し番号 (Method 3: a consecutive index for a specific quantity)

位置の関数として表される特定の1つの量(例: 速度の\(x\)成分) が定義される位置に対して一意な整数番号を割り当てたもので、 \((x,y,z)\)座標(3成分)を1つの整数で表現する(図4)。
Integer indices uniquely assigned to all definition points of a specific quantity that is a function of the place (e.g., the \(x\)-component of the velocity); the location of each point is expressed by a single integer that represents the three components of a \((x,y,z)\)-coordinate (Fig. 4).

この方法の利点はその量の定義点を1つの制御変数で無駄なく表現できることである。そのため位置の関数を表す配列の要素番号として主に用いる。
An advantage of this method is that it can express all definition points of the quantity with a single control variable without excess. Therefore, this method is mainly used as the index of an array that defines a function of location.

この方法による位置の表現で用いる変数には2種類のサブバージョンが存在する。
There are two subversions of variables to express the locations by this method.

1つ目(方法3-1) は特定の速度成分の定義点のみの通し番号igv、特定の応力成分の定義点のみの通し番号igtであり、これらはループ制御変数として用いられる。例えば「速度成分\(V_0^0\)の全ての定義点について以下の処理を実行せよ」というようなループでは変数igvを使用する。
One subversion (method 3-1) includes the consecutive index igv for the definition points of a velocity component, and the consecutive index igt for the definition points of a stress component. They are used as a control variable for loops. For example, an order “repeat the following processings for all definition points of a velocity component \(V_0^0\)” is implemented using the variable igv.
もう1つのサブバージョン(方法3-2) は位置を表す配列要素番号としての用途であり、これには「index_」から始まる変数名が用いられる。例えばindex_Vはある速度成分の位置に対応する配列要素番号を表す。
The other subversion (method 3-2) is used to express an array index that represents a location. For this purpose, variable names beginning with “index_” are used. For example, index_V is used for the array index that represents a location for a velocity component.
この違いは理解しづらいが、処理の順番の違いと思えば良い。まずループ制御変数としてigv, igt(方法3-1)を与え、ループ内で対応する半格子点の通し番号ig(方法2)を取得し、更にそのigを用いて各変数の配列要素番号index_∗(方法3-2)を取得するという処理の流れになる。つまりigv, igt → ig → index_∗ という一方向で処理が行われる。 igvとindex_∗の両方が必要なのは速度定義点に関するループ内でも応力など他の変数が必要になるためである。
This difference can be understood in a processing flow. First, igv or igt (method 3-1) is given as a control variable for a loop. Within the loop, the consecutive index for the corresponding half-grid node ig (method 2) is obtained. Then using ig, the array indices index_∗ (method 3-2) for various variables are obtained. Therefore the processing flow is always oriented from igv or igt through ig to index_∗. Both igv and index_∗ are needed because stress and some other variables are needed in a loop with regard to the velocity definition points.

図4. 特定の量のための通し番号の概念。 (a) 速度の水平成分用の例、(b)応力の非対角成分用の例。
Fig. 4. The concept of a consecutive index for a specific quantity: examples for (a) the horizontal component of the velocity, (b) the off-diagonal components of the stress.

◆3種類の方法で表現した整数間の変換 (Conversion of integers among the three representation methods)

●方法1→2 (Method 1 → 2)

半格子点の方向別インデックス(方法1)を同じ地点を表す半格子点の通し番号(方法2)に変換するには関数形式マクロ waterPML_n2ig を用いる。
Use a function-type macro waterPML_n2ig for the conversion from the directional indices of a half-grid node (method 1) to the consecutive index at the same node (method 2).

●方法2→1 (Method 2 → 1)

半格子点の通し番号(方法2)を同じ地点を表す半格子点の方向別インデックス(方法1)に変換するには関数形式マクロ waterPML_ig2n0, waterPML_ig2n1, waterPML_ig2n2 を用いる。
Use function-type macros waterPML_ig2n0, waterPML_ig2n1, and waterPML_ig2n2 for the conversion from the consecutive index of a half-grid node (method 2) to the unidirectional indices at the same node (method 1).

●方法2→3 (Method 2 → 3)

半格子点の通し番号(方法2)を同じ地点における特定の量のための通し番号(配列用; 方法3-2)に変換するには struct grid型構造体を用いる。
Use a struct grid-type structure for the conversion from the consecutive index of a half-grid node (method 2) to the consecutive index for a specific quantity at the same location (for arrays; method 3-2).

waterPMLコマンドのソースコード全体を用いて 1つのstruct grid型構造体を用いており、その変数名をstgとしている。また、構造体のメンバは具体的に以下のように与えている。
The waterPML command uses a single struct grid-type structure of a name stg throughout the entire code. Members of this structure are defined as follows.

メンバ Member	waterPMLコマンドにおける意味 Meaning in the waterPML command
Nv	異なる場所で定義された量の数であり、 \(12+3\times\)地震波動ソースの個数となる。 The number of quantities defined at different locations and equals to \(12+3\times\)the number of seismic sources.
Ng	半格子点の総数。 The total number of half-grid nodes.
index[iv][ig]	ig番目の半格子点(方法2)における iv番目の量のための通し番号(方法3-2)。 The consecutive index of the ivth quantity (method 3-2) at the igth half-grid node (method 2).

ここでivは変数番号と呼ばれ、具体的な定義は後述する。
Here, iv is called an index of variable, whose definition is described later.

●方法3→2 (Method 3 → 2)

特定の量のための通し番号(ループ制御変数; 方法3-1)を同じ地点における半格子点の通し番号(方法2)に変換するには int型2次元配列igv2ig, igt2igを用いる。
Use int-type 2-D arrays igv2ig and igt2ig for the conversion from the consecutive index for a specific quantity (for loop; method 3-1) to the consecutive index at the same half-grid node (method 2).

igv2ig:
第1引数を速度成分番号i(\(0=x\), \(1=y\), \(2=z\))、第2引数を特定の量のための通し番号igv(方法3-1)とする2次元配列で、 igv2ig[i][igv]の値はその地点における半格子点の通し番号igとなる。
A 2-D array with the 1st argument being the velocity component index i (\(0=x\), \(1=y\), \(2=z\)) and the 2nd argument being the consecutive index igv for a specific variable (method 3-1). The value of igv2ig[i][igv] equals the consecutive index ig at the same half-grid node.

その配列要素数を並べたint型1次元配列Nigvが別途定義される。 Nigv[i]はigv2ig[i]の配列要素数である。
To express the number of elements of this array, another int-type 1-D array Nigv is defined: Nigv[i] indicates the number of array elements of igv2ig[i].
igt2ig:
第1引数を応力成分番号iT (\(0=xx\), \(1=yy\), \(2=zz\), \(3=xy=yx\), \(4=yz=zy\), \(5=zx=xz\))、第2引数を特定の量のための通し番号igt(方法3-1)とする2次元配列で、 igt2ig[iT][igt]の値はその地点における半格子点の通し番号igとなる。
A 2-D array with the 1st argument being the stress component index iT (\(0=xx\), \(1=yy\), \(2=zz\), \(3=xy=yx\), \(4=yz=zy\), \(5=zx=xz\)) and the 2nd argument being the consecutive index igt for a specific variable (method 3-1). The value of igt2ig[iT][igt] equals the consecutive index ig at the same half-grid node.

その配列要素数を並べたint型1次元配列Nigtが別途定義される。 Nigt[iT]はigt2ig[iT]の配列要素数である。
To express the number of elements of this array, another int-type 1-D array Nigt is defined: Nigt[iT] indicates the number of array elements of igt2ig[iT].

◆変数番号 (Index of variable)

方法2→3の変換に登場した変数番号ivは下表のように定義される。同じ場所で定義される量に対しては同じivが用いられる。
The index of variable iv, which was introduced in the conversion from method 2 to 3, is defined in the Table below. Note that the same iv value is used for quantities defined at the same location.

量 Quantity	定義位置¹ Definition locations¹							iv
量 Quantity	c	x	y	z	xy	yz	zx	iv
密度\(\rho\) Density \(\rho\)	○	○	○	○	○	○	○	0
ラメ定数\(\lambda\) Lame constant \(\lambda\)
ラメ定数\(\mu\) Lame constant \(\mu\)
PML領域における波動の吸収の強さ\(\alpha^k\) (\(k=0,1,2\)) (計算式の(14)式) The strength of wave absorption in the PML volume \(\alpha^k\) (\(k=0,1,2\)) (Eq. 14 of formula)
差分方程式の係数\(c_0^k\) (\(k=0,1,2\)) ( 離散化後の計算式の(7)式) The coefficients \(c_0^k\) (\(k=0,1,2\)) of difference equations (Eq. 7 of discretized formula)
速度の\(x\)成分\(V_0^k\) (\(k=0,1,2\)) The \(x\)-component of velocity \(V_0^k\) (\(k=0,1,2\))		○						1
速度の\(y\)成分\(V_1^k\) (\(k=0,1,2\)) The \(y\)-component of velocity \(V_1^k\) (\(k=0,1,2\))			○					2
速度の\(z\)成分\(V_2^k\) (\(k=0,1,2\)) The \(z\)-component of velocity \(V_2^k\) (\(k=0,1,2\))				○				3
応力の\(xx\)成分\(\tau_{xx}^l\) (\(l=0,1,2\)) The \(xx\)-component of stress \(\tau_{xx}^l\) (\(l=0,1,2\))	○							4
応力の\(yy\)成分\(\tau_{yy}^l\) (\(l=0,1,2\)) The \(yy\)-component of stress \(\tau_{yy}^l\) (\(l=0,1,2\))
応力の\(zz\)成分\(\tau_{zz}^l\) (\(l=0,1,2\)) The \(zz\)-component of stress \(\tau_{zz}^l\) (\(l=0,1,2\))
応力の\(xy\)成分\(\tau_{xy}^l\) および\(yx\)成分\(\tau_{yx}^l\) (\(l=0,1,2\)) The \(xy\)- and \(yx\)-components of stress \(\tau_{xy}^l\), \(\tau_{yx}^l\) (\(l=0,1,2\))					○			5
応力の\(yz\)成分\(\tau_{yz}^l\) および\(zy\)成分\(\tau_{zy}^l\) (\(l=0,1,2\)) The \(yz\)- and \(zy\)-components of stress \(\tau_{yz}^l\), \(\tau_{zy}^l\) (\(l=0,1,2\))						○		6
応力の\(zx\)成分\(\tau_{zx}^l\) および\(xz\)成分\(\tau_{xz}^l\) (\(l=0,1,2\)) The \(zx\)- and \(xz\)-components of stress \(\tau_{zx}^l\), \(\tau_{xz}^l\) (\(l=0,1,2\))							○	7
観測点リスト Station list	○²							8
差分方程式の係数\(c_1^k\) (\(k=0,1,2\)) ( 離散化後の計算式の(8)式) The coefficients \(c_1^k\) (\(k=0,1,2\)) of difference equations (Eq. 8 of discretized formula)		○	○	○				9
差分方程式の係数\(c_2^{ijpl}\) (\(i,j,p,l=0,1,2; i=j, p=l\)) ( 離散化後の計算式の(14)式) The coefficients \(c_2^{ijpl}\) (\(i,j,p,l=0,1,2; i=j, p=l\)) of difference equations (Eq. 14 of discretized formula)	○							10
差分方程式の係数\(c_2^{ijpl}\) (\(i,j,p,l=0,1,2; i\neq j\)) ( 離散化後の計算式の(14)式) The coefficients \(c_2^{ijpl}\) (\(i,j,p,l=0,1,2; i\neq j\)) of difference equations (Eq. 14 of discretized formula)					○	○	○	11
is番目(\(is=0,1,2,\cdots\))の地震波動ソースによる等価体積力の\(x\)成分\(f_0\) The \(x\)-component of equivalent body force \(f_0\) from isth seismic source (\(is=0,1,2,\cdots\))		○³						\(12+3\times is\)
is番目(\(is=0,1,2,\cdots\))の地震波動ソースによる等価体積力の\(y\)成分\(f_1\) The \(y\)-component of equivalent body force \(f_1\) from isth seismic source (\(is=0,1,2,\cdots\))			○³					\(13+3\times is\)
is番目(\(is=0,1,2,\cdots\))の地震波動ソースによる等価体積力の\(z\)成分\(f_2\) The \(z\)-component of equivalent body force \(f_2\) from isth seismic source (\(is=0,1,2,\cdots\))				○³				\(14+3\times is\)

○が変数の定義位置を表す。スペースの都合上、以下のように略記した。
Circles indicate the definition points of the variable. To save the space, the following abbreviations are used.
- c:
  格子セル中心点(\(\posx_c\))。
  The centers of grid cells (\(\posx_c\)).
- x:
  \(x\)軸に直交する格子セル境界面の中心点 (\(\posx_c+\Delta\posx_0/2\))。
  The centers of the boundary planes of grid cells perpendicular to the \(x\)-axis (\(\posx_c+\Delta\posx_0/2\)).
- y:
  \(y\)軸に直交する格子セル境界面の中心点 (\(\posx_c+\Delta\posx_1/2\))。
  The centers of the boundary planes of grid cells perpendicular to the \(y\)-axis (\(\posx_c+\Delta\posx_1/2\)).
- z:
  \(z\)軸に直交する格子セル境界面の中心点 (\(\posx_c+\Delta\posx_2/2\))。
  The centers of the boundary planes of grid cells perpendicular to the \(z\)-axis (\(\posx_c+\Delta\posx_2/2\)).
- xy:
  \(x\), \(y\)軸に直交する格子セルの辺の中心点 (\(\posx_c+\Delta\posx_0/2+\Delta\posx_1/2\))。
  The centers of the boundary edges of grid cells perpendicular to the \(x\)- and \(y\)-axes (\(\posx_c+\Delta\posx_0/2+\Delta\posx_1/2\)).
- yz:
  \(y\), \(z\)軸に直交する格子セルの辺の中心点 (\(\posx_c+\Delta\posx_1/2+\Delta\posx_2/2\))。
  The centers of the boundary edges of grid cells perpendicular to the \(y\)- and \(z\)-axes (\(\posx_c+\Delta\posx_1/2+\Delta\posx_2/2\)).
- zx:
  \(z\), \(x\)軸に直交する格子セルの辺の中心点 (\(\posx_c+\Delta\posx_2/2+\Delta\posx_0/2\))。
  The centers of the boundary edges of grid cells perpendicular to the \(z\)- and \(x\)-axes (\(\posx_c+\Delta\posx_2/2+\Delta\posx_0/2\)).
観測点を含むセルのみ。
Only the cells that consist of a station.
等価体積力がノンゼロとなる地点のみ。
Only the locations where the equivalent body force is not zero.

戻る (Go back)
次へ (Next)
マニュアルのトップページに戻る (Go back to the top page of the documentation)