%--------------------------------------------------------------
\documentclass[ja=standard,xelatex,9pt]{bxjsarticle}
%--------------------------------------------------------------

%--------------------------------------------------------------
%packages%{
% #region
\usepackage{bigints}
\usepackage{changepage}
\usepackage[utf8]{inputenc}
\usepackage{lipsum}
\usepackage[bottom]{footmisc}
\usepackage{subcaption}
\usepackage{zxjatype}
\usepackage[hiragino]{zxjafont} 
\usepackage{tcolorbox}
\usepackage{amsmath}
\usepackage{nccmath}
\usepackage{amsthm}
\usepackage{empheq}
\usepackage{amssymb}
\usepackage{bm}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{tocloft}
\usepackage{ascmac}
\usepackage{latexsym}
\usepackage[version=4]{mhchem}
\usepackage{listings, xcolor}
\usepackage{mathtools}
\usepackage[colorlinks=true, linkcolor=blue , citecolor = blue]{hyperref}
\usepackage{array}
\usepackage{url}
\usepackage[english]{babel}
\usepackage[square,numbers]{natbib}
\usepackage{physics}
\usepackage{caption}
\usepackage{ulem}
\usepackage[section]{placeins}
\usepackage{xeCJKfntef}
\usepackage{longtable}
\usepackage{float}
\usepackage{multicol}
\usepackage{mathcomp}
\usepackage{IEEEtrantools}
\usepackage{flafter}
\usepackage{enumitem}
\usepackage[nospread]{cuted}
\setlength{\stripsep}{5pt plus 1pt minus 1pt}
% #endregion
%--------------------------------------------------------------%}

