% 表題  移流方程式の差分解法: MPDATA スキーム  1 次元 MPDATA
%
% 履歴  1997/11/21 小高正嗣
%       1997/11/25 小高正嗣
%       1998/01/20 小高正嗣
%       1998/05/11 小高正嗣
%
\section{1 次元 MPDATA}

     MPDATA は上流差分スキームで生じた拡散を適当な方法で補正するスキームで
     ある. 以下ではまず上流差分スキームとその誤差について解説し, その
     後 MPDATA の導出を行なう.
     
     \subsection{上流差分スキーム}

     解く方程式は以下の移流方程式である.

     \begin{equation}
       \DP{\psi }{t} + \Ddiv (\Dvect{v}\psi ) = 0. \Deqlab{eq: 1}
     \end{equation}

     ここで $\psi $ を正のスカラー量, $\Dvect{v}$ は速度である. 簡単のため
     定常で非発散な速度場を仮定する. 速度場が非定常または非発散の場合につい
     ては後述する. 物理量の格子配置はスタガードにとる.
     
     1 次元の場合\Deqref{eq: 1}式は,

     \begin{equation}
       \DP{\psi }{t} + \DP{}{x}(u\psi )= 0, \Deqlab{eq: 2}
     \end{equation}

     となる. 空間方向を上流差分を用いて差分化する. $x=i\Delta x,
     t=n\Delta t$ での $\psi $ を $\psi _{i}^{n}$ と表すことにすると,
     \Deqref{eq: 2}式は以下のように差分化される.

     \begin{equation}
