関数main(FDM_1D_vertical.c)の動作

(Behaviour of function main in FDM_1D_vertical.c)

関数mainではパラメータの読み込みの後、 速度場と応力場を初期化(全ての場所で0.0に設定)し、時間発展の計算を行う。 一番の骨格となるのはコード中の下記の部分である。
Function main first reads parameters, then initialize the velocity and stress fields (as 0.0 everywhere), and then calculate the temporal evolution. The most essential portion of the code is shown below.
for(k=1;k<parameters.Nt;k++){
    for(i=0;i<parameters.Nz-1;i++){
         T[i]+=V2tau_coefficients[i]∗(V[i+1]-V[i]);
    }
    for(i=1;i<parameters.Nz;i++){
         V[i]+=tau2V_coefficient∗(T[i]-T[i-1]);
    }
    V[0]=waveform_use.value[k];
    waveform_output.value[k]=V[parameters.Nz-1];
 }

ここで\(k\)は時刻番号、\(i\)は\(z\)方向の格子点番号を表す。 また配列Vが速度場、配列Tが応力場を表している。 これらは位置\(z\)と時間の関数であるが、 複数の時刻での値を同時に保持する必要はないので\(z\)のみの関数として 1次元配列で与えている。
Here, \(k\) represents an index of time, and \(i\) does an index of grid node in the \(z\)-direction. Arrays V and T represent velocity and stress fields. Although they are functions of location \(z\) and time. However, the values of velocity and stress at different times need not be preserved simultaneously. Therefore the values are expressed as a function of \(z\) alone, in a form of 1-D arrays.
上記のコードで用いている時間発展方程式は 離散化で述べた (5)(6)式、あるいは(9)(10)式である。 ある時刻での速度場を用いてその\(\Delta t/2\)だけ後の時刻での応力場を求める式は (6)(10)式であり、以下の形になっている。
\[\begin{eqnarray} \tau’(z_o+(i+1/2)\Delta z,(k+1/2)\Delta t) &=& \tau’(z_o+(i+1/2)\Delta z,(k-1/2)\Delta t) \nonumber \\ & & +\frac{\Delta t}{\Delta z}V_s(z_o+(i+1/2)\Delta z)^2 \nonumber \\ & & \left[ V(z_o+(i+1)\Delta z,k\Delta t) \right. \nonumber \\ & & \left. -V(z_o+i\Delta z,k\Delta t) \right] \label{eq.tau.evolution} \end{eqnarray}\] また、ある時刻での応力場を用いてその\(\Delta t/2\)だけ後の時刻の速度場を求める式は (5)(9)式であり、以下の形になっている。
\[\begin{eqnarray} V(z_o+i\Delta z,(k+1)\Delta t) &=& V(z_o+i\Delta z,k\Delta t) \nonumber \\ & & +\frac{\Delta t}{\Delta z}\left[ \tau’(z_o+(i+1/2)\Delta z,(k+1/2)\Delta t) \right. \nonumber \\ & & \left. -\tau’(z_o+(i-1/2)\Delta z,(k+1/2)\Delta t) \right] \label{eq.V.evolution} \end{eqnarray}\] これらの時間発展方程式における物理量とプログラム中の変数の対応関係は 下表のようになる。
The code above uses temporal evolution equations of eqs. (5) and (6), or eqs. (9) and (10), described in Discretization. Equations to calculate the stress field at a time \(\Delta t/2\) later from the velocity field at a given time are expressed by eqs. (6) and (10) in that document, which have the form of eq. (\ref{eq.tau.evolution}) above. Equations to calculate the velocity field at a time \(\Delta t/2\) later from the stress field at a given time are expressed by eqs. (5) and (9) in that document, which have the form of eq. (\ref{eq.V.evolution}) above. Relations between the physical quantities in these temporal evolution equations and the variables in the program are summarized in the table below.
プログラム中の変数
Variable in the program 		物理量
Physical quantity
V[i] 位置\(z=z_o+i\Delta z\)での速度 \(V(z_o+i\Delta z,t)\)。
The velocity \(V(z_o+i\Delta z,t)\) at a location \(z=z_o+i\Delta z\).
T[i] 位置\(z=z_o+(i+1/2)\Delta z\)での応力 \(\tau’(z_o+(i+1/2)\Delta z,t)\)。
The stress \(\tau’(z_o+(i+1/2)\Delta z,t)\) at a location \(z=z_o+(i+1/2)\Delta z\).
V2tau_coefficients[i] (\ref{eq.tau.evolution})式に登場する係数 \(\frac{\Delta t}{\Delta z}V_s(z_o+(i+1/2)\Delta z)^2\)。
The coefficient \(\frac{\Delta t}{\Delta z}V_s(z_o+(i+1/2)\Delta z)^2\) which appears in eq. (\ref{eq.tau.evolution}).
tau2V_coefficient (\ref{eq.V.evolution})式に登場する係数 \(\frac{\Delta t}{\Delta z}\)。
The coefficient \(\frac{\Delta t}{\Delta z}\) which appears in eq. (\ref{eq.V.evolution}).

