#はじめに
##何をするか
二次元非圧縮 Navier-Stokes 方程式解いてこんなアニメ作りたい…二次元非圧縮 Navier-Stokes 方程式解いてこんなアニメ作りたくない？？？
ということで作りましょう！！！
苦労もしましたが非常に勉強になり、とても楽しかったので「数値計算に興味あって多少は勉強したけど何をやったら分からない…」という過去の自分と同じ状態の人に是非挑戦してもらいたいと思い、記事としてまとめます。
##どうやるか
[使った言語]
C++ で計算し、Python の matplotlib で描画します。
計算手法としては、計算の流れが明確で、各パートごとに構築したコードがそのまま使える Fractional Step 法を採用しました。
また、東工大さんの資料を非常に参考にさせていただきました。
##この記事の目標とスタンス
Navier-Stokes 方程式は与えられたものとして、それをどう計算するかにフォーカスした記事にします。「簡単な数値計算はやったことがある人が一つ一つの計算過程を抵抗なく理解できる」程度の詳しさを目指し、流体力学の知識が無い人も計算を通して流体の振る舞いを実感できるような記事にできたら良いなと思います。
（できれば maplotlib でのアニメーション作成も解説したいです）
また、全部書いてから投稿だと一生投稿しなそうだったので、逐次更新とさせていただきます。初回は計算フローの大枠の解説とコードの提示までとします。
##エイリアスについて
まず、vector型の宣言時には、以下のエイリアス宣言を使っています（横着）
vector_macro.cpp
using vd = std::vector<double>;
using vvd = std::vector<vd>;
#計算の流れ
解くべき方程式がどんなものなのかを確認して、計算手順を説明します。
##どんな方程式を解くのか
さて、今回解くべき微分方程式は二次元 Navier-Stokes 方程式を無次元化し、非圧縮の条件を課した次の連立非線形微分方程式です。
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \\
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{\partial p}{\partial x} + \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -\frac{\partial p}{\partial y} + \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right)
こんなのどこから手を付けて良いかも分からん………って感じなので、分解して一つ一つ見ていきます。
第一式は連続の式と言われ、ここに非圧縮の条件が入っています。これと圧力項を忘れて
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right)
としたものを二次元 Burgers 方程式と言い、衝撃波を記述するモデルになっているようです。これを解くことが中間目標になります。
これでもまだ馴染みがないのでどんどん分解します。u, vについての非線形方程式は嫌なのであるスカラー場 f についての方程式
\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} + v\frac{\partial f}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \right)
を考えると、これは二次元移流拡散方程式と言います。
さらに分解して
\frac{\partial f}{\partial t} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \right)
\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} = 0
などとすれば、やっと馴染みのある式になったのではないでしょうか。
第一式は拡散方程式、第二式は移流方程式でこれらは初期条件によっては容易に解析的に解くことができます。
そこで、最初はこれらの方程式を解くコードを作るところからスタートします。そして、Fractional Step 法ではここで作ったコードがそのまま本体のコードで使えます。
##計算手順
###Burger 方程式の計算
まず連続の式と圧力項を忘れたBurgers 方程式：
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right)
を解きます。
これは、二次元拡散方程式とx、y 方向の１次元移流方程式：
\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} = 0\\
\frac{\partial f}{\partial t} + v\frac{\partial f}{\partial y} = 0\\
\frac{\partial f}{\partial t} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \right)\\
をタイムステップ幅 dt で時間前進を用いて解く関数がそれぞれ
advection_diffusion.cpp
void x_advection(vvd &f, vvd &fn, vvd &u, double dt);
void y_advection(vvd &f, vvd &fn, vvd &v, double dt);
void diffusion(vvd &f, vvd &fn, double dt);
などと出来ている場合、Fractional Step 法で
- u, v の x 方向への移流を計算
- u, v の y 方向への移流を計算
- 拡散項の計算
と分解して順番に解くことができます。
以下に擬似コードを示します。
do{
{//x方向に移流
x_advection(u, un, u, dt);
x_advection(v, vn, u, dt);
boundary(un, vn);
u = un; v = vn;
}
{//y方向に移流
y_advection(u, un, v, dt);
y_advection(v, vn, v, dt);
boundary(un, vn);
u = un; v = vn;
}
{//拡散
diffusion(u, un, dt);
diffusion(v, vn, dt);
boundary(un, vn);
u = un; v = vn;
}
t += dt;
} while ();
###速度の修正
時間 n ステップまで計算が終了しているときに Burgers 方程式に時間前進差分を適用することで得られる速度を $ u^* $, $ v^* $とすると、これらは当然（一般には）連続の式は満たしていません。
\frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \neq 0
が、恐らくこれらは（十分小さい時間ステップ幅 dt を設定すれば）求めたい真の速度 $ u^{n+1} $, $ v^{n+1} $ にある程度近いはずです。そこでこれを修正することを考えます。
$ u^* $, $ v^* $ が満たす式：
\frac{u^* - u^n}{\Delta t} + u^n\frac{\partial u^n}{\partial x} + v^n\frac{\partial u^n}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2} \right)\\
\frac{v^* - v^n}{\Delta t} + u^n\frac{\partial v^n}{\partial x} + v^n\frac{\partial v^n}{\partial y} = \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 v^n}{\partial x^2} + \frac{\partial^2 v^n}{\partial y^2} \right)
と、$ u^{n+1} $, $ v^{n+1} $ が満たすべき式：
\frac{u^{n+1} - u^n}{\Delta t} + u^n\frac{\partial u^n}{\partial x} + v^n\frac{\partial u^n}{\partial y} = -\frac{\partial p^n}{\partial x} + \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2} \right)\\
\frac{v^{n+1} - v^n}{\Delta t} + u^n\frac{\partial v^n}{\partial x} + v^n\frac{\partial v^n}{\partial y} = -\frac{\partial p^n}{\partial y} + \frac{1}{\mathrm{Re}}\left( \frac{\partial^2 v^n}{\partial x^2} + \frac{\partial^2 v^n}{\partial y^2} \right)
との差を取ることで
u^{n+1} = u^* - \frac{\partial p^n}{\partial x}\Delta t\\
v^{n+1} = v^* - \frac{\partial p^n}{\partial y}\Delta t
を得ます。したがって圧力場 $p^n$ が分かれば真の速度が分かることになります。
この $p^n$ を決定するために、上式を $u^{n+1}$, $v^{n+1}$ に関する連続の式に代入します。
\frac{\partial^2 p^n}{\partial x^2} + \frac{\partial^2 p^n}{\partial y^2} = \frac{1}{\Delta t}\left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right) \equiv s(x, y)
右辺のソース項は既知の量なので、この圧力に関する poisson 方程式を解けば速度の補正量が分かり真の速度が計算できる、というわけです。
「とりあえず近い解を持ってきて、それが連続の式を満たすように自己矛盾のない圧力場を持ってきて補正する」という感じです。
物理的には「連続の式を満たさない部分は圧力に押し付ける」ようなイメージだと解釈しています。
この連続の式の使い方、面白くないですか？？僕は非常に好きです
###計算手順まとめ
以上から、Fractional Step 法を用いて二次元非圧縮 Navier-Stokes 方程式を解く計算手順は
- 計算条件、初期状態などを決定
- Burgers 方程式を解いて $u^* $, $v^*$ を求める。
- x方向への移流の計算
- y方向への移流の計算
- 拡散項の計算
- Poisson 方程式を解いて速度を補正する
- $u^* $, $v^*$ からソース項 s(x, y) を計算
- poisson 方程式を解く
- 速度の補正
- u, v から rotation などを計算し出力
となります。したがって、後に示すコードの全貌のうち、本質的な部分のみを抜き出すと以下のようになります。
summury.cpp
#include <vector>
using vd = std::vector<double>;
using vvd = std::vector<vd>;
/******************************計算条件******************************/
const int nx = 200+5;
const int ny = 200+5;
const double Lx = 1.0;
const double Ly = 1.0;
const double dx = Lx/double(nx-5);
const double dy = Ly/double(ny-5);
// 計算終了時間
const double ENDTIME = 2.0;
/*********************************************************************/
/******************************関数******************************/
// 初期状態を決定
void initial(vvd &u, vvd &v);
// un, vn に境界条件を課す
void boundary(vvd &un, vvd &vn);
// Burgers eq を解く
void x_advection(vvd &f, vvd &fn, vvd &u, double dt);
void y_advection(vvd &f, vvd &fn, vvd &v, double dt);
void diffusion(vvd &f, vvd &fn, double dt);
// 仮の流速 u*, v* からpoisson 方程式の右辺 s を求める
void get_s(vvd &u, vvd &v, vvd &s, double dt);
// poisson eq を解く
int poisson2d(vvd &f, vvd &s);
// 圧力から速度を補正
void correction(vvd &u, vvd &v, vvd &p, double dt);
void output_all(/*u, v, rot など*/);
/****************************************************************/
int main(){
vvd u(ny,vd(nx, 0.0)), un(ny, vd(nx, 0.0));
vvd v(ny,vd(nx, 0.0)), vn(ny, vd(nx, 0.0));
vvd p(ny,vd(nx, 0.0)), s(ny, vd(nx, 0.0));
double t, dt; //dt は適宜決定
initial(u, v);
do{
{// x方向に移流
x_advection(u, un, u, dt);
x_advection(v, vn, u, dt);
boundary(un, vn);
u = un; v = vn;
}
{// y方向に移流
y_advection(u, un, v, dt);
y_advection(v, vn, v, dt);
boundary(un, vn);
u = un; v = vn;
}
{// 拡散
diffusion(u, un, dt);
diffusion(v, vn, dt);
boundary(un, vn);
u = un; v = vn;
}
{// poisson 方程式と結合して速度の修正
get_s(u, v, s, dt);
poisson2d(p, s);
correction(u, v, p, dt);
boundary(u, v);
}
t += dt;
output_all(/*u, v, rot など*/);
} while (t < ENDTIME);
return 0;
}
#各パートの解説
随所にアニメーションのコードの解説も挟みつつ
- 計算領域・境界条件
- 移流方程式
- 拡散方程式
- 移流拡散方程式
- Burgers 方程式
- Poisson 方程式
の順に解説します。（予定）
##計算領域・境界条件（更新予定）
##移流方程式（更新予定）
Cubic セミ・ラグランジュ法で解きます
##拡散方程式（更新予定）
時間前進差分、空間中心差分で解きます
##移流拡散方程式（更新予定）
Fractional Step 法で上記の二つを単純に組み合わせることで解けます
##Burgers 方程式（更新予定）
移流拡散方程式の引数を変えるだけで解けます
冒頭で示したアニメーションと同じ初期条件で、Burgers 方程式を解くと次のようになります。
非圧縮の条件が無いだけで全然違いますねえ（非常に面白い）
##Poisson 方程式（更新予定）
SOR 法で解きます
#Navier-Stokes 方程式
今まで作ったコードをもとに、Burgers 方程式と Poisson 方程式を結合した完全なコードを示します
##計算コード
NS.cpp
/* 最終更新：2020/02/18
[機能]
・2D Navier-Stokes 方程式を非圧縮条件のもとで Fractional Step 法で計算する。
・50 行目あたりのパラメータをいじって計算条件を設定できる
・67 行目あたりで出力する時間の間隔、長さを調整できる
・INITIAL の値を変えることで初期条件を選べる
・BOUNDARY で境界条件を選択できるが fixed は試作中
・IMCOMPRESSIBLE を 0 にすると 連続の式を切って比較することができる
[出力形式]
・相対パス "./data/condition.txt" 中に以下の形式で計算条件を出力
NX NY FRAMES
Lx Ly Re mu
ただし、NX, NY は出力する格子点の数で FRAMES は出力した回数
・"./data/u.txt" に以下の形式で u の計算結果を出力
time1
（時刻 time1 での配列 u ）
time2
（時刻 time2 での配列 u ）
（以下略）
・v, p, rot, div についても同様
*/
#define _USE_MATH_DEFINES
#include <cstdio>
#include <cmath>
#include <time.h>
#include <vector>
using vd = std::vector<double>;
using vvd = std::vector<vd>;
/******************************CONFIG******************************/
// 0:diagonal 1:random 2:left 3:sin 4:static
const int INITIAL = 0;
// 0:fixed, 1:periodic (fixedは試作中)
const int BOUNDARY = 1;
// 0:2D burgers, 1: imcompressible
const int IMCOMPRESSIBLE = 1;
/*******************************************************************/
/******************************計算条件******************************/
// nx, ny は100くらいなら1分程度で終わる 200 くらいが精度良さそうだけど7分くらいかかる
const int nx = 100+5;
const int ny = 100+5;
const double Lx = 1.0;
const double Ly = 1.0;
const double dx = Lx/double(nx-5);
const double dy = Ly/double(ny-5);
const double Re = 10000.0;
//計算の安定性を決めるファクターμ, mu > 0.25 だと計算が爆発する
const double mu = 0.20;
/*
以下出力枚数の調整用
時刻 0 ~ endtime の間のプロファイルを時間 DT ごとに出力させるための変数達
TIME に配列 (DT, 2DT, ..., ENDTIME) をsetする
*/
const double T_EPS = 1.0e-10;
const double DT = 0.02;
const double ENDTIME = 10.0;
vd TIME;
//出力時刻を set する関数
void TIME_set();
/*************************************************************************/
/*************************poisson eq 関連の関数・定数*************************/
const double OMEGA = 1.8; // sor法の加速係数
const double P_EPS = 1.0e-4;
// 仮の流速 u*, v* からpoisson 方程式の右辺 s を求める
void get_s(vvd &u, vvd &v, vvd &s, double dt);
// SOR法で1ステップだけf更新。誤差|f - fn|の最大値を返す
double sor(vvd &f, vvd &s);
// sor をたくさん呼び出して poisson 方程式を解く。返り値は反復回数
int poisson2d(vvd &f, vvd &s);
// 残差の最大値を計算。sor反復中に呼び出すことで精度の上昇を見れる。
double residual(vvd &f, vvd &s);
// 圧力から速度を補正
void correction(vvd &u, vvd &v, vvd &p, double dt);
/****************************************************************************/
/******************************burgersを解く関数******************************/
// 初期状態を決定
void initial(vvd &u, vvd &v);
// fn に境界条件を課す
void boundary(vvd &un, vvd &vn);
void diffusion(vvd &f, vvd &fn, double dt);
void x_advection(vvd &f, vvd &fn, vvd &u, double dt);
void y_advection(vvd &f, vvd &fn, vvd &v, double dt);
void rotation(vvd &u, vvd &v, vvd &rot);
void divergence(vvd &u, vvd &v, vvd &div);
/****************************************************************************/
/*********************************ファイル関連*********************************/
FILE *condition_fp = fopen("data/condition.txt","w");
FILE *u_fp = fopen("data/u.txt", "w");
FILE *v_fp = fopen("data/v.txt", "w");
FILE *div_fp = fopen("data/div.txt", "w");
FILE *rot_fp = fopen("data/rot.txt", "w");
FILE *p_fp = fopen("data/p.txt", "w");
// その時刻における f[jy][jx] の値をファイルにアウトプット
void output(vvd &f, double t, FILE *data_fp);
// u, v, p, rot, div をまとめてoutput
void output_all(vvd &u, vvd &v, vvd &p, vvd &rot, vvd &div, double t);
/*****************************************************************************/
int main(){
double dt, t = 0.0;
int ti = 0, icnt = 0; //TIMEのindex
int iteration;
vvd u(ny,vd(nx, 0.0)), un(ny, vd(nx, 0.0));
vvd v(ny,vd(nx, 0.0)), vn(ny, vd(nx, 0.0));
vvd rot(ny,vd(nx, 0.0)), div(ny, vd(nx, 0.0));
vvd p(ny,vd(nx, 0.0)), s(ny, vd(nx, 0.0));
clock_t start_t, end_t;
start_t = time(NULL);
initial(u, v);
boundary(u, v);
TIME_set();
dt = fmin(0.2*fmin(dx,dy)/1.0, mu*fmin(dx*dx, dy*dy)*Re);
printf("NX:%d NY:%d\nRe:%f mu:%f\n", nx, ny, Re, mu);
printf("dt:%.10f\n", dt);
output_all(u, v, p, rot, div, 0.0); icnt++;
do{
//printf("time:%f\n", t);
{// x方向に移流
x_advection(u, un, u, dt);
x_advection(v, vn, u, dt);
boundary(un, vn);
u = un; v = vn;
}
{// y方向に移流
y_advection(u, un, v, dt);
y_advection(v, vn, v, dt);
boundary(un, vn);
u = un; v = vn;
}
{// 拡散
diffusion(u, un, dt);
diffusion(v, vn, dt);
boundary(un, vn);
u = un; v = vn;
}
// poisson 方程式と結合して速度の修正（ここで非圧縮の条件が入る）
if(IMCOMPRESSIBLE == 1){
get_s(u, v, s, dt);
iteration = poisson2d(p, s);
if(iteration == -1){
printf("poisson does not converged when time = %f\n", t);
return 0;
}
correction(u, v, p, dt);
boundary(u, v);
}
t += dt;
// 全ての更新が終わったら出力
if(ti < TIME.size() && t > TIME[ti] - T_EPS){
output_all(u, v, p, rot, div, t);
ti++; icnt++;
}
} while (t < ENDTIME + dt);
printf("number of pictures:%d\n", icnt);
fprintf(condition_fp, "%d %d %d\n%f %f %f %f\n", nx-4, ny-4, icnt, Lx, Ly, Re, mu);
fclose(condition_fp);
fclose(u_fp);
fclose(v_fp);
fclose(div_fp);
fclose(rot_fp);
fclose(p_fp);
end_t = time(NULL);
printf("This calculatioin took %ld second \n", end_t - start_t);
return 0;
}
void initial(vvd &u, vvd &v){
if(INITIAL == 0){
for(int jy = 0; jy < ny; jy++) {
for(int jx = 0; jx < nx; jx++) {
if(0.3*nx < jx && jx < 0.7*nx && 0.3*ny < jy && jy < 0.7*ny){
u[jy][jx] = 1.0;
v[jy][jx] = 1.0;
}
}
}
}
if(INITIAL == 1){
//全部の点の50%くらいをランダムに選ぶ
for(int i = 0; i < nx*ny*0.5; i++) {
u[rand()%ny][rand()%nx] = 1.0;
v[rand()%ny][rand()%nx] = 1.0;
}
}
if(INITIAL == 2){
for(int jy = 0; jy < ny; jy++) {
for(int jx = 0; jx < nx; jx++) {
if(0.3*ny < jy && jy < 0.7*ny && 0.3*nx < jx && jx < 0.7*nx){
u[jy][jx] = -1.0;
v[jy][jx] = 0.0;
}
}
}
}
if(INITIAL == 3){
double x, y, kx = 2.0*M_PI, ky = 2.0*M_PI;
for(int jy = 0; jy < ny; jy++) {
for(int jx = 0; jx < nx; jx++) {
x = dx*(double)(jx-2);
y = dy*(double)(jy-2);
u[jy][jx] = -cos(kx*x)*sin(ky*y)/kx;
v[jy][jx] = sin(kx*x)*cos(ky*y)/ky;
x += 0.3; y+= 0.7;
u[jy][jx] = -0.6*cos(2.0*kx*x)*sin(2.0*ky*y)/kx;
v[jy][jx] = 0.6*sin(2.0*kx*x)*cos(2.0*ky*y)/ky;
}
}
}
u[ny/2][nx/2] = 1.0;
//check用
// for(int jy = 0; jy < ny; jy++) {
// for(int jx = 0; jx < nx; jx++) {
// if(jx < nx-1){
// printf("%.1f ", f[jy][jx]);
// }
// else{
// printf("%.1f\n", f[jy][jx]);
// }
// }
// }
return;
}
void output(vvd &f, double t, FILE *data_fp){
fprintf(data_fp,"%f\n", t);
//端っこの境界条件のためのダミーの格子点は出力しない
for(int jy = 2; jy <= ny-3; jy++) {
for(int jx = 2; jx <= nx-3; jx++) {
if(jx < nx-3){
fprintf(data_fp, "%f ", f[jy][jx]);
}
else{
fprintf(data_fp, "%f\n", f[jy][jx]);
}
}
}
}
void output_all(vvd &u, vvd &v, vvd &p, vvd &rot, vvd &div, double t){
divergence(u, v, div);
rotation(u, v, rot);
output(u, t, u_fp);
output(v, t, v_fp);
output(div, t, div_fp);
output(rot, t, rot_fp);
output(p, t, p_fp);
}
void diffusion(vvd &f, vvd &fn, double dt){
//境界は更新しない
for(int jy = 2; jy < ny-2; jy++) {
for(int jx = 2; jx < nx-2; jx++) {
fn[jy][jx] = f[jy][jx] + dt *
( (f[jy][jx+1] - 2.0*f[jy][jx] + f[jy][jx-1])/dx/dx
+ (f[jy+1][jx] - 2.0*f[jy][jx] + f[jy-1][jx])/dy/dy
)/Re;
}
}
return;
}
void boundary(vvd &un, vvd &vn){
if(BOUNDARY == 0){
for(int jy=0 ; jy < ny; jy++){
un[jy][nx-1] = un[jy][nx-3];// 流入境界を導入中
un[jy][nx-2] = un[jy][nx-3];
un[jy][0] = 1.0;
un[jy][1] = 1.0;
vn[jy][nx-1] = vn[jy][nx-3];
vn[jy][nx-2] = vn[jy][nx-3];
vn[jy][1] = 0.0;
vn[jy][0] = 0.0;
}
for(int jx=0 ; jx < nx; jx++){
un[ny-1][jx] = 1.0;
un[ny-2][jx] = 1.0;
un[1][jx] = 1.0;
un[0][jx] = 1.0;
vn[ny-1][jx] = 0.0;
vn[ny-2][jx] = 0.0;
vn[1][jx] = 0.0;
vn[0][jx] = 0.0;
}
}
if(BOUNDARY == 1){
for(int jy=0 ; jy < ny; jy++){
un[jy][nx-1] = un[jy][3];
un[jy][nx-2] = un[jy][2];
un[jy][1] = un[jy][nx-3];
un[jy][0] = un[jy][nx-4];
vn[jy][nx-1] = vn[jy][3];
vn[jy][nx-2] = vn[jy][2];
vn[jy][1] = vn[jy][nx-3];
vn[jy][0] = vn[jy][nx-4];
}
for(int jx=0 ; jx < nx; jx++){
un[ny-1][jx] = un[3][jx];
un[ny-2][jx] = un[2][jx];
un[1][jx] = un[ny-3][jx];
un[0][jx] = un[ny-4][jx];
vn[ny-1][jx] = vn[3][jx];
vn[ny-2][jx] = vn[2][jx];
vn[1][jx] = vn[ny-3][jx];
vn[0][jx] = vn[ny-4][jx];
}
}
//check 用
// printf("boundary check\n");
// for(int jy = 0; jy < ny; jy++) {
// for(int jx = 0; jx < nx; jx++) {
// if(jx < nx-1){
// printf("%.1f ", vn[jy][jx]);
// }
// else{
// printf("%.1f\n", vn[jy][jx]);
// }
// }
// }
return;
}
void x_advection(vvd &f, vvd &fn, vvd &u, double dt){
double a,b,c,z;
for (int jy = 2; jy < ny - 2; jy++){
for (int jx = 2; jx < nx - 2; jx++){
if (u[jy][jx] > 0.0){
a = (f[jy][jx+1] - 3.0*f[jy][jx] + 3.0*f[jy][jx-1] - f[jy][jx-2]) / (6.0*dx*dx*dx);
b = (f[jy][jx+1] - 2.0*f[jy][jx] + f[jy][jx-1]) / (2.0*dx*dx);
c = (2.0*f[jy][jx+1] + 3.0*f[jy][jx] - 6.0*f[jy][jx-1] + f[jy][jx-2]) / (6.0*dx);
}
else{
a = (f[jy][jx+2] - 3.0*f[jy][jx+1] + 3.0*f[jy][jx] - f[jy][jx-1]) / (6.0*dx*dx*dx);
b = (f[jy][jx+1] - 2.0*f[jy][jx] + f[jy][jx-1]) / (2.0*dx*dx);
c = (-f[jy][jx+2] + 6.0*f[jy][jx+1] - 3.0*f[jy][jx] - 2.0*f[jy][jx-1]) / (6.0*dx);
}
z = -u[jy][jx]*dt;
fn[jy][jx] = a*z*z*z + b*z*z + c*z + f[jy][jx];
}
}
}
void y_advection(vvd &f, vvd &fn, vvd &v, double dt){
double a,b,c,z;
for (int jy = 2; jy < ny - 2; jy++){
for (int jx = 2; jx < nx - 2; jx++){
if (v[jy][jx] > 0.0){
a = (f[jy+1][jx] - 3.0*f[jy][jx] + 3.0*f[jy-1][jx] - f[jy-2][jx]) / (6.0*dy*dy*dy);
b = (f[jy+1][jx] - 2.0*f[jy][jx] + f[jy-1][jx]) / (2.0*dy*dy);
c = (2.0*f[jy+1][jx]+ 3.0*f[jy][jx] - 6.0*f[jy-1][jx] + f[jy-2][jx]) / (6.0*dy);
}
else{
a = (f[jy+2][jx] - 3.0*f[jy+1][jx] + 3.0*f[jy][jx] - f[jy-1][jx]) / (6.0*dy*dy*dy);
b = (f[jy+1][jx] - 2.0*f[jy][jx] + f[jy-1][jx]) / (2.0*dy*dy);
c = (-f[jy+2][jx] + 6.0*f[jy+1][jx] - 3.0*f[jy][jx] - 2.0*f[jy-1][jx]) / (6.0*dy);
}
z = -v[jy][jx]*dt;
fn[jy][jx] = a*z*z*z + b*z*z + c*z + f[jy][jx];
}
}
}
void TIME_set(){
double tmp = DT;
while(tmp < ENDTIME + T_EPS){
TIME.push_back(tmp);
tmp += DT;
}
}
void rotation(vvd &u, vvd &v, vvd &rot){
for(int jy = 1; jy < ny-1; jy++) {
for(int jx = 1; jx < nx-1; jx++) {
rot[jy][jx] = 0.5*(v[jy+1][jx+1] - v[jy+1][jx]+ v[jy][jx+1] - v[jy][jx])/dx
- 0.5*(u[jy+1][jx+1] - u[jy][jx+1] + u[jy+1][jx] - u[jy][jx])/dy;
}
}
}
void divergence(vvd &u, vvd &v, vvd &div){
for (int jy = 1; jy < ny-1; jy++){
for (int jx = 1; jx < nx-1; jx++){
div[jy][jx] = 0.5*(u[jy+1][jx+1] - u[jy+1][jx] + u[jy][jx+1] - u[jy][jx])/dx
+ 0.5*(v[jy+1][jx+1] - v[jy][jx+1] + v[jy+1][jx] - v[jy][jx])/dy;
}
}
}
void get_s(vvd &u, vvd &v, vvd &s, double dt){
for (int jy = 2; jy < ny-2; jy++){
for (int jx = 2; jx < nx-2; jx++){
s[jy][jx] = 0.5*(u[jy+1][jx+1] - u[jy+1][jx] + u[jy][jx+1] - u[jy][jx])/dx
+ 0.5*(v[jy+1][jx+1] - v[jy][jx+1] + v[jy+1][jx] - v[jy][jx])/dy;
s[jy][jx] /= dt;
}
}
}
double residual(vvd &f, vvd &s){
double res, rmax = 0.0;
//各格子点の(d^2f/dx^2 + d^2f/dy^2 - s)を計算
for (int jy = 1; jy < ny - 1; jy++){
for (int jx = 1; jx < nx - 1; jx++){
res = (f[jy][jx + 1] - 2.0 * f[jy][jx] + f[jy][jx - 1]) / (dx * dx)
+ (f[jy + 1][jx] - 2.0 * f[jy][jx] + f[jy - 1][jx]) / (dy * dy)
- s[jy][jx];
rmax = fmax(res, rmax);
}
}
return rmax;
}
double sor(vvd &f, vvd &s){
double fn, err = 0.0;
for (int jy = 1; jy < ny - 1; jy++){
for (int jx = 1; jx < nx - 1; jx++){
fn = ((f[jy][jx + 1] + f[jy][jx - 1]) / (dx * dx)
+ (f[jy + 1][jx] + f[jy - 1][jx]) / (dy * dy) - s[jy][jx])
* 0.5 * dx * dx * dy * dy / (dx * dx + dy * dy);
err = fmax(fabs(fn - f[jy][jx]), err);
f[jy][jx] = (1.0 - OMEGA) * f[jy][jx] + OMEGA * fn;
}
}
return err;
}
int poisson2d(vvd &f, vvd &s){
int icnt = 0, imax = 99999;
double err;
// 反復
while(icnt++ < imax){
//fの更新と更新前後の差の最大値確認
err = sor(f, s);
//更新してもほとんど変わらないようなら終了
if(P_EPS > err) return icnt;
}
// imax 回反復しても収束しないなら-1を返して終了
return -1;
}
void correction(vvd &u, vvd &v, vvd &p, double dt){
for (int j = 2; j < ny-2; j++){
for (int i = 2; i < nx-2; i++){
u[j][i] += -0.5*(p[j][i] - p[j][i-1] + p[j-1][i] - p[j-1][i-1])/dx*dt;
v[j][i] += -0.5*(p[j][i] - p[j-1][i] + p[j][i-1] - p[j-1][i-1])/dy*dt;
}
}
}
##アニメーション作成
animation.py
#最終更新： 2020/02/18
'''
[機能]
・相対パス "./data/" の中にあるファイルから計算データを読み込みアニメーションを作成
・20 行目あたりで好きな表示方法を選べる
・MP4 と GIF の場合は相対パス "./animes/" の中にアニメーションを出力
・P_DIV を 1 にすれば圧力と発散も表示できる
[注意]
・colorbar のスケールは 62 行目あたりの vmax で無理矢理合わせているから、初期条件によっては不適切なスケールになる
'''
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.animation import FuncAnimation
from matplotlib.colors import Normalize # Normalizeをimport
from IPython.display import HTML
###################### CONFIG ##################
# 動画保存方法の config, 使いたい保存方法を１にする
MP4 = 0
GIF = 1
HTML_SHOW = 0
PLT = 0
# 圧力と発散を表示させるかどうか、1なら表示
P_DIV = 0
###################### CONFIG ##################
################### PARAMETER ##################
TITLE = '2D Navier-Stokes equation'
cfp = open('data/condition.txt')
NX, NY, FRAMES = map(int, cfp.readline().split())
Lx, Ly, Re, mu = map(float, cfp.readline().split())
cfp.close
# quiver で描画するときに領域を何区間に分割するか（あんまり多いと見辛いので
# NX-1. Ny-1 の約数じゃないと壊れる
QNX = 20
QNY = 20 * int(Ly/Lx)
################### PARAMETER ##################
####################################描画のための関数####################################
def get_data(f):
time = float(f.readline())
data=[]
for _ in range(NY):
tmp = list(map(float, f.readline().split()))
data.append(tmp)
return data, time
def init_im(ax, f, title):
data, _ = get_data(f)
ax.tick_params(labelsize=14)
ax.set_title(title, fontsize=20)
if(title == "divergence"):
vmax = 0.1
elif(title == "rotation"):
vmax = 20
elif(title == "pressure"):
vmax = 0.5
im = ax.imshow(data, extent=(0,Lx,0,Ly), origin="lower", animated=True,
cmap='jet', norm=Normalize(vmin=-vmax, vmax=vmax))
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = fig.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=12)
return im
def im_reset(im, f):
data, time = get_data(f)
im.set_data(data)
return time
def get_UV():
_ = float(f_u.readline())
U = []
for _ in range(NY):
tmp = list(map(float, f_u.readline().split()))
tmp = tmp[::NX//QNX]
U.append(tmp)
U = U[::NY//QNY]
_ = float(f_v.readline())
V = []
for _ in range(NY):
tmp = list(map(float, f_v.readline().split()))
tmp = tmp[::NX//QNX]
V.append(tmp)
V = V[::NY//QNY]
return U, V
def init_quiver(ax):
X, Y = np.meshgrid(np.linspace(0, Lx, QNX+1), np.linspace(0, Ly, QNY+1))
U, V = get_UV()
C = np.hypot(U, V)
ax.set_title("flow vector", fontsize=20)
ax.set_aspect("equal")
ax.tick_params(labelsize=14)
im = ax.quiver(X, Y, U, V, C, scale_units="width", cmap="jet")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = fig.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=12)
return im
def quiver_reset(im):
U, V = get_UV()
C = np.hypot(U, V)
im.set_UVC(U, V, C)
######################################################################################
f_u = open('data/u.txt')
f_v = open('data/v.txt')
f_div = open('data/div.txt')
f_rot = open('data/rot.txt')
f_p = open('data/p.txt')
if(P_DIV == 0):
fig = plt.figure(figsize=(10, 6))
else:
fig = plt.figure(figsize=(10, 10))
fig.suptitle(TITLE, fontsize=20)
#fig.subplots_adjust(bottom=0.10)
if(P_DIV == 1):
ax_q = fig.add_subplot(221)
ax_rot = fig.add_subplot(222)
ax_p = fig.add_subplot(223)
ax_div = fig.add_subplot(224)
else:
ax_q = fig.add_subplot(121)
ax_rot = fig.add_subplot(122)
fig.text(0, 0.01, "Re="+str(Re),
backgroundcolor="black",color="white", size=20)
time_text = fig.text(0, 0.99, '', size=20, color="white", horizontalalignment='left',
verticalalignment='top', backgroundcolor='black')
################### アニメの初期画像 ###################
# im_uの初期化
im_q = init_quiver(ax_q)
im_rot = init_im(ax_rot, f_rot, "rotation")
if(P_DIV == 1):
im_p = init_im(ax_p, f_p, "pressure")
im_div = init_im(ax_div, f_div, "divergence")
time_text.set_text("time = 0.00")
######################################################
plt.tight_layout()
# 更新用関数
def update(cnt):
# 最初は更新せずに初期の画像を出力させる
if cnt == 0:
return
print("cnt:",cnt)
# 変更箇所だけreset
quiver_reset(im_q)
time = im_reset(im_rot, f_rot)
if(P_DIV == 1):
im_reset(im_p, f_p)
im_reset(im_div, f_div)
time_text.set_text("time = %.2f"%time)
ani = FuncAnimation(fig, update, repeat=True, interval=20, frames=FRAMES)
# オプション frames 何枚の画像か，つまり関数を呼び出す回数
# mp4 ファイルとして保存
if(MP4 == 1):
ani.save("animes/imcompressible_anime.mp4", writer="ffmpeg", fps=20)
# gif ファイルとして保存
if(GIF == 1):
ani.save("animes/imcompressible_anime.gif", writer="pillow", fps=20)
# HTML上で表示
if(HTML_SHOW == 1):
HTML(ani.to_jshtml())
if(PLT == 1):
plt.show()
f_u.close
f_v.close
f_div.close
f_rot.close
f_p.close
#おわりに
数値計算は楽しい✌('ω'✌ )三✌('ω')✌三( ✌'ω')✌
いっぱい手を動かそう！
#謝辞
計算テーマの構想を一緒に考え、アニメーションのコードの基盤を作ってくれた、ど(^-^)/君とmonmon君に感謝いたします。