%--------------------------------------------------------------
%define%{
% #region

%%skip,space関係
\newcommand{\vs}[1]{\vskip#1\baselineskip}
\newcommand{\parlskip}[1]{\par\setlength{\leftskip}{#1pt}}

%equation環境
\newcommand{\eqn}[1]{\begin{equation*}{#1}\end{equation*}}
\newcommand{\neqn}[1]{\begin{equation}{#1}\end{equation}}
\newcommand{\leqn}[2]{\begin{fleqn}[#1pt]\begin{equation*}{#2}\end{equation*}\end{fleqn}[#1pt]}
\newcommand{\leqnarray}[2]{\begin{fleqn}[#1pt]\begin{eqnarray*}{#2}\end{eqnarray*}\end{fleqn}}
\newcommand{\eqns}[2]{\begin{equation*}\begin{array}{#1}#2\end{array}\end{equation*}}
\newcommand{\neqns}[2]{\begin{IEEEeqnarray}{#1}#2\end{IEEEeqnarray}}
\newcommand{\neqnarray}[1]{\begin{eqnarray}#1\end{eqnarray}}

%定理・定義・証明等の記述（数学）
\newcommand{\setting}[1]{\vs{0.3}\underline{\textbf{Setting}}\vs{0.3}\parlskip{20}{#1}\parlskip{0}}
\newcommand{\thm}[1]{\vs{0.3}\underline{\textbf{Thm}}\vs{0.3}\parlskip{20}{#1}\parlskip{0}}
\newcommand{\Def}[1]{\vs{0.3}\underline{\textbf{Def}}\vs{0.3}\parlskip{20}{#1}\parlskip{0}}
\newcommand{\rem}[1]{\vs{0.5}\large{\underline{{\textit{remark}}}}\normalsize\vs{0.5}{#1}}
\newcommand{\Proof}[1]{\vs{0.5}\large{\underline{\textit{Proof}}}\normalsize\vs{0.5}{#1}\setcounter{equation}{0}}
\newcommand{\free}[2]{\vs{0.5}\large{\underline{\textit{#1}}}\normalsize\vs{0.5}{#2}\setcounter{equation}{0}}

%tinysections
\newcommand{\tinysection}[3]{\vs{0.5}\large\underline{\textbf{(#1)：#2}}\normalsize\vs{0.5}{#3}}
\newcommand{\tinytinysection}[4]{\vs{0.3}\normalsize{\underline{\textbf{(#1.#2)：#3}}}\vs{0.3}{#4}}
\newcommand{\tinytinytinysection}[5]{\vs{0.2}\normalsize{\underline{(#1.#2.#3)：#4}}\vs{0.2}{#5}}

%微分
\newcommand{\dif}[2]{\dfrac{d#1}{d#2}}
\newcommand{\difop}[1]{\dfrac{d}{d#1}}
\newcommand{\lagdif}[2]{\frac{D#1}{D#2}}

%偏微分
\newcommand{\pdif}[2]{\dfrac{\partial{#1}}{\partial#2}}
\newcommand{\pdifop}[1]{\dfrac{\partial}{\partial#1}}
\newcommand{\pdifl}[1]{\partial_{#1}}
\newcommand{\pdiftwice}[2]{\dfrac{\partial^{2}{#1}}{\partial {#2} ^{2}}}

%積分
\newcommand{\dint}{\displaystyle\int}
\newcommand{\diint}{\displaystyle\iint}
\newcommand{\diiint}{\displaystyle\iiint}

%総和
\newcommand{\sumn}[2]{\sum_{n=#1}^{#2}}
\newcommand{\sumi}[2]{\sum_{i=#1}^{#2}}
\newcommand{\sumj}[2]{\sum_{j=#1}^{#2}}
\newcommand{\sumk}[2]{\sum_{k=#1}^{#2}}
\newcommand{\suml}[2]{\sum_{l=#1}^{#2}}
\newcommand{\summ}[2]{\sum_{m=#1}^{#2}}
\newcommand{\sump}[2]{\sum_{p=#1}^{#2}}
\newcommand{\sumq}[2]{\sum_{q=#1}^{#2}}
\newcommand{\sumfr}[3]{\sum_{#1=#2}^{#3}}
\newcommand{\dsum}{\displaystyle\sum}

%ベクトル
\newcommand{\vect}[1]{\mathbf{#1}}
\newcommand{\vecp}[2]{\mathbf{#1}^{#2}}
\newcommand{\vecn}[1]{\mathbf{#1}^n}
\newcommand{\vecm}[1]{\mathbf{#1}^m}
\newcommand{\pvect}[1]{\begin{pmatrix}#1\end{pmatrix}}

%括弧
\newcommand{\bracketnormal}[1]{\left(#1\right)}
\newcommand{\curlybracket}[1]{\left\{#1\right\}}

%ε-δ論法
\newcommand{\epsilondelta}[4]{\forall\;\epsilon>0,\exists\delta\;\; s.t, \;\; 0<\left| #1-#2 \right|<\delta \Rightarrow \left| #3(#1)-#4\right|<\epsilon}

%その他
\newcommand{\realpower}[1]{\mathbb{R}^{#1}}
\newcommand{\Real}{\mathbb{R}}
\newcommand{\Bf}[1]{\mathbf{#1}}
\newcommand{\limit}[2]{\lim_{#1 \to #2}}
\newcommand{\absdefined}[1]{\left| #1 \right|}
\newcommand{\map}[5]{#1:#2\to#3;\;\;#4\mapsto#5}
\newcommand{\inner}[1]{\left\langle #1 \right\rangle}
\newcommand{\vertline}[2]{\vs{0}\begin{center}\rule{#1\textwidth}{#2pt}\end{center}\vs{0}}
\newcommand{\coord}[2]{\left(#1\;,\;#2\right)}
\newcommand{\scipow}[1]{\times 10^{#1}}
\newcommand{\labeqn}[1]{\label{eqn:#1}}
\newcommand{\labfig}[1]{\label{fig:#1}}
\newcommand{\refeqn}[1]{（\ref{eqn:#1}）}
\newcommand{\reffig}[1]{Figure\ref{fig:#1}}

\usepackage{xcolor}
\usepackage{xparse}
\usepackage{titlesec}

\titlespacing*{\section}{0pt}{1.0\baselineskip}{0.5\baselineskip}
\titlespacing*{\subsection}{0pt}{0.8\baselineskip}{0.4\baselineskip}
\titlespacing*{\subsubsection}{0pt}{0.6\baselineskip}{0.3\baselineskip}



% #endregion
%--------------------------------------------------------------%}

%--------------------------------------------------------------
%title,author,etc%{
% #region
\title{コンピュータグラフィクスアニメーション課題}
\author{05242628 三田村彰大}
\date{}
% #endregion
%--------------------------------------------------------------%}

\begin{document}

%--------------------------------------------------------------
% #region
%code_area_settings%{
\lstset{
    basicstyle = {\ttfamily\small}, % 基本的なフォントスタイル
    frame = {tbrl}, % 枠線の枠線。t: top, b: bottom, r: right, l: left
    breaklines = true, % 長い行の改行
    numbers = left, % 行番号の表示。left, right, none
    showspaces = false, % スペースの表示
    showstringspaces = false, % 文字列中のスペースの表示
    showtabs = false, %　タブの表示
    keywordstyle = \color{blue}, % キーワードのスタイル。intやwhileなど
    commentstyle = {\color[HTML]{1AB91A}}, % コメントのスタイル
    identifierstyle = \color{black}, % 識別子のスタイル　関数名や変数名
    stringstyle = \color{brown}, % 文字列のスタイル
    captionpos = t % キャプションの位置 t: 上、b: 下
}
\captionsetup*[lstlisting]{labelformat=empty}
% #endregion
%--------------------------------------------------------------%}

%--------------------------------------------------------------
%document_settings%{
% #region

\captionsetup*[figure]{labelfont=bf,font=small}
\captionsetup*[table]{labelfont=bf,font=footnotesize}
\captionsetup*[subfigure]{position=bottom,font=footnotesize}
\renewcommand{\arraystretch}{1.15}

\setlength{\dbltextfloatsep}{10pt plus 2pt minus 2pt}
\setlength{\dblfloatsep}{8pt plus 2pt minus 2pt}

\setcounter{tocdepth}{3}
\setlength{\parindent}{0pt}
\newcommand{\setdisplaymathskip}{%
\setlength{\abovedisplayskip}{10pt plus 1pt minus 2pt}
\setlength{\belowdisplayskip}{10pt plus 1pt minus 2pt}
\setlength{\abovedisplayshortskip}{6pt plus 1pt minus 1pt}
\setlength{\belowdisplayshortskip}{6pt plus 1pt minus 1pt}
\setlength{\jot}{2pt}
}
\setlength{\jot}{3pt}
\newgeometry{top=18truemm,bottom=20truemm,left=14truemm,right=14truemm}
\renewcommand{\arraystretch}{1.5}
\bibliographystyle{abbrvnat}
\setlength{\columnsep}{5mm}
\setlength{\columnseprule}{0pt}
%asm_settings
\theoremstyle{definition}
\newtheorem{prop}{proposition}[section]
%subfigure_settings
\captionsetup*[subfigure]{position=bottom}
\captionsetup*[figure]{labelfont=bf}
\captionsetup*[table]{labelfont=bf}
%footnote_settings
\setlength{\footnotesep}{0.4\baselineskip} % 脚注どうしの間隔
\makeatletter
\renewcommand\@makefntext[1]{%
  % ---- 脚注本文の段落設定（同一脚注内） ----
  \setlength{\parindent}{0pt}%
  \setlength{\parskip}{0.4\baselineskip}%
  % ---- ぶら下げ＋マーク表示 ----
  \noindent
  \hb@xt@1.8em{\hss\@makefnmark}%
  \hspace{0.4em}%
  \hangindent=2.2em%
  \hangafter=1%
  #1%
}
\renewcommand{\footnoterule}{%
  \kern -3pt
  \hrule width 1.0\linewidth
  \kern 5pt
}
\makeatother
\makeatletter
\renewcommand{\bibsection}{} %"References" 見出しを無効化
\makeatother
\makeatletter
\renewenvironment{proof}[1][\proofname]{%
  \par\pushQED{\qed}\normalfont
  \trivlist
  \item[\hskip\labelsep\bfseries #1\@addpunct{.}]%
  \ignorespaces
}{%
  \popQED\endtrivlist\@endpefalse
}
\makeatother
% #endregion
\renewcommand{\abstractname}{}


%--------------------------------------------------------------%}


%--------------------------------------------------------------
%maketitle_makehedding_tableofcontents%{
% #region


\twocolumn[{
\maketitle
\begin{adjustwidth}{1cm}{1cm}
\begin{abstract}
　\textbf{概要}　粒子法、特にMPS法を実装し、実際に流体シミュレーションを行った。また異なる手法の比較という観点から、半陰解法（semi-implicit）および陽解法（explicit）を用いた場合についてそれぞれ計算を行い挙動の違いを観察した。その上で、粒子法が波面のダイナミックな動きを再現することができることや、境界条件の設定を非常に自由に行えることなどを体感した。また、陽解法が数値不安定であることや、安定性を達成できても非圧縮性をうまく再現できないことなどを実際に計算を可視化することで理解した。
\end{abstract}
\vs{0.7}
\noindent
\begin{minipage}{0.9\linewidth}
\raggedright
\begin{tabular}{@{}l@{\quad}l@{}}
※コード置き場： & \url{https://github.com/ramutami/mps_particle_fluid}\\[-0.3em]
※可視化動画置き場： & \url{https://shota-mitamura.com/works/ryuushi/ryuushhou/}\\[-0.3em]
※ AIとの会話ログ：&\url{https://chatgpt.com/g/g-p-6a31a79525788191ac9f0113ba5c1b07-yang-jie-fa-li-zi-fa/project}
\end{tabular}
\end{minipage}
\end{adjustwidth}
\vs{1.5}
}]
\setdisplaymathskip

% #endregion
%--------------------------------------------------------------%}

%start_document

\section{はじめに}
　流体のシミュレーションの主な手法としては、空間を格子に分割し各格子に物理量を割り当てることで流れを計算するEuler的な手法と、流体を粒子の集まりとして表現し各粒子の位置を計算点と考えるLagrange的な手法がある。前者は一般に格子法などと呼ばれ、後者は一般に粒子法などと呼ばれる。\cite{intro}
\vs{0}
　理論が整備されており微分の離散化が比較的単純な格子法に比べると、粒子法は理論がまだ成熟しておらずまた精度や数値安定性に問題を抱える手法である。しかし、格子法が苦手としている波面の表現や、連続体の変形・破断などは表現においては粒子法は非常に優れている。（\reffig{1}参照。）そのためゲームやCGなど、物理的な正確性よりも見栄えの美しさなどが優先される場面においては粒子法の利用が非常に盛んとなっている。
\vs{0.5}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/firsttry_velocity0063.png}
  \caption{粒子法が得意とする波面の表現の一例。今回の数値計算の結果を示している。}
  \label{fig:1}
\end{figure}

\FloatBarrier
　さて、授業では粒子法の手法の一つとしてSPH法が紹介されたが、今回はMPS法というこれとは異なる手法を用いた粒子法を実装することを考える。両者の違いとしては、SPH法が圧縮性流体の計算手法として開発された手法であるのに対しMPS法は非圧縮流体の計算手法として開発されたという経緯がある。そのためSPH法が安定性に優れる反面物理的に正確な挙動を表現することを苦手としているのに対し、MPS法は非圧縮性を保ちながら計算するためSPH法よりかは厳密であり、工学的な応用が多い（参考：\cite{intro2}）\footnote{SPHは宇宙物理やゲーム、MPSは製品開発におけるシミュレーションなどの工学分野、といった雰囲気の棲み分けがあるらしい…？}。
\vs{0}
　その上で、今回は\cite{koshizuka2014}を参考にFortranを用いてMPS法を用いた水柱崩壊を実装し、粒子法の性質について観察することを試みた。また、計算の途中で必要となる圧力をimplicitに解いた場合とexplicitに解いた場合でどのような性質の違いがあるのかについても観察した。
\vs{1}
　コードは\href{https://github.com/ramutami/mps_particle_fluid}{github}にアップロードしており、本レポートには書ききれなかった詳細な理論背景をREADMEにまとめている。また、\href{https://shota-mitamura.com/works/ryuushi/ryuushhou/}{可視化ページ}に計算結果の動画を置いている。

\section*{※参考図書の利用・スパコンの利用・AIの利用について}
　今回、MPS法の実装にあたって丸善出版より発売されている『粒子法入門：流体シミュレーションの基礎から並列計算と可視化まで』\cite{koshizuka2014}および付随しているサンプルコードを大いに参考にした。ただしコードは基本的に一から自分の手で書いている。\footnote{
  もともとC言語で書かれたサンプルコードをFortranに書き写す遊びを半年ほど前に行っており、今回はそれを再開する形で実装を行った。
}
\vs{0}
　また、今回計算に当たっては東大のWisteria Odysseyスーパーコンピュータシステムを利用した。たまたま工学部の計算機系の授業をとっていた関係で学生一人一人に計算資源（8ノード・48コア・一計算15分）が割り当てられていたため、そちらを利用した。\footnote{
  利用にあたっては授業用計算資源の趣旨に反しない範囲であることを意識した。件の工学部の授業でも粒子法は一つの大きなトピックとして扱われており実際にその並列化方法をコードを動かして学ぶことが求められていたため、実質的に授業の内容が被っていると判断し計算資源を利用した。
}
　また今回粒子法の実装にあたってAIをデバッグや並列化手法の提案、サンプルコードの解読などに利用した。またレポートの誤字脱字のチェックや内容の査読にもAIを用いている。詳しくは本文冒頭の会話ログのリンクを参照。

\section{理論}
　流体のシミュレーションにおいて基本となるのは次に示すナビエ・ストークス方程式である。
\neqn{
\dfrac{D\boldsymbol{u}}{Dt}  = -\dfrac{1}{\rho}\nabla P +\nu\nabla^2\mathbf{u}+\boldsymbol{g}\labeqn{1}
}
ただし以下では流体の流れを$\boldsymbol{u}$で表すことにする。また、上式において$P$は圧力、$\rho$は密度、$\nu$は流体の粘性を決める動粘性係数、$\boldsymbol{g}$は重力加速度である。以下では、粒子法、特にMPS法において\refeqn{1}がどのように数値的に解かれるかを見ていく。以下の内容は基本的に\cite{koshizuka2014}を参考にしている。

\subsection{粒子法の概要}
　最初に述べたように、粒子法においては流体を粒子の集まりとして表現し各粒子$i$の位置において物理量を評価する。よって、$i$番目の粒子位置$\boldsymbol{r}_i$においてナビエストークス方程式の右辺を評価し、それを用いて次のように流速$\boldsymbol{u}$を更新する。

\neqn{
\boldsymbol{u}^{k+1} _i = \boldsymbol{u}^k _i + \Delta t\cdot\left(-\dfrac{1}{\rho}\nabla P +\nu\nabla^2\mathbf{u}+\boldsymbol{g} \right)^{k} _i \labeqn{2}
}

ただし$k,k+1$は時間ステップを表す。また、流速が求まれば粒子の位置は$\boldsymbol{r}^{k+1} _i= \boldsymbol{r}^k _i + \Delta t\cdot  \boldsymbol{u}^{k}$として更新することができる。
　よって、\refeqn{1}式の離散化をどのように行うかということが本質であり、MPSやSPHといった手法の違いを理解する上でも重要である。

\subsection{MPS法における微分演算子の計算スキーム}
　ナビエストークス方程式の右辺を離散化するということはつまり微分演算子$\nabla,\nabla^2$を離散化するということである。MPS法においては、粒子からの距離に応じて減衰する重み関数を用いてある粒子の近傍に存在する粒子からの微分への寄与を足し込むことでこれらの値を計算する。具体的には、着目している粒子からの距離が近ければ近いほど大きくなり、同時に一定距離$R_e$より外側では重みが$0$となるような重みを考える。（\reffig{2}参照。）
\vs{0}


\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/weight.png}
  \caption{重み関数。赤点が着目粒子であり、影響半径$R_e$より内側の粒子の情報がその近さに応じて足し込まれる。}
  \label{fig:2}
\end{figure}

\vs{0}
　以下にMPS法で用いられる重み関数の一例を示す。この重み関数はあくまで一例であり、今回の場合は半陰解法を用いた手法で用いられているが例えば陽解法においては別の重み関数を用いることになる。（詳しくはコード置き場の\href{https://github.com/ramutami/mps_particle_fluid/blob/main/README.md}{README}を参照。）

\neqn{
  w \left( r\right) = \left\{ 
  \begin{aligned}
  &\dfrac{r _e}{r}-1&\quad\left(r<r _e\right)\\[2pt]
  &\; \; 0 &\quad\left(r>r _e\right)
\end{aligned} \right.
}

\vs{0}
　その上で、この重み関数を用いて微分演算子は次のように離散化することができる。

\neqns{l}{
\nabla\phi = \dfrac{d}{n^0}\displaystyle\sum _{i\neq j}\dfrac{\phi(\boldsymbol{X} _{j})-\phi(\boldsymbol{X} _{i})}{|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|}\dfrac{\boldsymbol{X} _{j}-\boldsymbol{X} _{i}}{|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|}w(|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|)\labeqn{4}\\[5pt]
\nabla^2 \phi= \dfrac{2d}{\lambda _0 n^0}\displaystyle\sum _{i\neq j}(\phi(\boldsymbol{X} _{j})-\phi(\boldsymbol{X} _{i}))w(|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|)\labeqn{5}
}

ただし$\phi$は着目している粒子$i$での物理量を表し、$\boldsymbol{X}_i,\boldsymbol{X}_j$は$i$番目、$j$番目の粒子の位置を表す。また$d$は空間の次元（例えば３次元の水柱の崩壊を計算するなら$d=3$）、$n^0$は初期状態における流体の数密度、$\lambda^0$は$n^0$を用いて計算されるパラメタ\footnote{$\lambda^0$は具体的には次のようん計算される。
  \neqn{
    \lambda _0 = \dfrac{1}{n^0}\displaystyle\sum _{i\neq j}|\boldsymbol{X} _j -\boldsymbol{X} _i|^2 w(|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|)
  }
}である。（詳しい導出は\cite{koshizuka2014}を参照。）
\vs{0}
　ここでは詳しい導出は行わないが、本質的に大事なのはそれぞれの微分演算子が実質的に「影響半径内に存在する粒子からの寄与の重み付きの足し込み」によって計算されているということである。\footnote{
  $w(|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|)$が影響半径の外側で$0$となるので$\dsum_{j\neq i}$は実質的に影響半径内の粒子についての足しあわせと考えることができる。
}これが後ほど並列化による計算高速化を考える上で重要となる。
\vs{0}
　また、\refeqn{4}\refeqn{5}式には初期状態での流体の数密度$n^0$が現れている。今一般の時刻における粒子$i$近傍の数密度$n_i$は

\neqn{
n_i = \displaystyle\sum _{i\neq j}{w(|\boldsymbol{X} _{j}-\boldsymbol{X} _{i}|)}
}

によって計算されるが、\refeqn{4}\refeqn{5}の導出にあたっては$n_i\sim n^0$という近似を用いている。これは

\begin{center}
  「流体の密度は初期基準値からずれない」
\end{center}

という非圧縮性をMPS法が織り込んでいることを意味する。実際、流体の密度変化を解消するように圧力で粒子を押し戻すという計算によってMPS法は非圧縮性を実現している。
\vs{1}
　さて、微分演算子の離散化スキームがもとまったので、我々が\refeqn{2}を計算するために知りたかった、「粒子$i$の位置においてのナビエストークス方程式の右辺の値」を計算することができる。ただし、MPS法においては先にまずナビエストークス方程式の右辺第二項（粘性項）と右辺第三項（重力項）のみを計算する。そしてそれらを用いて流速を更新したのち、改めて右辺第一項（圧力項）を計算して流速を再び更新する。式に示すと次のよう。

\neqns{l}{
\dfrac{\boldsymbol{u}^{k+\frac{1}{2}} _i-\boldsymbol{u}^k _i}{\Delta t} = \nu \left(\nabla^2 \boldsymbol{u}\right)^{k} _i + \boldsymbol{g}\labeqn{8}\\[5pt]
\dfrac{\boldsymbol{u}^{k+1} _i-\boldsymbol{u}^{k+\frac{1}{2}} _i}{\Delta t} = -\dfrac{1}{\rho^0}\left(\nabla P\right)^{k+\frac{1}{2}} _i
\labeqn{9}
}

　この時、粘性項（と重力項）は粒子の時刻$k$での情報から計算することができるが、圧力項を計算するには時刻$k$における圧力を知っていなければならない。しかし流体方程式は圧力の時間発展を記述する式を持たないので、前の時刻の圧力をステップ更新して現在の時刻における圧力を求めるといった手法を取ることができない。
\vs{0}
　すなわち、現在時刻における流速や密度の情報を使って都度その時刻における圧力を計算する必要がある。今回はその手法として次に述べる半陰解法と陽解法を利用する。

\subsection{圧力項の計算：半陰解法}
　半陰解法とは、圧力をポアソン方程式の解として求める手法である。今\refeqn{8}\refeqn{9}式と、非圧縮条件

\neqn{
  \dfrac{D\rho_i}{Dt} + \rho_i \nabla \cdot \boldsymbol{u}_i = 0
}

などを用いることで次を得ることができる。\footnote{
  詳しい導出は誌面の都合上割愛するが、\href{https://github.com/ramutami/mps_particle_fluid/blob/main/README.md}{README}にまとめてある。
}

\neqn{
  -\dfrac{1}{\rho _0}\left(\nabla^2 P \right)^{k+\frac{1}{2}} _i=\dfrac{1}{\Delta t^2}\dfrac{n _i^{k+\frac{1}{2}} -n^0}{n^0}\labeqn{12}
}

よってこのポアソン方程式を解けば圧力$P$を求めることができる。この時境界条件として、液面に存在する粒子の圧力を$0$とする。（\reffig{3}参照。）

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/surface.png}
  \caption{ポアソン方程式の境界条件に用いられる自由表面粒子の図。粒子数密度$n_i$が一定値を下回った粒子を液面に存在する自由表面粒子と判定する。（図中赤色で示された粒子が自由表面粒子。）}
  \label{fig:3}
\end{figure}

\subsection{圧力項の計算：陽解法}
　半陰解法はポアソン方程式を解くことを必要としており、計算コストが重くなりがちである。そこでよりexplicitな方法として圧力$P$を方程式の解ではなく直接計算する手法が知られており、陽解法と呼ばれる。
\vs{0}
　今、音速の定義 $c ^2=\partial P/ \partial \rho$ に従って圧力を $P _i\sim P^0+c^2 \left(\rho _i - \rho^0\right)$ テイラー展開し、 $\rho _i \sim \rho _0\frac{n _i}{n^0}$ 、 $P^0=0$ と考えることによって

\neqn{
P^{k+\frac{1}{2}} _i = c^2\dfrac{\rho^0}{n^0}\left(n _i^{k+\frac{1}{2}}-n^0\right) \labeqn{122}
}

という式を得ることができる。よってこの式を用いれば、方程式を解くことなく粒子数密度の情報から直接粒子$i$の位置における圧力を計算することが可能である。ただしこの陽解法を用いることによる計算コストの削減は数値安定性や非圧縮性の再現度とトレードオフの関係にあり、後々実装結果を議論する際に注目される。

\section{実装}
　ここまででMPS法の基本となる理論の解説を行ってきた。以下では実際に実装する上でプログラムがどのように構築されるかについて簡単に解説する。

\subsection{プログラムの概要}
　今回用いるプログラムの概要を\reffig{4}に示す。

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.7\columnwidth]{./image_tex/program.png}
  \caption{プログラムの概要}
  \label{fig:4}
\end{figure}

\vs{0}
　まず初めに粒子の初期配置を行う。この時、壁も流体と同じように粒子で表現する。この時この壁粒子に流体粒子が接近すると粒子数密度が上昇し非圧縮条件を満たすように圧力による押し戻しが発生することにより、「流体粒子を跳ね返す」という壁として満たしてほしい性質が満たされることになる。（\reffig{5}参照。）

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.7\columnwidth]{./image_tex/init.png}
  \caption{水柱崩壊における粒子の初期配置。黒が壁粒子、青が流体粒子。}
  \label{fig:5}
\end{figure}

\vs{0}
　粒子の初期配置が決定したら実際の計算に進む。まず初め先に示した$\nabla^2$の離散化スキームに従って粘性項を計算し、重力項による加速と足し合わせて得られる加速度で粒子の一回目の位置更新を行う。その後、半陰解法・陽解法のいずれかの手法により圧力を計算し、先に示した$\nabla$の離散化スキームに従って圧力項を計算する。その後、この圧力項による加速を元に粒子位置を今一度更新すれば、一ステップの計算が終了する。詳しいプログラムの内容については\href{https://github.com/ramutami/mps_particle_fluid/blob/main/README.md}{README}を参照。

\subsection{半陰解法における共役勾配法を用いた行列計算}
　さて、半陰解法においては\refeqn{12}に示されたポアソン方程式を解くわけであるが、先ほど示した$\nabla^2$に関する離散化スキームを用いればこれは次のように離散化できる。

\neqn{
  -\dfrac{1}{\rho _0}\dfrac{2d}{\lambda^0 n^0}\displaystyle\sum _{j\neq i}\left(P _j-P _i\right)w(|\boldsymbol{r} _j -\boldsymbol{r} _i |)=\dfrac{1}{\Delta t^2}\dfrac{n _i -n^0}{n^0}
}

よってポアソン方程式は$A\boldsymbol{P} = \boldsymbol{b}$という形の行列方程式に帰着されることがわかる。特に今、$A$が正定値実対称行列となっていることを示すことができるので\footnote{
  \href{https://github.com/ramutami/mps_particle_fluid/blob/main/README.md}{README}参照。
}、共役勾配法によって$P$を求めることができる。
\vs{0}
　この時、「近傍粒子を辿っていっても自由表面粒子に辿りつかないような内部粒子」が存在してしまうと$A$の正定値性は弱くなってしまうことがわかり、安定に計算するためには行列の正定値性を強くするような工夫が必要となる。今回の場合においては

\neqn{
  A\boldsymbol{P}=\boldsymbol{b}\; \leadsto\;\left(A+\dfrac{\kappa}{\Delta t^2}I\right)\boldsymbol{P} = \gamma \boldsymbol{b}
}

のようにして$A$の対角成分を大きくする工夫を行っている。

\subsection{近傍粒子探索・陽解法・並列化による計算量の削減}
　さて、半陰解法において共役勾配法を用いることができるとわかったので、一番計算が重くなる可能性のあったポアソン方程式の求解は粒子数$N$に対し$\mathcal{O}(kN^2)$となる。よって、それ以外の粘性項の計算や数密度の計算などに現れる$\sum_{i\neq j}$が$i$及び$j$の二重ループとして現れ$\mathcal{O}(N^2)$となることを考えると、ポアソン方程式の計算が計算量を決定することがわかる。
\vs{0}
　ただしこの時、先に述べたように$\sum_{i\neq j}$はほとんど全ての場合において重み関数をかけながら足し込むので、影響半径の外側で重みが$0$となることを踏まえると実質的な計算は着目粒子$i$についてのループといくつかの近傍粒子$j$についての足し合わせになり$\mathcal{O}(N)$となる。よって適切に近傍粒子のリストを保持することができれば、ポアソン方程式の求解の計算量は$\mathcal{O}(kN)$となり、これが支配的な計算オーダーになる。よって近傍粒子探索によって全体の計算量を抑えることができる。\footnote{
  ただし$\mathcal{O}(N^2)$の計算を行わないようにするには、「初めに各粒子の近傍粒子をリスト化するルーチン」を$i,j$について走査して$\mathcal{O}(N^2)$で計算してはいけない。実際には、空間をパケット化することによってこの近傍粒子探索を$\mathcal{O}(N)$で実現することができることが知られている。
}
\vs{0}
　いずれにせよ、半陰解法において計算のボトルネックとなるのはポアソン方程式の共役勾配法による求解の$\mathcal{O}(kN)$であり、計算量の削減にはこの圧力計算を軽量化する必要がある。この時、\refeqn{122}式に示したexplicitな方法によって圧力を求めれば、この圧力計算を$\mathcal{O}(N)$で実現することが可能になり、陽解法を用いることにより計算量が削減できる。
\vs{1}
　※ただし今回の計算においては自分のキャパの限界もあり近傍粒子探索を実装することができなかった。すなわち圧力の計算以外の諸計算が$\mathcal{O}(N^2)$で実現されており、陽解法を用いた計算の軽量化を確認することはあまり望めなかった。
\vs{1}
　また、今回は計算を東大のWisteria Odysseyスーパーコンピュータシステムを利用して行っており、8ノード48コアを用いることによって大幅な計算高速化を達成している。具体的にはMPIとopenMPを用いてFortranコードの並列化を行った。

\subsection{陽解法の不安定性と圧縮性}
　\cite{koshizuka2014}によれば、MPS法における数値安定性条件は次のようになる。

\neqns{ccl}{
\dfrac{u _\text{max}\Delta t}{l^0} &<& 0.2\labeqn{141}\\[4pt]
\dfrac{c\Delta t}{l^0} &<& 1.0\labeqn{142}
}

上が半陰解法と陽解法両方が満たさなければいけない数値安定性条件で、下は陽解法のみが満たさなければいけない数値安定性条件である。ただし上式において$l^0$は初期粒子間距離である。
\vs{0}
　この時、物理的な音速$c$の大きさと粒子の小ささを考えると、下の条件を満たすために陽解法においては$\Delta t$は相当小さな値を取らなければならない。言い換えれば、陽解法においては計算コストが低くなる代わりに数値安定性が悪くなってしまう。
\vs{0}
　このようにして生じる数値不安定性を回避する手段の一つとして、音速$c$を実際の値よりも小さくとるというものがある。すなわち$c=u_\text{max}/0.2$などと置くことで、\refeqn{141}が満たされるときに\refeqn{142}が自動的に満たされるようにしてしまう。これにより$\Delta t$を小さくとらずとも陽解法の解は\refeqn{142}を満たし安定となることができる。
\vs{0}
　ただしこの場合、音速を本来の物理的な値よりも小さく見積もっていることになってしまい、これは実質的に圧縮性を大きく見積もっていることに相当する。すなわち非圧縮条件が崩れ解が非圧縮解からズレることになる。よって陽解法による計算コストの削減は「数値不安定性」か「非圧縮条件からのずれ」のいずれかを必ず代償として受け入れなければならないことを意味する。

\section{実装結果}
　実際にMPS法により水柱崩壊を計算した動画を\href{https://shota-mitamura.com/works/ryuushi/ryuushhou/}{動画ページ}に示している。

\subsection{半陰解法}
　まず初めに半陰解法による二次元水柱崩壊の計算を行った。（\reffig{7}参照。）

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/haninn.png}
  \caption{半陰解法を用いた二次元水柱崩壊の数値計算結果。色は速度場を示している。}
  \label{fig:7}
\end{figure}
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/nensei.png}
  \caption{高粘性流体における半陰解法を用いた二次元水柱崩壊の数値計算結果。色は圧力場を示している。}
  \label{fig:8}
\end{figure}
\vs{0}
　その結果、粒子法による計算では波面の波飛沫が非常にダイナミックに再現できること、８ノード48コアのスパコン計算を用いても2分近くかかり、適切な近傍粒子探査による計算量削減を実施しないと粒子法の計算は非常に重くなってしまうことなどがわかった。
\vs{0}
　また\reffig{8}に示すように高い粘性を持つ流体でも水柱崩壊の計算を行った結果、粒子法は流体の粘性をよく捉えることができることがわかった。（物理的な正確性は別として、流体の「粘り気をよく表現することがわかった。」）

\subsection{陽解法}
　次に、陽解法による二次元水柱崩壊の計算を行った。
\vs{0}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/hassan.png}
  \caption{陽解法を用いた二次元水柱崩壊の計算結果。数値不安定により発散してしまっている。}
  \label{fig:9}
\end{figure}
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/asshuku.png}
  \caption{安定性条件を満たすように時間ステップ幅を小さくして計算した図。流体粒子同士が密に重なり合って圧縮されてしまっていることがわかる。}
  \label{fig:10}
\end{figure}

　その結果、\reffig{9}に示すように数値不安定により計算が爆発してしまうことがわかった。より詳しく動画を見るとわかるが、これは流体の重力によって加速し続け速度が一定値を超えた瞬間起こっていることがわかる。すなわち先に示した不安定条件が効いていることがわかる。また、この際の計算時間は1分半程度に抑えられており先ほどよりも早くなっていた。よって陽解法が計算量の削減と数値安定性のトレードオフの関係にあることを実際に確認できた。
\vs{0}
　また、時間ステップ幅を小さく取ることにより陽解法における数値不安定性を抑えた実験では流体が圧縮や膨張をしてしまうことがわかった。（\reffig{10}参照。）すなわち、陽解法において安定性を実現しようとすると今度は物理的な精度がトレードオフされてしまうということがわかり、先の章で見た通り陽解法の計算量削減が必ず何がしかの悪影響とトレードオフになってしまっていることがわかった。

\subsection{３次元の水柱落下}
　最後に三次元の水柱落下の計算を行った。\reffig{11}、\reffig{12}参照。

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/3d1.png}
  \caption{初め空中にあった水柱が浅い水槽に落下する様子。（１）}
  \label{fig:11}
\end{figure}
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.65\columnwidth]{./image_tex/3d2.png}
  \caption{初め空中にあった水柱が浅い水槽に落下する様子。（２）}
  \label{fig:12}
\end{figure}

\vs{0}
　特に今回、水槽を記述するための粒子数が次元の増加とともに大きく増えること、それに伴い２次元から３次元に変わると爆発的に計算時間が伸びることなどがわかった。先の二次元の水柱崩壊では2分程度で20000ステップの計算が終わったが、３次元の場合においては20000ステップの計算の終了までにおおよそ2時間近い時間がかかることがわかった。\footnote{
  今回は計算資源の都合上15分で計算できる最初の数ステップだけ計算を行った。
}

\section{結果とまとめ}
　今回、粒子法の中でも特にMPS法を実装し、二次元と三次元での水柱崩壊・水柱落下のシミュレーションを行った。その結果として、

\begin{itemize}
  \item 粒子法は波面の運動をダイナミックに表現することができる。
  \item 粒子法は粒子数の増大とともに計算コストが重くなってしまいがちであり、近傍粒子探索や並列化を効率よく用いる必要がある。
  \item 陽解法は計算コストが小さくて済む代わりに不安定であり、その不安定性を抑えると今度は非圧縮性の喪失という形で別の問題が現れてしまう。
\end{itemize}

といったことがわかった。特に陽解法の計算コスト削減と諸々のトレードオフは実験前から予想できていたことであり、それを実際に確認できたことは非常に良い成果であった。
\vs{0}
　また陽解法による計算コストの削減については、今回はポアソン方程式の求解の$\mathcal{O}(kN^2)$が陽解法の$\mathcal{O}(N)$に置き換わった程度であり、計算が全体として$\mathcal{O}(N^2)$に支配されていた為体感できる計算量削減は小さかった。しかし、適切な近傍粒子探索を実装して半陰解法の支配的な計算オーダーを$\mathcal{O}(kN)$にすることができれば、陽解法によって計算オーダーが$\mathcal{O}(N)$に落ちた時の変化をより強く感じることができると期待される。

\section{おわりに}
　今回の粒子法の実装は、サンプルコード付きの教科書があったとはいえそれなりに骨の折れる作業であった。特にMPIで並列化する際に計算のどのタイミングでデータを送受信しなければならないかを考えるのが非常に難しかった。結果的にコードを完成させるまでに（休み休みであったとはいえ）半年ほどかかってしまったが、最終的にこのように授業の課題に還元することができ、良い思い出になったと感じている。\footnote{
  期日を授業日の木曜だと勘違いしており、気付いてから慌ててレポートを始めた結果１時間ほど遅れての提出になってしまいました。本当に申し訳ございません。
}

\section*{参考文献}

\begin{thebibliography}{9}
\bibitem{koshizuka2014}
越塚誠一，柴田和也，室谷浩平，
『粒子法入門：流体シミュレーションの基礎から並列計算と可視化まで』，
丸善出版，2014．

\bibitem{intro}
\url{https://www.jstage.jst.go.jp/article/jjws/80/6/80_534/_pdf}

\bibitem{intro2}
\url{https://www.jepoc.or.jp/tecinfo/library.php?_w=Library&_x=detail&library_id=436}

\end{thebibliography}
\end{document}



