% 表題   DCPAM5  コード解説 -- 力学過程 -- セミインプリシット時間積分
%
% 履歴
%\Drireki{2009/10/05  高橋 芳幸}
%\Drireki{2009/07/11  森川 靖大}
%\Drireki{2009/02/25  森川 靖大}
%


\section{セミインプリシット時間積分}
\label{sec:セミインプリシット時間積分}

注: 本節は書き換え中であるため, 支配方程式系とその離散化の文章とは対応していない. 

\subsection{セミインプリシット時間積分の概要}
\label{subsec:セミインプリシット時間積分の概要}

セミインプリシット法については, 支配方程式系とその離散化ドキュメント
の 3.5 節に解説があるので
参照のこと. 時刻 $t-\Delta t$ から $t+\Delta t$ へのセミインプリシットの
計算手順は以下のようにまとめられる. 

\begin{enumerate}
 \item 発散項 $\Dvect{D}$ のみに関する単一のセミインプリシット方程式, 
   \begin{align}
    \Deqlab{code-smiimp:発散項のみに関する単一のセミインプリシット方程式}
      \tilde{\underline{M}}^{m}_{n}\,
      \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}
      = \tilde{\Dvect{f}}^{m}_{n}
   \end{align}
   を解く. ここで
   \[
     \dS \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} \equiv 
      \frac{1}{2}
        \left(
          \tilde{\Dvect{D}}^{m,t+\Delta t}_{n}
        + \tilde{\Dvect{D}}^{m,t-\Delta t}_{n} \right)
   \]
   である. 

 \item $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$
       と時刻 $t$ における時間変化の非重力波
       (NG)項から各変数の $t+\Delta t$ 
       での値を求める.

\end{enumerate}

以上の手順の説明は \ref{s:semiimp詳細} 節で行う.

\subsection{セミインプリシット時間積分の詳細}\label{s:semiimp詳細}

\ref{subsec:セミインプリシット時間積分の概要}節で述べたように,
セミインプリシット時間積分
は2つのステップに分けられる. これを詳しく書くと以下のようになる.

\noindent
● 第1段階: $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ の計算\label{semiimp手続き1}

\Deqref{code-smiimp:発散項のみに関する単一のセミインプリシット方程式}
を解くのが第1段階
である. ここでは支配方程式系とその離散化ドキュメント
の 3.5 節と同様, 太字で鉛直に離散化したベクトルを表
し, 下線で行列を表す. $\Dvect{D}$ は発散で,
%
\begin{align}
 \Deqlab{code-smiimp:Dbar^tの式}
  \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} \equiv 
      \frac{1}{2} (\tilde{\Dvect{D}}^{m,t+\Delta t}_{n}
    + \tilde{\Dvect{D}}^{m,t-\Delta t}_{n})
\end{align}
%
である. 行列 $\tilde{\underline{M}}^{m}_{n}$ は
%
%\\-----------------------------------ここから消す予定
%
\begin{align}
%\Deqlab{code-smiimp:行列Mの式}
 \begin{split}
  \tilde{\underline{M}}^{m}_{n} &\equiv
    ( 1            -2\Delta t             \tilde{\cal D}_{H,n}^{m}     )
    ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
  \\
 & \qquad
             - ( \Delta t )^{2}
                 \left\{   \underline{W} \ \underline{h} 
                        + ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} )
                          \Dvect{G} \Dvect{C}^{T} \right\}
                          \left( - \frac{n(n+1)}{a^2} \right)
 \end{split}
\end{align}
%
%\\-----------------------------------ここまで消す予定
%%
%\begin{align}
%\Deqlab{code-smiimp:行列Mの式}
% \begin{split}
%  \tilde{\underline{M}}^{m}_{n} &\equiv
%                ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
%  \\
% & \qquad
%             - ( \Delta t )^{2}
%                 \left( - \frac{n(n+1)}{a^2} \right)
%                 \left\{
%                    \underline{W}
%                    ( \underline{I} -2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )\
%^{-1}
%                    \underline{h}
%                  + \Dvect{G} \Dvect{C}^{T}
%                 \right\}
% \end{split}
%\end{align}
%%
%であり, ベクトル $\Dvect{f}$ は,
%%
%\\-----------------------------------ここから消す予定
%
\begin{align}
 \Deqlab{code-smiimp:ベクトルfの式}
 \begin{split}
  \tilde{\Dvect{f}}^{m}_{n}
 & =   ( 1            -2\Delta t             \tilde{\cal D}_{H,n}^{m}     )
       ( \underline{I}- \Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
           \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}
       + ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} ) \Delta t 
           \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG}  
  \\
 & \qquad \quad
       -  \Delta t \left( - \frac{n(n+1)}{a^2} \right) 
                   \Biggl\{ 
                            ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} ) \tilde{\Dvect{\Phi}}_{s,n}^{m}
                   \Biggr.
  \\
 & \qquad  \hspace{20mm}
                          + \underline{W} 
                            \left[ ( 1-\Delta t \tilde{\cal D}_{H,n}^{m} ) 
                                    \tilde{\Dvect{T}}^{m,t-\Delta t}_{n}
                                  + \Delta t 
                                      \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t}     
                                      \right)^{\rm NG} \right]
  \\
  & \qquad
                   \left.  \hspace*{20mm} 
                          + ( 1-2\Delta t \tilde{\cal D}_{H,n}^{m} ) \Dvect{G} 
                            \left[ \tilde{\pi}^{m,t-\Delta t}_{n}
                                  + \Delta t
                                     \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}  
                            \right]
                   \right\}
 \end{split}
