diff --git a/docs_part2/vanilaCavityFlowImplicit_JP.mlx b/docs_part2/vanilaCavityFlowImplicit_JP.mlx index 2d0fd79..9ae98d2 100644 Binary files a/docs_part2/vanilaCavityFlowImplicit_JP.mlx and b/docs_part2/vanilaCavityFlowImplicit_JP.mlx differ diff --git a/docs_part4/NSSolverValidation.md b/docs_part4/NSSolverValidation.md new file mode 100644 index 0000000..6c61be0 --- /dev/null +++ b/docs_part4/NSSolverValidation.md @@ -0,0 +1,454 @@ +# 非圧縮性 Navier-Stokes 方程式の数値解法４：ソルバー精度検証 +# はじめに + + +Navier-Stokes 方程式を数値的に解くシリーズ、第４回目です。 + + + + - [非圧縮性 Navier-Stokes 方程式の数値解法１：導入編](https://qiita.com/eigs/items/628a2aeb3cb9eef91093) + - [非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法](https://qiita.com/eigs/items/5a62d9aff7d83af55099) + - [非圧縮性 Navier-Stokes 方程式の数値解法３：陰解法の解き方比較](https://qiita.com/eigs/items/81ceaf8d0abc426d29d7) + - **非圧縮性 Navier-Stokes 方程式の数値解法４：ソルバー精度検証（この記事）** + + + +ここまで非圧縮性 Navier-Stokes 方程式を解くための手法・部品についていろいろ紹介してきましたが、そもそも Navier-Stokes 方程式が解けているのか？という点について（ようやく・・）検証してみます。 + + + +# 本記事の内容 + + +ある外力項を加えた Navier-Stokes 方程式の数値解と解析解と比較することで、 + + + + - 時間積分の精度（1次精度） + - 空間微分の精度（2次精度） + + + +を確認します。以下結果ですが、、想定通りですね。 + + + + +![image_0.png](NSSolverValidation_images/image_0.png) + + +## 実行環境 + + - MATLAB R2020a + - Signal Processing Toolbox（ポワソン方程式の解法に `dct` を使う場合 ） + - Symbolic Math Toolbox（なくても大丈夫。あると嬉しい） + + +# Method of Manufactured Solution + + +Navier-Stokes 方程式が ”ある解” を満たすような外力項を追加で加えて数値解の誤差を検証します。まず外力項を加えた Navier-Stokes 方程式と質量保存の式を考えます。 + + + + + + + +ここで速度（）と圧力（）の式を代入すれば必要な外力項 も分かりますね。もちろん"ある解"は質量保存の式を満たすものある必要があります。 + + + + +解析解と数値解の誤差は以下の式で評価します。 + + + + + + + +この方法、学生時は Method of Manufactured Solution など呼んだりしていましたが、他の名前をご存知の方コメントください。 + + + +# 減衰する 2 次元渦 + + +ここではテストケースとして以下の時間とともに減衰する 2 次元渦を使います。（レイノルズ数）が小さい（粘性が高い）と減衰も早くなりますし、 が領域内の渦の数を決めます。ここでは領域を とします。 + + + + + + + + +![image_1.png](NSSolverValidation_images/image_1.png) + + + + +図：等高線図（速さ）と矢印（流れの方向） + + + + + +以下で実際に検証しますが、質量保存の式は満たしています。また、Navier-Stokes 方程式に代入すると外力項は以下のように求まります。 + + + + + + + + + と を追加して Navier-Stokes 方程式を解けば、上の と が求まるという話。そしてシミュレーション結果と比較して精度を検証するというのが狙いです。 + + + +## 外力項の確認 + + +インストール済みであれば Symbolic Math Toolbox で計算してみますが、ない場合は結果だけ使用します。 + + + +```matlab:Code +addpath('../functions/'); +clear, close all + +if license('checkout','symbolic_toolbox') % Symbolic Math Toolbox がある場合 + syms usol(x,y,t) vsol(x,y,t) psol(x,y,t) fu(x,y,t) fv(x,y,t) Re a + + % 解析解の定義 + usol(x,y,t) = -sin(a*x)*cos(a*y)*exp(-2*t/Re); + vsol(x,y,t) = cos(a*x)*sin(a*y)*exp(-2*t/Re); + psol(x,y,t) = 1/4*(cos(2*a*x) + sin(2*a*y))*exp(-4*t/Re); + + % x 成分（Navier-Stokes） + fu = diff(usol,t) + usol*diff(usol,x) + vsol*diff(usol,y) + diff(psol,x)... + - 1/Re*(diff(usol,x,2)+diff(usol,y,2)); + + % y 成分（Navier-Stokes） + fv = diff(vsol,t) + usol*diff(vsol,x) + vsol*diff(vsol,y) + diff(psol,y) ... + - 1/Re*(diff(vsol,x,2)+diff(vsol,y,2)); + + % 式を簡略化 + fu = simplify(fu); + fv = simplify(fv,'Steps',1000); % 手間がかかるのでより入念に簡略化 + + % 計算用に MATLAB 関数にしておきます。 + usolFunc = matlabFunction(usol); + vsolFunc = matlabFunction(vsol); + psolFunc = matlabFunction(psol); + fuFunc = matlabFunction(fu); + fvFunc = matlabFunction(fv); + +else % ない場合 + usolFunc = @(x,y,t,Re,a) -exp((t.*-2.0)./Re).*cos(a.*y).*sin(a.*x); + vsolFunc = @(x,y,t,Re,a) exp((t.*-2.0)./Re).*cos(a.*x).*sin(a.*y); + psolFunc = @(x,y,t,Re,a) exp((t.*-4.0)./Re).*(cos(a.*x.*2.0)./4.0+sin(a.*y.*2.0)./4.0); + fuFunc = @(x,y,t,Re,a) -(exp((t.*-2.0)./Re).*cos(a.*y).*sin(a.*x).*(a.^2.*2.0-2.0))./Re; + fvFunc = @(x,y,t,Re,a) (a.*exp((t.*-4.0)./Re).*(cos(a.*y.*2.0)+sin(a.*y.*2.0)))./2.0+(exp((t.*-2.0)./Re).*cos(a.*x).*sin(a.*y).*(a.^2-1.0).*2.0)./Re; +end +``` + +## `質量保存の式確認` + + +念のため質量保存の式も確認しておきます。大丈夫ですね。 + + + +```matlab:Code +if license('checkout','symbolic_toolbox') + % 質量保存の式確認 + divergence = diff(usol,x) + diff(vsol,y) +end +``` + +divergence(x, y, t) = + + +# コードの変更箇所 + + +[非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法](https://qiita.com/eigs/items/5a62d9aff7d83af55099) まではキャビティ流れを想定していましたので、境界条件として上辺の速度だけが非ゼロでした。ここではまた任意の外力項に加えて、任意の速度を境界条件に設定できるよう変更を加えています。 + + + + +詳細は `updateVelocityField_CNAB.m` や `updateVelocityField_RK3.m` に譲りますが、大きな変更箇所は以下の2点。 + + + + +```matlab:Code(Display) +Lubc(1,:) = velbc1.uLeft(2:end-1)/dx^2; +Lubc(end,:) = velbc1.uRight(2:end-1)/dx^2; +Lubc(:,1) = 2*velbc1.uBottom(2:end-1)/dy^2; +Lubc(:,end) = 2*velbc1.uTop(2:end-1)/dy^2; +``` + + + +1つ目は拡散項の陰解法を解くにあたっての境界条件、そしてここは 1 ステップ先の速度場を使用する必要があります。 + + + + + + + +の と の項です。 + + + + +```matlab:Code(Display) +% 3-1. get derivative for u +Nu = (uuce(2:end,:) - uuce(1:end-1,:))/dx; +Nu = Nu + (uvco(2:end-1,2:end) - uvco(2:end-1,1:end-1))/dy; +% forcing +Nu = Nu - uForce; + +% 3-2. get derivative for v +Nv = (vvce(:,2:end) - vvce(:,1:end-1))/dy; +Nv = Nv + (uvco(2:end,2:end-1) - uvco(1:end-1,2:end-1))/dx; +% forcing +Nv = Nv - vForce; +``` + + + +2つ目は外力項です。こちらは非線形項に加えておきます。 + + +# 精度検証の条件設定 + + +早速検証を行いましょう。上でも書きましたが、ここでは + + + +```matlab:Code +a = 2*pi; +Re = 100; +``` + + + +として計算します。 の値を大きくすると渦の数が増えて空間微分の誤差が増えそうな気がしますね。さて、空間微分と時間積分の精度それぞれ検証してみます。 + + + + +解析解との誤差を返す関数 `checkNSSolverError` としてまとめていますが、入力引数の数が多くて雑多な感じになっちゃいましたが、これはこれで実直で分かりやすいかもしれないのでこのままにしておきます。 + + + +```matlab:Code(Display) +N = 64, dt = 0.01, tEnd = 10; +checkNSSolverError(Re,a,N,dt,tEnd,fuFunc,fvFunc,usolFunc,vsolFunc,psolFunc,'RK3',true,true); +``` + + + +条件設定の時間ステップ `dt` で計算し `tEnd` での誤差を返します。その際に `fuFunc, fvFunc` などは上で計算した外力項、`usolFunc` などは解析解です。 + + + + +最後の `true, true` は可視化設定です（1つ目：Figure 作成、2つ目: GIF作成）。 + + + + +![image_2.png](NSSolverValidation_images/image_2.png) + + + + +といった具合で比較します。速さの等高線、矢印は速度方向を表します。渦が減衰している様子が分かります。 + + +# 空間微分精度検証 + + +まずは空間微分の精度。 + + + + +ここではできるだけ時間積分による誤差は小さくしておきたいので、3 段階のルンゲ・クッタ法を使用して `t = 0.01` における誤差をグリッド数を変えて（`Nx = Ny = 16` から `512` まで）確認します。確認のため 2 つの時間ステップサイズを使います。 + + + +```matlab:Code +tEnd = 0.01; + +Ns = [16,32,64,128,256,512]; +dts = [1e-3, 1e-4]; +L2error_spacial = zeros(length(Ns),length(dts)); + +for jj=1:length(dts) + dt = dts(jj); + for ii=1:length(Ns) + N = Ns(ii); + L2error_spacial(ii,jj) = checkNSSolverError(Re,a,N,dt,tEnd,fuFunc,fvFunc,... + usolFunc,vsolFunc,psolFunc,'RK3',false,false); + end +end +``` + + +```text:Output +(N,dt) = (16,0.001) error: 0.019005 at 0.01 CFL: 0.01631 +(N,dt) = (32,0.001) error: 0.0047481 at 0.01 CFL: 0.032155 +(N,dt) = (64,0.001) error: 0.0012257 at 0.01 CFL: 0.06408 +(N,dt) = (128,0.001) error: 0.00032862 at 0.01 CFL: 0.12806 +(N,dt) = (256,0.001) error: 0.00012523 at 0.01 CFL: 0.25609 +(N,dt) = (512,0.001) error: 9.9427e-05 at 0.01 CFL: 0.51217 +(N,dt) = (16,0.0001) error: 0.019003 at 0.0101 CFL: 0.0016296 +(N,dt) = (32,0.0001) error: 0.0047471 at 0.0101 CFL: 0.0032127 +(N,dt) = (64,0.0001) error: 0.001223 at 0.0101 CFL: 0.0064034 +(N,dt) = (128,0.0001) error: 0.00031469 at 0.0101 CFL: 0.012799 +(N,dt) = (256,0.0001) error: 7.9397e-05 at 0.0101 CFL: 0.025596 +(N,dt) = (512,0.0001) error: 2.2161e-05 at 0.0101 CFL: 0.051193 +``` + +## 何次精度だったかな？ + + +プロットしてみます。 + + + +```matlab:Code +%% ここからプロット +figure +loglog(Ns,L2error_spacial,'Marker','o','LineWidth',4); + +% 見栄えの微調整 +handle_axes = gca; +handle_axes.FontSize = 14; +box(handle_axes,'on'); +grid(handle_axes,'on'); + +% legend を作成 +legend({'dt = 1e-3','dt = 1e-4'},'Location','best'); +grid on +xlabel('Number of Grid') +ylabel('Relative error') +title('Error vs Number of Grid') +hold off + +annotation('line',[0.24 0.62],[0.64 0.25],'LineWidth',3) +annotation('textbox',[0.3 0.3 0.2 0.1],... + 'String',{'- dx^2'},... + 'FontSize',18,'EdgeColor',[1 1 1]); +``` + + +![figure_0.png](NSSolverValidation_images/figure_0.png) + + + +無事に 2 次精度が確認できますね。`dt = 1e-3 `の時グリッド数 `N = 256, 512` で誤差が多めに出ているのは、時間積分に由来する誤差の影響と見ています。 + + +# 時間積分精度検証 + + +上の結果から、時間積分の精度を検証する場合には、空間微分による誤差で埋もれてしまわないように `Nx = Ny = 512` で試すとよさそうですね。ただ、`dt = 1e-4` まで小さくすると空間微分誤差に支配されていそうです。 + + + + +この辺の考察をもとに `dt` を `dt = 1e-3` から `dt = 1e-4` の間で時間ステップサイズを変え、 `t = 0.01` における誤差を確認します。ルンゲ・クッタ法（RK3）だけでなく、クランクニコルソン・アダムスバッシュフォース法（CNAB）も試します。 + + + +```matlab:Code +N = 512; +tEnd = 0.01; +dts = 1e-3*[2,1,1/2,1/4]; +L2error_temporal = zeros(length(dts),2); + +for ii=1:length(dts) + dt = dts(ii); + L2error_temporal(ii,1) = checkNSSolverError(Re,a,N,dt,tEnd,fuFunc,fvFunc,... + usolFunc,vsolFunc,psolFunc,'CNAB',false,false); + L2error_temporal(ii,2) = checkNSSolverError(Re,a,N,dt,tEnd,fuFunc,fvFunc,... + usolFunc,vsolFunc,psolFunc,'RK3',false,false); +end +``` + + +```text:Output +(N,dt) = (512,0.002) error: 0.00090377 at 0.01 CFL: 1.0278 +(N,dt) = (512,0.002) error: 0.00019527 at 0.01 CFL: 1.025 +(N,dt) = (512,0.001) error: 0.00024862 at 0.01 CFL: 0.51276 +(N,dt) = (512,0.001) error: 9.9427e-05 at 0.01 CFL: 0.51217 +(N,dt) = (512,0.0005) error: 0.00012595 at 0.01 CFL: 0.25616 +(N,dt) = (512,0.0005) error: 5.3331e-05 at 0.01 CFL: 0.25602 +(N,dt) = (512,0.00025) error: 6.5312e-05 at 0.01 CFL: 0.12803 +(N,dt) = (512,0.00025) error: 3.1945e-05 at 0.01 CFL: 0.12799 +``` + + +## 何次精度だったかな？ + + +プロットしてみます。 + + + +```matlab:Code +%% ここからプロット +figure +loglog(dts,L2error_temporal,'Marker','o','LineWidth',4); + +% 見栄えの微調整 +handle_axes = gca; +handle_axes.FontSize = 14; +box(handle_axes,'on'); +grid(handle_axes,'on'); + +% legend を作成 +legend({'CNAB','RK3'},'Location','best'); +grid on +xlabel('Time step size') +ylabel('Relative Error') +title('Error vs Time step size') +hold off + +annotation('line',[0.32 0.75],[0.24 0.69],'LineWidth',3) +annotation('textbox',[0.5 0.35 0.2 0.1],... + 'String',{'dt'},... + 'FontSize',18,'EdgeColor',[1 1 1]); +xlim([1e-4,1e-2]) +ylim([1e-5,1e-3]) +``` + + +![figure_1.png](NSSolverValidation_images/figure_1.png) + + + +1 次精度ですね。同じ時間ステップサイズならルンゲ・クッタ法の方が誤差が少ない結果になりました。ただ 3 段階のルンゲ・クッタが 1 次精度というのは予定通りですが少し寂しい。 + + +# まとめ + + +Method of Manufactured Solution を使って空間微分の 2 次精度、時間積分の 1 次精度を確認しました。空間微分の精度を高次化するのは単純に手間ですし、ステンシルが広がることで境界の取り扱いや離散サイン変換が使えるか否かの再検証など必要。ただ、グリッド数を10倍にしても誤差が 1/100 にしかならないのは、3 次元への展開も考えると多少心もとない気もします。 + + + + +3 段階のルンゲ・クッタで 1 次精度！？という気はしますが、これは定式化の段階（LU分解）の段階で織り込み済み（参照：[非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法](https://qiita.com/eigs/items/5a62d9aff7d83af55099)）という理解です。こうすれば高次化できますよ、というヒントあればコメントお待ちしております。 + + diff --git a/docs_part4/NSSolverValidation.mlx b/docs_part4/NSSolverValidation.mlx new file mode 100644 index 0000000..d85e19e Binary files /dev/null and b/docs_part4/NSSolverValidation.mlx differ diff --git a/docs_part4/NSSolverValidation_images/figure_0.png b/docs_part4/NSSolverValidation_images/figure_0.png new file mode 100644 index 0000000..765472a Binary files /dev/null and b/docs_part4/NSSolverValidation_images/figure_0.png differ diff --git a/docs_part4/NSSolverValidation_images/figure_1.png b/docs_part4/NSSolverValidation_images/figure_1.png new file mode 100644 index 0000000..1d0cbdd Binary files /dev/null and b/docs_part4/NSSolverValidation_images/figure_1.png differ diff --git a/docs_part4/NSSolverValidation_images/image_0.png b/docs_part4/NSSolverValidation_images/image_0.png new file mode 100644 index 0000000..3b43b07 Binary files /dev/null and b/docs_part4/NSSolverValidation_images/image_0.png differ diff --git a/docs_part4/NSSolverValidation_images/image_1.png b/docs_part4/NSSolverValidation_images/image_1.png new file mode 100644 index 0000000..1b74d51 Binary files /dev/null and b/docs_part4/NSSolverValidation_images/image_1.png differ diff --git a/docs_part4/NSSolverValidation_images/image_2.png b/docs_part4/NSSolverValidation_images/image_2.png new file mode 100644 index 0000000..34cae4b Binary files /dev/null and b/docs_part4/NSSolverValidation_images/image_2.png differ diff --git a/docs_part4/checkNSSolverError.m b/docs_part4/checkNSSolverError.m new file mode 100644 index 0000000..c8cd573 --- /dev/null +++ b/docs_part4/checkNSSolverError.m @@ -0,0 +1,222 @@ +function [L2error, h_fig] = checkNSSolverError(Re,a,N,dt,tEnd,fuFunc,fvFunc,usolFunc,vsolFunc,psolFunc,timeScheme, visFlag, recordGIF) + +% Animation Setting +% recordGIF = false; % set to true if you need GIF +% visFlag = false; % set to true to visualize the flow +visRate = 4; % downsample rate of the data for quiver +recordRate = 20; +filename = 'errorCheck.gif'; % Specify the output file name + +% Simulation Setting +Lx = 1; Ly = 1; % domain size +Nx = N; Ny = N; % Number of grids + +% Grid size (Equispaced) +dx = Lx/Nx; +dy = Ly/Ny; +% Coordinate of each grid (cell center) +xce = ((1:Nx)-0.5)*dx; +yce = ((1:Ny)-0.5)*dy; +% Coordinate of each grid (cell corner) +xco = (0:Nx)*dx; +yco = (0:Ny)*dy; + +% cell center mesh grid +[Xce,Yce] = meshgrid(xce,yce); +%% Initialize the memory for flow field. +u = nan(Nx+1,Ny+2); % velocity in x direction (u) +v = nan(Nx+2,Ny+1); % velocity in y direction (v) +p = nan(Nx+1,Ny+1); % velocity in y direction (v) + +% For the current boundary condition +velbc.uTop = nan(size(u,1),1); +velbc.uBottom = nan(size(u,1),1); +velbc.uLeft = nan(1,size(u,2)); +velbc.uRight = nan(1,size(u,2)); +velbc.vTop = nan(size(v,1),1); +velbc.vBottom = nan(size(v,1),1); +velbc.vLeft = nan(1,size(v,2)); +velbc.vRight = nan(1,size(v,2)); + +% For the boundary conditoin at one time step ahead +velbc1 = velbc; + +%% Initialization of the flow field +t = 0; +[uX,uY] = meshgrid(xco(2:end-1),yce); +usol = usolFunc(uX', uY', t, Re, a); + +[vX,vY] = meshgrid(xce,yco(2:end-1)); +vsol = vsolFunc(vX', vY', t, Re, a); + +[pX,pY] = meshgrid(xce,yce); +p = psolFunc(pX', pY', t, Re, a); +u(2:end-1,2:end-1) = usol; +v(2:end-1,2:end-1) = vsol; + +% to be consistent with the velocity arrays, u, v +% Top, Buttom: column vector +% Left, Right: row vector +velbc = getVelocityBoudanryCondition(velbc, usolFunc, vsolFunc, xco, yco, xce, yce, t, Re, a); + +u(:,1) = 2*velbc1.uBottom - u(:,2); v(:,1) = velbc1.vBottom; %bottom +u(:,end) = 2*velbc1.uTop - u(:,end-1); v(:,end) = velbc1.vTop; %top +u(1,:) = velbc1.uLeft; v(1,:) = 2*velbc1.vLeft - v(2,:); %left +u(end,:) = velbc1.uRight; v(end,:) = 2*velbc1.vRight - v(end-1,:); %right + +% divergence check +inflow = sum(velbc.vBottom(2:end-1))*dx + sum(velbc.uLeft(2:end-1))*dy; +outflow = sum(velbc.vTop(2:end-1))*dx + sum(velbc.uRight(2:end-1))*dy; +assert(abs(inflow - outflow) < eps, "Inflow flux must match the outflow flux.") + + +if visFlag + %% Contour plot + h_fig = figure; + h_fig.Position = [350 150 700 300]; + ha1 = subplot(1,2,1); + ha2 = subplot(1,2,2); + + % Simulation Results + % get velocity at the cell center (for visualization) + uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center + vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center + [~,h_abs] = contourf(ha1, Xce',Yce',sqrt(uce.^2+vce.^2)); % contour + + % Analytical Results + % get velocity at the cell center (for visualization) + usol = usolFunc(Xce', Yce', t, Re, a); + vsol = vsolFunc(Xce', Yce', t, Re, a); + [~,h_absSol] = contourf(ha2, Xce',Yce',sqrt(usol.^2+vsol.^2)); % contour + + hold(ha1,'on'); + hold(ha2,'on'); + + % Some cosmetics + h_abs.LevelList = linspace(-1,1,24); + h_absSol.LevelList = linspace(-1,1,24); + ha1.CLim = [0,1]; + ha2.CLim = [0,1]; + ha1.FontSize = 14; + ha2.FontSize = 14; + ha1.XTick = [];ha1.YTick = []; + ha2.XTick = [];ha2.YTick = []; + + %% Quiver plot + % Downsample the dataid = downsampledj + xced = xce(1:visRate:end); + yced = yce(1:visRate:end); + [Xced,Yced] = meshgrid(xced, yced); + + % Simulation Results + uced = uce(1:visRate:end,1:visRate:end); + vced = vce(1:visRate:end,1:visRate:end); + h_quiver = quiver(ha1, Xced',Yced',uced,vced,1,'Color',[0,0,0]); + + % Analytical Results + usold = usolFunc(Xced', Yced', t, Re, a); + vsold = vsolFunc(Xced', Yced', t, Re, a); + h_quiverSol = quiver(ha2, Xced',Yced',usold,vsold,1,'Color',[0,0,0]); + + hold(ha1,'off'); + xlim(ha1,[0 Lx]); ylim(ha1,[0 Ly]); + hold(ha2,'off'); + xlim(ha2,[0 Lx]); ylim(ha2,[0 Ly]); + + + title(ha1,'Simulation Results'); + title(ha2,'Analytical Results'); + axis(ha1,'equal'); + axis(ha2,'equal'); +end + +%% Start the Simulation +% Just to make it a little fun, the velocity of the top lid reverses at the +% 1000 steps (out of nt steps). + +initialFrame = true; +index = 1; +while t < tEnd + + velbc = getVelocityBoudanryCondition(velbc1, usolFunc, vsolFunc, xco, yco, xce, yce, t, Re, a); + velbc1 = getVelocityBoudanryCondition(velbc1, usolFunc, vsolFunc, xco, yco, xce, yce, t+dt, Re, a); + + uForce = fuFunc(uX', uY', t, Re, a); + vForce = fvFunc(vX', vY', t, Re, a); + % psol = psolFunc(Xce', Yce', t, Re, a); + + % Update the velocity field (uses dct) + switch timeScheme + case 'Euler' + [u,v,p] = updateVelocityField_Euler(u,v,[],Nx,Ny,dx,dy,Re,dt,velbc,'dct',uForce,vForce,velbc1); + case 'CNAB' + [u,v,p] = updateVelocityField_CNAB( u,v,Nx,Ny,dx,dy,Re,dt,velbc,'dct',uForce,vForce,velbc1); + case 'RK3' + [u,v,p] = updateVelocityField_RK3(u,v,Nx,Ny,dx,dy,Re,dt,velbc,'dct',uForce,vForce,velbc1); + otherwise + error('timeScheme is not recognized'); + end + + t = t + dt; + + % get velocity at the cell center (for visualization) + uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center + vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center + usol = usolFunc(Xce', Yce', t, Re, a); + vsol = vsolFunc(Xce', Yce', t, Re, a); + + uSim = [uce(:); vce(:)]; + uSol = [usol(:); vsol(:)]; + L2error = norm(uSim-uSol)/norm(uSol); + + CFL = max([u(:)/dx; v(:)/dy])*dt; +% disp("error: " + num2str(L2error) + " at " + num2str(t) + " CFL: " + num2str(CFL)); + + % Update the plot at every recordRate steps + if visFlag && mod(index,recordRate) == 0 + + % update plot (downsample) + h_quiver.UData = uce(1:visRate:end,1:visRate:end); + h_quiver.VData = vce(1:visRate:end,1:visRate:end); + h_abs.ZData = sqrt(uce.^2+vce.^2); + + h_quiverSol.UData = usol(1:visRate:end,1:visRate:end); + h_quiverSol.VData = vsol(1:visRate:end,1:visRate:end); + h_absSol.ZData = sqrt(usol.^2+vsol.^2); + + drawnow + + if recordGIF + frame = getframe(gcf); %#ok % Figure ʂ[r[t[i\́jƂăLv` + tmp = frame2im(frame); % 摜ɕύX + [A,map] = rgb2ind(tmp,256); % RGB -> CfbNX摜 + if initialFrame + imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.1); + initialFrame = false; + else + imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);% 摜Ayh + end + end + end + + index = index + 1; + + +end + disp("(N,dt) = (" + num2str(N) + "," + num2str(dt) +") error: " + num2str(L2error) + " at " + num2str(t) + " CFL: " + num2str(CFL)); + +end + +function velbc = getVelocityBoudanryCondition(velbc, usolFunc, vsolFunc, xco, yco, xce, yce, t, Re, a) + + +velbc.uTop = usolFunc(xco, 1, t, Re, a)'; +velbc.uBottom = usolFunc(xco, 0, t, Re, a)'; +velbc.uLeft(2:end-1) = usolFunc(0, yce, t, Re, a); +velbc.uRight(2:end-1) = usolFunc(1, yce, t, Re, a); +velbc.vTop(2:end-1) = vsolFunc(xce, 1, t, Re, a)'; +velbc.vBottom(2:end-1) = vsolFunc(xce, 0, t, Re, a)'; +velbc.vLeft = vsolFunc(0, yco, t, Re, a); +velbc.vRight = vsolFunc(1, yco, t, Re, a); + +end \ No newline at end of file diff --git a/docs_part4/errorCheck.gif b/docs_part4/errorCheck.gif new file mode 100644 index 0000000..b2de0d8 Binary files /dev/null and b/docs_part4/errorCheck.gif differ diff --git a/functions/updateVelocityField_CNAB.m b/functions/updateVelocityField_CNAB.m index c3cdd28..a7ca490 100644 --- a/functions/updateVelocityField_CNAB.m +++ b/functions/updateVelocityField_CNAB.m @@ -105,22 +105,21 @@ error("Specified method: " + method + " is not supported." + ... "It should be either direct or dct"); end -% correction to get the final velocity -p = dp; switch method case 'dct_p' - u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; - v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; - u(end,2:end-1) = u(end,2:end-1) + 2*p(end,:)/dx; %right - v(2:end-1,end) = v(2:end-1,end) + 2*p(:,end)/dy; %right - v(2:end-1,1) = v(2:end-1,1) - 2*p(:,1)/dy; %right + u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (dp(2:end,:)-dp(1:end-1,:))/dx; + v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (dp(:,2:end)-dp(:,1:end-1))/dy; + u(end,2:end-1) = u(end,2:end-1) + 2*dp(end,:)/dx; %right + v(2:end-1,end) = v(2:end-1,end) + 2*dp(:,end)/dy; %top + v(2:end-1,1) = v(2:end-1,1) - 2*dp(:,1)/dy; %bottom otherwise - u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; - v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; + u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (dp(2:end,:)-dp(1:end-1,:))/dx; + v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (dp(:,2:end)-dp(:,1:end-1))/dy; end - +% correction to get the final velocity +p = dp; % check the divergence % b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ... diff --git a/functions/updateVelocityField_RK3.m b/functions/updateVelocityField_RK3.m index 6b96d19..339e63a 100644 --- a/functions/updateVelocityField_RK3.m +++ b/functions/updateVelocityField_RK3.m @@ -116,28 +116,29 @@ % by using the discrete cosine transformiRTCϊgpj % Note: Signal Processing Toolbox required dp = solvePoissonEquation_2dDCT(b,nx,ny,dx,dy); - case 'dct_p' + case 'dct_p' % by using the discrete cosine transformiRTCϊgpj % Note: Signal Processing Toolbox required dp = solvePoissonEquation_2dDCT_p(b,nx,ny,dx,dy); - + otherwise error("Specified method: " + method + " is not supported." + ... "It should be either direct or dct"); end -% correction to get the final velocity -p = dp; -% u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; -% v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; +% correction to get the final velocity switch method case 'dct_p' - u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; - v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; - u(end,2:end-1) = u(end,2:end-1) + 2*p(end,:)/dx; %right + u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (dp(2:end,:)-dp(1:end-1,:))/dx; + v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (dp(:,2:end)-dp(:,1:end-1))/dy; + u(end,2:end-1) = u(end,2:end-1) + 2*dp(end,:)/dx; %right + v(2:end-1,end) = v(2:end-1,end) + 2*dp(:,end)/dy; %top + v(2:end-1,1) = v(2:end-1,1) - 2*dp(:,1)/dy; %bottom otherwise - u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (p(2:end,:)-p(1:end-1,:))/dx; - v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (p(:,2:end)-p(:,1:end-1))/dy; + u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - (dp(2:end,:)-dp(1:end-1,:))/dx; + v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - (dp(:,2:end)-dp(:,1:end-1))/dy; end +p = dp; + end \ No newline at end of file