Update doc part 1 and added solvers for velocity using sine transform.

docs_part1/vanilaCavityFlow_JP.md

@@ -115,7 +115,7 @@ CFD に関わる方にとっては基本的なことかと思いますが、自
 この手法自体は Harlow and Welch (1965) [2] や Chorin (1968) [3] などで提案されたと紹介されることが多く、MAC (Marker and Cell method) 法とも呼ばれます。圧力項を残したバージョンは SMAC (Simplified MAC) とか。この辺の歴史や呼び名は正直はっきり理解できていないのですが、詳しい方コメントください！時間発展の高精度化・安定化を図った Kim and Moin (1985) [4] の手法も実用的です。
 
 
-## 圧力の境界条件
+## 仮の速度場の境界条件
 
 
 ちなみにこの手法は 
@@ -131,7 +131,7 @@ CFD に関わる方にとっては基本的なことかと思いますが、自
 
 
 
-個人的には Perot (1993) [5] や Chang (2002) [7] で紹介されている、離散化した Navier-Stokes 方程式の LU 分解をベースにした議論が気に入っているので、これをここでは採用します。**圧力の境界条件は不要！**と鮮やかな見解を提示しています。
+個人的には Perot (1993) [5] や Chang (2002) [7] で紹介されている、離散化した Navier-Stokes 方程式の LU 分解をベースにした議論が気に入っているので、これをここでは採用します。**仮の速度場の境界条件なんて考える必要もない！**と鮮やかな見解を提示しています。
 
 
 # 離散化手法

functions/getIntermediateU_xyRK3_dst.m

@@ -0,0 +1,30 @@
+% Copyright 2020 The MathWorks, Inc.
+function xu = getIntermediateU_xyRK3_dst(u, b, dt, Re, nx, ny, dx, dy, id)
+
+kx = [1:nx-1]';
+ax = pi*kx/(nx);
+mwx = 2*(cos(ax)-1)/dx^2;% DST-I
+
+ky = [1:ny]';
+ay = pi*ky/(ny);
+mwy = 2*(cos(ay)-1)/dy^2; % DST-II
+
+mw = mwx+mwy'; % Modified Wavenumber
+
+% 各 column 毎に変換されるので、
+% column 方向（x 方向に）dst1, 転置して
+% y 方向に dst2、転置し返して元に戻す処理
+rhshat = mydst2(mydst1(b)')';
+
+switch id
+    case 1
+        uhat = rhshat./(1-(4/15)*dt/(Re)*mw);
+    case 2
+        uhat = rhshat./(1-(1/15)*dt/(Re)*mw);
+    case 3
+        uhat = rhshat./(1-(1/6)*dt/(Re)*mw);
+end
+
+xu = mydst1(mydst3(uhat')');
+
+end

functions/getIntermediateV_xyRK3_dst.m

@@ -0,0 +1,32 @@
+% Copyright 2020 The MathWorks, Inc.
+function xv = getIntermediateV_xyRK3_dst(v, b, dt, Re, nx, ny, dx, dy,id)
+
+
+kx = [1:nx]';
+ax = pi*kx/(nx);
+mwx = 2*(cos(ax)-1)/dx^2;% DST-I
+
+ky = [1:ny-1]';
+ay = pi*ky/(ny);
+mwy = 2*(cos(ay)-1)/dy^2; % DST-II
+
+mw = mwx+mwy'; % Modified Wavenumber
+
+% 各 column 毎に変換されるので、
+% column 方向（x 方向に）dst2, 転置して
+% y 方向に dst1、転置し返して元に戻す処理
+rhshat = mydst1(mydst2(b)')';
+
+switch id
+    case 1
+        vhat = rhshat./(1-(4/15)*dt/(Re)*mw);
+    case 2
+        vhat = rhshat./(1-(1/15)*dt/(Re)*mw);
+    case 3
+        vhat = rhshat./(1-(1/6)*dt/(Re)*mw);
+end
+
+xv = mydst3(mydst1(vhat')');
+
+end
+

functions/updateVelocityField_RK3.m

@@ -81,10 +81,12 @@
 
 % Implicit treatment for xy direction
 b = u(2:end-1,2:end-1) - dt*(gamma(id)*Nu + zeta(id)*Nu_old - alpha(id)/(Re)*(Lux+Luy+Lubc));
-xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id);
+% xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id);
+xu = getIntermediateU_xyRK3_dst(u, b, dt, Re, nx, ny, dx, dy, id);
 
 b = v(2:end-1,2:end-1) - dt*(gamma(id)*Nv + zeta(id)*Nv_old - alpha(id)/(Re)*(Lvx+Lvy+Lvbc));
-xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy, id);
+% xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy, id);
+xv = getIntermediateV_xyRK3_dst(v, b, dt, Re, nx, ny, dx, dy, id);
 
 u(2:end-1,2:end-1) = xu;
 v(2:end-1,2:end-1) = xv;