\end{align}
%
%\\-----------------------------------ここまで消す予定
%%
%\begin{align}
% \Deqlab{code-smiimp:ベクトルfの式}
% \begin{split}
%  \tilde{\Dvect{f}}^{m}_{n}
% &   =  ( \underline{I}- \Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} )
%           \tilde{\Dvect{D}}^{m,t-\Delta t}_{n}
%       + \Delta t \left( \DP{\tilde{\Dvect{D}}^{m}_{n}}{t} \right)^{\rm NG}
%  \\
% & \qquad \quad
%       -  \Delta t \left( - \frac{n(n+1)}{a^2} \right)
%             \Biggl[
%                \tilde{\Dvect{\Phi}}_{s,n}^{m}
%             \Biggr.
%  \\
% & \qquad  \hspace{20mm}
%              + \underline{W}
%                ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )^{-1}
%                \left\{
%                       ( 1-\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m} )
%                       \tilde{\Dvect{T}}^{m,t-\Delta t}_{n}
%                     + \Delta t \left( \DP{\tilde{\Dvect{T}}^{m}_{n}}{t} \right)^{\rm NG}
%                \right\}
%  \\
%  & \qquad
%                \left.  \hspace*{20mm}
%                       + \Dvect{G}
%                         \left\{ \tilde{\pi}^{m,t-\Delta t}_{n}
%                               + \Delta t
%                                  \left( \DP{\tilde{\pi}^{m}_{n}}{t} \right)^{\rm NG}
%                         \right\}
%                \right]
% \end{split}
%\end{align}
%
である. 

$\tilde{\underline{M}}^{m}_{n}$ はサブルーチン {\tt SemiImplMatrix} で計算され, LU分解される. 
$\tilde{\underline{M}}^{m}_{n}$ は時間刻みの大きさが変らない限り設定し直す必要がない.
$\tilde{\Dvect{f}}^{m}_{n}$ の計算と
\Deqref{code-smiimp:発散項のみに関する単一のセミインプリシット方程式}
を解く作業
はサブルーチン {\tt TimeIntegration} で行われる. 

\noindent
● 第2段階: 時間積分 \\
時間積分も {\tt TimeIntegration} が行う.
第1段階で $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$ が
求まったので発散の時間積分は容易である:
%
\begin{align}
\Deqlab{code-smiimp:Dの時間積分}
 \tilde{\Dvect{D}}^{m,t+\Delta t}_{n} =
    2 \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t} - \tilde{\Dvect{D}}^{m,t-\Delta t}_{n} .
\end{align}
%
他の物理量についても, $\overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}$
と時刻 $t$ における時間変化率の非重力波(NG)項から,
$t+\Delta t$ での値が求まる:
%
\begin{align}
\Deqlab{code-smiimp:D以外の時間積分}
 \tilde{\cal X}^{m,t+\Delta t}_n =
   (\tilde{\gamma}_{{\cal X},n}^{m})^{-1}
   \left\{
        \tilde{\cal X}^{m,t-\Delta t}_n
      + 2\Delta t \left[ \left( \DP{\tilde{\cal X}^{m}_{n}}{t} \right)^{\rm NG}
                        + \tilde{\cal G}_{{\cal X},n}^{m}
                          \overline{\tilde{\Dvect{D}}^{m}_{n}}^{t}
                  \right]
   \right\} \, .
\end{align}
%
ここで,
%
\renewcommand{\arraystretch}{1.5}
\begin{equation}
 \begin{array}{cccccc}
  \dS \tilde{\cal X}^{m}_{n}
     & = & \dS \tilde{\pi}^{m}_{n},
         & \dS \tilde{\Dvect{\zeta}}^{m}_{n},
	 & \dS \tilde{\Dvect{T}}^{m}_{n},
	 & \dS \tilde{\Dvect{q}}^{m}_{n}, \\
  \dS \left(\DP{\tilde{\cal X}^{m}_{n}}{t}\right)^{\rm NG}
     & = &\dS \left(\DP{\tilde{\pi}^{m}_{n}}{t}\right)^{\rm NG},
         & \dS \left(\DP{\tilde{\Dvect{\zeta}}^{m}_{n}}{t}\right)^{\rm NG}, 
	 & \dS \left(\DP{\tilde{\Dvect{T}}^{m}_{n}}{t}\right)^{\rm NG}, 
	 & \dS \left(\DP{\tilde{\Dvect{q}}^{m}_{n}}{t}\right)^{\rm NG},\\
  \dS \tilde{\gamma}_{{\cal X},n}^{m}
     & = & \dS 1, & \dS ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{M} }_{n}^{m} ), 
         & \dS ( \underline{I}-2\Delta t \underline{ \tilde{\cal D}_{H} }_{n}^{m}),
	 & \dS (1-2\Delta t \tilde{\cal D}_{q,n}^{m}),  \\
  \dS \tilde{\cal G}_{{\cal X},n}^{m}
     & = & \dS -\Dvect{C}^T,
         & \dS \underline{0},
	 & \dS -\underline{h},
	 & \dS \underline{0}. 
 \end{array}
\end{equation}
%
\renewcommand{\arraystretch}{1}

以下の各節で, 具体的な計算手順とプログラムソースとの対応を述べる.