%       \psi _{i}^{n+1} = \psi _{i}^{n} - 
       \DP{\psi }{t} = -
       \{F(\psi _{i}^{n},\psi _{i+1}^{n},u_{i+\frac{1}{2}}^{n}) - 
       F(\psi _{i-1}^{n},\psi _{i}^{n},u_{i-\frac{1}{2}}^{n})\}. \Deqlab{eq: 3}
     \end{equation}

     ただし,
     \begin{equation}
       F(\psi _{i}^{n},\psi _{i+1}^{n},u_{i+\frac{1}{2}}^{n})
       = [(u_{i+\frac{1}{2}}^{n}+|u_{i+\frac{1}{2}}^{n}|)\psi _{i}^{n} + 
          (u_{i+\frac{1}{2}}^{n}-|u_{i+\frac{1}{2}}^{n}|)\psi _{i+1}^{n}]
          \frac{1}{2\Delta x},
          \Deqlab{eq: 4}
     \end{equation}

     等である. 

     \subsubsection*{上流差分スキームの誤差}

     上流差分スキーム\Deqref{eq: 3}式の誤差を評価する. まず
     \Deqref{eq: 3}式を以下のように変形する.

     \begin{eqnarray}
       \DP{psi }{t}
       &=& 
         -\frac{1}{\Delta x}\left(u_{i+\frac{1}{2}}^{n}\frac{
         \psi _{i}^{n}+\psi _{i+1}^{n}}{2} - u_{i-\frac{1}{2}}^{n}\frac{
         \psi _{i}^{n}+\psi _{i-1}^{n}}{2}\right) \nonumber \\
       &&  -\frac{1}{2\Delta x}\left[|u_{i+\frac{1}{2}}^{n}|
         (\psi _{i}^{n} - \psi _{i+1}^{n}) - |u_{i-\frac{1}{2}}^{n}|
         (\psi _{i-1}^{n} - \psi _{i}^{n})\right] \Deqlab{eq: 5}
     \end{eqnarray}

     ここで $u_{i+\frac{1}{2}}^{n}, u_{i-\frac{1}{2}}^{n}, \psi _
     {i+1}^{n}, \psi _{i-1}^{n}$ を $u_{i}^{n}, \psi _{i}^{n}$ の周り
     のテーラー展開で表現する. $|u_{i+\frac{1}{2}}^{n}|$ 等は
     $|u|_{i+\frac{1}{2}}^{n}$ と考えて $|u|_{i}^{n}$ の周りの展開として
     表現する.

     \begin{eqnarray*}
       u_{i+\frac{1}{2}}^{n}&=&u_{i}^{n}+\left.\DP{u}{x}\right|_{i}^{n}
                  \frac{\Delta x}{2}
                  + \frac{1}{2!}\left.\DP[2]{u}{x}\right|_{i}^{n}
                  (\frac{\Delta x}{2})^{2}
                  + \cdot \cdot \cdot , \\
       u_{i-\frac{1}{2}}^{n}&=&u_{i}^{n}-\left.\DP{u}{x}\right|_{i}^{n}
                  \frac{\Delta x}{2}
                  + \frac{1}{2!}\left.\DP[2]{u}{x}\right|_{i}^{n}
                  (\frac{\Delta x}{2})^{2}
                  + \cdot \cdot \cdot , \\
       |u|_{i+\frac{1}{2}}^{n}&=&|u|_{i}^{n}+\left.\DP{|u|}{x}\right|_{i}^{n}
                  \frac{\Delta x}{2}
                  + \frac{1}{2!}\left.\DP[2]{|u|}{x}\right|_{i}^{n}
                  (\frac{\Delta x}{2})^{2}
                  + \cdot \cdot \cdot , \\
       |u|_{i-\frac{1}{2}}^{n}&=&|u|_{i}^{n}-\left.\DP{|u|}{x}\right|_{i}^{n}
                  \frac{\Delta x}{2}
                  + \frac{1}{2!}\left.\DP[2]{|u|}{x}\right|_{i}^{n}
                  (\frac{\Delta x}{2})^{2}
                  + \cdot \cdot \cdot , \\
       \psi _{i+1}^{n}&=&\psi _{i}^{n} + \left.\DP{\psi }{x}\right|_{i}^{n}
                  \Delta x
                  + \frac{1}{2!}\left.\DP[2]{\psi }{x}\right|_{i}^{n}
                  (\Delta x)^{2}
                  + \cdot \cdot \cdot , \\
       \psi _{i-1}^{n}&=&\psi _{i}^{n} - \left.\DP{\psi }{x}\right|_{i}^{n}
                  \Delta x
                  + \frac{1}{2!}\left.\DP[2]{\psi }{x}\right|_{i}^{n}
                  (\Delta x)^{2}
                  + \cdot \cdot \cdot . 
     \end{eqnarray*}

     それぞれの 2 次までをとって\Deqref{eq: 5}式の右辺へ代入する.
     右辺の$( \; )$ 内は

     \begin{eqnarray*}
       && -\frac{1}{\Delta x}\left[\left(u_{i}^{n}\psi _{i}^{n} 
         + \frac{u_{i}^{n}}{2}\left.\DP{\psi }{x}\right|_{i}^{n} \Delta x
         + \psi _{i}^{n}\left.\DP{u}{x}\right|_{i}^{n}\frac{\Delta x}{2}
          \right. \right. \\
       &&  + \left. \frac{u_{i}^{n}}{2}\left.\DP[2]{\psi }{x}\right|_{i}^{n}
          (\Delta x)^{2}
          + \left. \DP{u}{x}\right| _{i}^{n}
           \left. \DP{\psi }{x}\right| _{i}^{n}
         \frac{(\Delta x)^{2}}{2}
         + \frac{\psi _{n}^{i}}{2}\left.\DP[2]{u}{x}\right|_{i}^{n}
         \left(\frac{\Delta x}{2}\right)^{2}
         + O(\Delta x^{3}) \right) \\
       && -\left(u_{i}^{n}\psi _{i}^{n} 
         - \frac{1}{2}u_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n} \Delta x
         - \psi _{i}^{n}\left.\DP{u}{x}\right|_{i}^{n}
         \frac{\Delta x}{2} \right. \\
       && +\left. \left. 
         \frac{u_{i}^{n}}{2}\left.\DP[2]{\psi }{x}\right|_{i}^{n}
          (\Delta x)^{2}
       + \left.\DP{u}{x}\right|_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n}
         \frac{(\Delta x)^{2}}{2}
         + \frac{\psi _{n}^{i}}{2}\left.\DP[2]{u}{x}\right|_{i}^{n}
         \left(\frac{\Delta x}{2}\right)^{2}
         + O(\Delta x^{3}) \right) \right] \\
       &=& -\frac{1}{\Delta x}\left[
         u_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n}\Delta x
         + \psi _{i}^{n}\left.\DP{u}{x}\right|_{i}^{n} \Delta x 
       + O(\Delta x^{3}) \right] \\
       &=& -\left.\DP{}{x}(u\psi )\right|_{i}^{n} 
       + O(\Delta x^{2}),
     \end{eqnarray*}

     となる. 右辺の $[ \; ]$ 内は整理すると

     \begin{eqnarray*}
       && -\frac{1}{2\Delta x}\left[\left(
         -|u_{i}|\left.\DP{\psi }{x}\right|_{i}^{n}\Delta x
         -\frac{|u_{i}|}{2}\left.\DP[2]{\psi }{x}\right|_{i}^{n}(\Delta x)^{2}
         -\left.\DP{|u|}{x}\right|_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n}
         \frac{(\Delta x)^{2}}{2} +O(\Delta x^{3}) \right)\right. \\
       && -\left.\left(
         -|u_{i}|\left.\DP{\psi }{x}\right|_{i}^{n}\Delta x
         +\frac{|u_{i}|}{2}\left.\DP[2]{\psi }{x}\right|_{i}^{n}(\Delta x)^{2}
         +\left.\DP{|u|}{x}\right|_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n}
         \frac{(\Delta x)^{2}}{2} + O(\Delta x^{3}) \right)\right] \\
       &=& \frac{1}{2\Delta x}\left[
         |u_{i}|\left.\DP[2]{\psi }{x}\right|_{i}^{n}(\Delta x)^{2}+
         \left.\DP{|u|}{x}\right|_{i}^{n}\left.\DP{\psi }{x}\right|_{i}^{n}
         \Delta x^{2} + O(\Delta x^{3})\right] \\
       &=& \DP{}{x}\left.\left(\frac{1}{2}|u|\Delta x
         \DP{\psi }{x}\right)\right|_{i}^{n} + O(\Delta x^{2}).
     \end{eqnarray*}

     よって\Deqref{eq: 5}式は, 

     \begin{equation}
       \DP{\psi }{t} =
       -\left.\DP{}{x}(u\psi )\right|_{i}^{n}
         + \DP{}{x}\left.\left(\frac{1}{2}|u|\Delta x\DP{\psi }{x}\right)
         \right|_{i}^{n} +O(\Delta x^{2}) .\Deqlab{eq: 6}
     \end{equation}

     となる. すなわち 1 次の正確度を持つ上流差分スキームは
     移流拡散方程式

     \[
       \DP{\psi }{t} =
       -\DP{}{x}(u\psi )
         + \DP{}{x}\left(K\DP{\psi }{x}\right)
         , \; K =\frac{1}{2}|u|\Delta x ,
     \]

     を 2 次の正確度で近似したものに等しいことがわかる.

     \subsection{時間方向の差分と誤差}

     さらに\Deqref{eq: 6}の時間微分を前進差分で差分する.

     \begin{equation}
       \frac{\psi _{i}^{n+1}-\psi _{i}^{n}}{\Delta t}=
       -\left.\DP{}{x}(u\psi )\right|_{i}^{n}
         + \DP{}{x}\left.\left(\frac{1}{2}|u|\Delta x\DP{\psi }{x}\right)
         \right|_{i}^{n} +O(\Delta x^{2}) .
         \Deqlab{eq: 7}
     \end{equation}
     
     この場合の誤差を評価する. 左辺の $\psi _{i}^{n+1}$ を $\psi _{i}^{n}$ 
     の周りのテーラー展開で表現する.

     \[
       \psi _{i}^{n+1}=\psi _{i}^{n} + \left.\DP{\psi }{t}\right|_{i}^{n}
                  \Delta t
                  + \frac{1}{2!}\left.\DP[2]{\psi }{t}\right|_{i}^{n}
                  (\Delta t)^{2}
                  + \cdot \cdot \cdot . 
     \]
     
     これの 2 次までをとって\Deqref{eq: 6}に代入する.
     ただし右辺第 3 項の $t$ の 2 階微分は\Deqref{eq: 2}式を用いて空間微分で
     置き換える. すなわち,

     \[
       \DP[2]{\psi }{t} 
       = - \DP{}{x}\left[  u \left( -\DP{}{x}(u\psi ) \right) \right]
       = \DP{}{x}\left[  u \DP{}{x}(u\psi ) \right] 
       = \DP{}{x}\left[  u^{2} \DP{\psi }{x} \right]
     \]

     とする. ただし $u$ は定常かつ非発散としている. 
     これより\Deqref{eq: 7}式は以下のようになる.

     \begin{equation}
       \left.\DP{\psi }{t}\right|_{i}^{n}
       = -\left.\DP{}{x}(u\psi )\right|_{i}^{n}
         + \DP{}{x}\left.\left[\frac{1}{2}(|u|\Delta x-\Delta t u^{2})
             \DP{\psi }{x}\right]\right|_{i}^{n} +O(\Delta x^{2},\Delta t^{2}).
         \Deqlab{eq: 8}
     \end{equation}
     
     よって\Deqref{eq: 2}式を時間, 空間方向にそれぞれ 1 次の正確度を持つ
     上流・前進差分スキームで差分表現したものは, 移流拡散方程式

     \[
       \DP{\psi }{t} + \DP{}{x}(u\psi ) =
         \DP{}{x}\left(K_{impl}\DP{\psi }{x}\right),  
     \]

     $(K_{impl} = 0.5[|u|\Delta x-\Delta t u^{2}])$を 2 次の正確度で近似した
     ものであることがわかる.

     \newpage
     \subsection{1 次元 MPDATA}
     
     MPDATA では\Deqref{eq: 8}式に現れた数値拡散の効果を補正することを考える.
     まず数値拡散による $\psi $ の変化

     \[
       \DP{}{x}\left(K_{impl}\DP{\psi }{x}\right),  
     \]

     が, 
     \begin{equation}
       - \DP{}{x}(u_{d}\psi ), \Deqlab{eq: 9}
     \end{equation}

     と表されるような拡散速度(diffusion velocity) $u_{d}$

     \begin{equation}
       u_{d}=\left\{
         \begin{array}{ccl}
           -\frac{K_{impl}}{\psi}\DP{\psi }{x}, & \mbox{if} &\psi >0 \\
           0, & \mbox{if} & \psi =0. 
         \end{array}
         \right. \Deqlab{eq: 10}
     \end{equation}

     の移流によるものであるとする. $u_{d}$ の移流による変化分を補正するため
     反拡散速度(antidiffusive velocity) $\tilde{u}$ 
     \begin{equation}
       \tilde{u}= -u_{d},  \Deqlab{eq: 11}
     \end{equation}

     を定義し, これを用いた上流差分\Deqref{eq: 3}をもう一度計算する.
     全体の手続きをまとめると,

     \begin{screen}
     \begin{enumerate}
       \item 普通に上流差分を計算する.
     \begin{equation}
       \psi _{i}^{*} = \psi _{i}^{n} - 
       \{F(\psi _{i}^{n},\psi _{i+1}^{n},u_{i+\frac{1}{2}}^{n}) - 
       F(\psi _{i-1}^{n},\psi _{i}^{n},u_{i-\frac{1}{2}}^{n})\}. 
       \Deqlab{eq: 12}
     \end{equation}
     
       \item $\tilde{u}$ を求め, それを用いて上流差分を計算する.
     \begin{equation}
       \psi _{i}^{n+1} = \psi _{i}^{*} - 
       \{F(\psi _{i}^{*},\psi _{i+1}^{*},\tilde{u}_{i+\frac{1}{2}}) - 
       F(\psi _{i-1}^{*},\psi _{i}^{*},\tilde{u}_{i-\frac{1}{2}})\}. 
       \Deqlab{eq: 13}
     \end{equation}

     ただし,
     \begin{equation}
       \tilde{u}_{i+\frac{1}{2}} = \frac{
       (|u_{i+\frac{1}{2}}^{n}|\Delta x - \Delta t u_{i+\frac{1}{2}}^{2})
       (\psi _{i+1}^{*}-\psi _{i}^{*})}
       {(\psi _{i}^{*}+\psi _{i+1}^{*}+\epsilon )\Delta x}. \Deqlab{eq: 14}
     \end{equation}
     \end{enumerate}
     \end{screen}

     となる. $\epsilon $ は 0 割を防ぐための微小量で Smolarkiewicz(1983) では 
     $10^{-15}$ としている. \Deqref{eq: 13}, \Deqref{eq: 14}式を用いた
     計算を繰り返すことで精度はよくなるが, Smolarkiewicz(1984) では 3 回以上
     繰り返しても計算量が増えるわりには精度が上がらなかったとしている.

     \newpage
     \subsection{スキームの適合性と安定性}

     スキーム\Deqref{eq: 12}, \Deqref{eq: 13}の適合性と安定性について
     考察する. まず適合性を調べる. \Deqref{eq: 12}式は $\Delta t, \Delta x
     \rightarrow 0$ の極限で\Deqref{eq: 2}式に一致する.  \Deqref{eq: 13}
     式だが $\Delta t, \Delta x \rightarrow 0$ の極限で\Deqref{eq: 14}式から,

     \[
       \tilde{u} \rightarrow 0, 
     \]

     となることがわかるので, \Deqref{eq: 13}は

     \[
       \DP{\psi }{t} = 0,
     \]

     となる. よって\Deqref{eq: 12}, \Deqref{eq: 13}のスキームは適合性がある
     ことがわかる.

     次に安定性を調べてみる. スキーム\Deqref{eq: 12}の安定条件は

     \begin{equation}
       \mbox{max}\left(\frac{|u_{i+\frac{1}{2}}|\Delta t}
         {\Delta x}\right) \leq 1, \Deqlab{eq: 15}
     \end{equation}
     
     である. スキーム\Deqref{eq: 13}も形は同じ上流差分なので安定条件は,

     \begin{equation}
       \mbox{max}\left(\frac{|\tilde{u}_{i+\frac{1}{2}}|\Delta t}
         {\Delta x}\right) \leq 1,
       \Deqlab{eq: 16}
     \end{equation}

     となる. $\tilde{u}_{1+\frac{1}{2}}$ に\Deqref{eq: 14}式を代入すると,

     \begin{equation}
       \mbox{max}\left[\underbrace{\frac{|u_{i+\frac{1}{2}}|\Delta t}
         {\Delta x}}_{A}\times \underbrace{
         \left( 1 - \frac{|u_{i+\frac{1}{2}}|\Delta t}{\Delta x}\right)}_{B}
         \times \underbrace{
         \frac{|\psi _{i+1}^{*}-\psi _{i}^{*}|}
       {(\psi _{i}^{*}+\psi _{i+1}^{*}+\epsilon )}}_{C}\right]\leq 1,
       \Deqlab{eq: 17}
     \end{equation}

     となる. \Deqref{eq: 15}式が成り立つとすれば $A, B$ は 1 以下である.
     $C$ も 1 以下となるので\Deqref{eq: 15}式が成り立てば\Deqref{eq: 17}
     式も成り立つことがわかる. 特に\Deqref{eq: 15}が成り立つとき $A\times B$ 
     の最大値は 0.25 であるので $C$ が 1 に近いときでもスキーム\Deqref{eq: 13}
     はかなり安定である.