Skip to content

Commit

Permalink
Restructored folder setup
Browse files Browse the repository at this point in the history
  • Loading branch information
minoue-xx committed Mar 16, 2020
1 parent 0fdaab4 commit e3aef1e
Show file tree
Hide file tree
Showing 22 changed files with 26 additions and 17 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
%

%% Setting up

clear, close all
Re = 500; % Reynolds number
nt = 2000; % max time steps
Lx = 1; Ly = 1; % domain size
Expand Down Expand Up @@ -36,7 +36,7 @@
% |recordRate| time steps.

visRate = 4; % downsample rate of the data for quiver
recordGIF = true; % set to true if you wanna make GIF
recordGIF = false; % set to true if you wanna make GIF
recordRate = 20;
filename = 'animation_sample.gif'; % Specify the output file name
%%
Expand Down Expand Up @@ -85,7 +85,7 @@
end

% Update the velocity field (uses dct)
[u,v] = updateVelocityField(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
[u,v] = updateVelocityField_Euler(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');

% Update the plot at every recordRate steps
if mod(ii,recordRate) == 0
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -338,7 +338,9 @@ Enough equations already. Let's put everything in MATLAB.


```matlab
clear, close all
clear
close all
addpath('../functions/');
```
## Setting up Computational Domain

Expand Down Expand Up @@ -554,7 +556,7 @@ Good.
# Animation of the Velocity Field


Making animation of the flow field is the fun part of CFD (Colorful Fluid Dynamics). In addition to the contour plot, let's create an arrow plot that represents the velocity using `quiver` function. The code below generates GIF. The numerical integration process discussed above is now in the function `updateVelocityField.m.`
Making animation of the flow field is the fun part of CFD (Colorful Fluid Dynamics). In addition to the contour plot, let's create an arrow plot that represents the velocity using `quiver` function. The code below generates GIF. The numerical integration process discussed above is now in the function `updateVelocityField_Euler.m.`



Expand Down Expand Up @@ -672,7 +674,7 @@ for ii = 1:2000
end
% Update the velocity field (uses dct)
[u,v] = updateVelocityField(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
[u,v] = updateVelocityField_Euler(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% Update the plot at every recordRate steps
if mod(ii,recordRate) == 0
Expand Down
Binary file not shown.
Original file line number Diff line number Diff line change
Expand Up @@ -419,7 +419,9 @@ p(1:end, 1:end)


```matlab
clear, close all
clear
close all
addpath('../functions/');
```
## 解析領域の設定

Expand Down Expand Up @@ -638,7 +640,7 @@ b = ((u(2:end,2:end-1)-u(1:end-1,2:end-1))/dx ...
disp(norm(b))
```
```
4.7044e-14
5.8843e-14
```


Expand All @@ -648,7 +650,7 @@ disp(norm(b))
# 流れ場のアニメーション表示


CFD は Colorful Fluid Dynamics の略とも言われるだけあって、可視化が数値流体の醍醐味ですよね。等高線図だけじゃなくて速度を矢印で表示させて GIF 出力します。上の計算処理は `updateVelocityField.m` として関数化しておきます。
CFD は Colorful Fluid Dynamics の略とも言われるだけあって、可視化が数値流体の醍醐味ですよね。等高線図だけじゃなくて速度を矢印で表示させて GIF 出力します。上の計算処理は `updateVelocityField_Euler.m` として関数化しておきます。



Expand Down Expand Up @@ -712,7 +714,7 @@ figure
[~,h_abs] = contourf(Xce',Yce',sqrt(uce.^2+vce.^2)); % 等高線図
```
```
Warning: Contour not rendered for constant ZData
警告: 等高線図は ZData が定数の場合はレンダリングされません
```
```matlab
hold on
Expand Down Expand Up @@ -771,7 +773,7 @@ for ii = 1:2000
end
% 速度場更新(コサイン変換使用)
[u,v] = updateVelocityField(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
[u,v] = updateVelocityField_Euler(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% 描画は recordRate 毎に実施
if mod(ii,recordRate) == 0
Expand Down
Binary file not shown.
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
phat(1,1) = 0;

% Inverse 2D DCT
p = idct2(phat); % Needs Image Processing Toolbox
% u = idct(idct(uhat)')'; % Same as above (Needs Signal Processing Toolbox instead)
% p = idct2(phat); % Needs Image Processing Toolbox
p = idct(idct(phat)')'; % Same as above (Needs Signal Processing Toolbox instead)

end
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
function dp = solvePoissonEquation_direct(b,nx,ny,dx,dy)
% Copyright 2020 The MathWorks, Inc.

persistent Abig
persistent Abig sizeAbig

if isempty(Abig) || any(size(Abig) ~= [nx*ny,nx*ny])
% So that you do not have to calculate Abig everytime
if isempty(Abig) || any(sizeAbig ~= [nx*ny,nx*ny])

% 行列 A の構築(メモリ節約のためスパース行列で定義)
% まず x 方向の微分に関する部分(三重対角行列)から定義します。
Expand Down Expand Up @@ -32,6 +34,9 @@
Abig(1,:) = 0;
Abig(1,1) = 1;

sizeAbig = size(Abig);
% Pre-decompose the matrix for fater inversion
Abig = decomposition(Abig);
end

% 右辺
Expand Down
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [u,v] = updateVelocityField(u,v,nx,ny,dx,dy,Re,dt,bctop,method)
function [u,v] = updateVelocityField_Euler(u,v,nx,ny,dx,dy,Re,dt,bctop,method)
% Copyright (c) 2020, The MathWorks, Inc.

% Apply boundary conditions:
Expand Down Expand Up @@ -59,7 +59,7 @@
% by using the discrete cosine transform(コサイン変換使用)
% Note: Signal Processing Toolbox required
dp = solvePoissonEquation_2dDCT(b,nx,ny,dx,dy);
case 'minres'
case 'iterative'
[dp,~] = minres(@(x) operatorDG(x,nx,ny,dx,dy),b(:),1e-4,300);
dp = reshape(dp,nx,ny);

Expand Down

0 comments on commit e3aef1e

Please sign in to comment.