また、\(z\)方向の格子点数を\(N_z\)として \(z=z_o\)は観測点位置、\(z=z_o+(N_z-1/2)\Delta z\)は地表を表すのであった。 このことと上の表からV[0]は観測点での速度、 T[Nz-1]は地表での応力を表すことが分かる。
Remember that \(z=z_o\) represents the station location, and \(z=z_o+(N_z-1/2)\Delta z\) does the ground surface, where \(N_z\) is the number of grid nodes in the \(z\)-direction. Then from the table below, it is understood that V[0] is the velocity at the station, and T[Nz-1] is the stress at the ground surface.
以上を踏まえて上記のコードを見ていこう。 まず紫色で示した部分(下記に再掲)に注目する。
We now proceed to look at the code above. First, pay attention to the purple part (reproduced below).
for(i=0;i<parameters.Nz-1;i++){
     T[i]+=V2tau_coefficients[i]∗(V[i+1]-V[i]);
}

ここでは位置\(z=z_o+(i+1)\Delta z\)での速度V[i+1]と \(z=z_o+i\Delta z\)での速度V[i]の差分を用いて 位置\(z=z_o+(i+1/2)\Delta z\)での応力の値T[i]を更新しており、 上記の(\ref{eq.tau.evolution})式に対応する。 時刻については引数に現れないが、速度場の時刻よりも\(\Delta t/2\)だけ後の時刻での 応力場を計算していることが(\ref{eq.tau.evolution})式との比較から分かる。 \(i=N_z-1\)について計算を行っていないのは地表のグリッドだからである。 時間発展ループを始める前に全ての場所での応力を0.0に初期化しているので、 値を更新しなければ応力は0となり、自由表面の境界条件が実現される。
Here, the stress T[i] at a location \(z=z_o+(i+1/2)\Delta z\) is updated based on the difference between the velocities V[i+1] at a location \(z=z_o+(i+1)\Delta z\) and V[i] at \(z=z_o+i\Delta z\), corresponding to eq. (\ref{eq.tau.evolution}). Although the time does not appear in the arguments, eq. (\ref{eq.tau.evolution}) shows that the stress at a time by \(\Delta t/2\) later than that of velocity is calculated. The calculation for \(i=N_z-1\) is skipped because it is the grid node on the ground surface. Before starting the temporal evolution loop, the stress at all the places are initialized to 0.0. Therefore the stress is kept zero if an update is skipped, thereby realizing the free surface boundary condition.
次に緑色で示した部分(下記に再掲)に注目しよう。
Next, let us focus on the green part (reproduced below).
for(i=1;i<parameters.Nz;i++){
     V[i]+=tau2V_coefficient∗(T[i]-T[i-1]);
}
V[0]=waveform_use.value[k];

ここでは位置\(z=z_o+(i+1/2)\Delta z\)での応力T[i]と \(z=z_o+(i-1/2)\Delta z\)での応力T[i-1]の差分を用いて 位置\(z=z_o+i\Delta t\)での速度の値V[i]を更新しており、 上記の(\ref{eq.V.evolution})式に対応する。 V[0]は観測点での速度を表すので 別扱いで実際に観測された波形の速度値waveform_use.value[k] を代入している。 これによって計算領域下端(\(z=z_o\))での境界条件も満たされる。
Here, the velocity V[i] at a location \(z=z_o+i\Delta z\) is updated based on the difference between the stresses T[i] at a location \(z=z_o+(i+1/2)\Delta z\) and T[i-1] at \(z=z_o+(i-1/2)\Delta z\), corresponding to eq. (\ref{eq.V.evolution}). For V[0], which represents the velocity at the station, the actual observed velocity value waveform_use.value[k] is inserted. This realizes the boundary condition at the lower end of the computational domain (\(z=z_o\)).
最後に青で示した部分(下記に掲載)に注目しよう。
Finally, let us focus on the blue part (reproduced below).
waveform_output.value[k]=V[parameters.Nz-1];

変数waveform_outputはstruct sequence型構造体で、 計算した波形を保存するための変数として用意している。 最終的に出力するのはこの波形である。 右辺のV[parameters.Nz-1]は 位置\(z=z_o+(N_z-1)\Delta z=z_s-\Delta z/2\)での速度を表す (\(z_s\)は地表の標高)。 本当は地表での速度波形が得られれば良いが、 staggered gridの配置の関係で得られない。 そこで地表に最も近い半グリッド下の格子点での速度値を 計算波形として保存している。
Here, the variable waveform_output is a struct sequence-type structure used to preserve the calculated waveform. V[parameters.Nz-1] in the right hand side represents the velocity at a location \(z=z_o+(N_z-1)\Delta z=z_s-\Delta z/2\), where \(z_s\) is the elevation of the ground surface. What is actually needed is the velocity waveform at the ground surface, which could, however, not be obtained due to the staggered grid configuration. Instead, the velocity value nearest to and by half grids below the ground surface is preserved as the calculated waveform.