Update README.md

README.md

@@ -4,23 +4,34 @@
 This repo provides a MATLAB example code for the lid-driven cavity flow where incompressible 
 Navier Stokes equation is numerically solved using a simple 2nd order finite difference scheme on a staggered grid system.
 
+![sample](./gif/animation_sample1e2and1e4.gif)
+(Left: Re = 100, Right: Re = 10,000)
 
-![sample](./gif/animation_sample.gif)
-The top lid changes its direction. The arrow denotes the velocity field and the contour denotes its magnitude.
+The top lid changes its direction. The arrow denotes the velocity field, and the contour denotes its magnitude.
 
 ## Part 1
 
+
 - Click [here](./docs_part1/vanilaCavityFlow_EN.md) for detailed documentation in English.
 - 日本語のドキュメントは[こちら](./docs_part1/vanilaCavityFlow_JP.md) から
 
-In order to focus on understanding basic idea of the numerical schemes, the method is kept premitive; explicit treatment of viscous term (the solution will diverge at low Reynolds number) and
-the time integration is Euler.
+The numerical scheme is kept primitive; the explicit treatment of viscous term (the solution diverges at low Reynolds number), and the time integration is Euler.
 
 まずは単純な手法でキャビティ流れのシミュレーションを実施します。
 
+## Part 2
+
+- Click [here](./docs_part2/vanilaCavityFlowImplicit_EN.md) for detailed documentation in English.
+- 日本語のドキュメントは[こちら](./docs_part2/vanilaCavityFlowImplicit_JP.md) から
+
+The implicit treatments for viscous terms are implemented at low Reynolds number, namely the Crank-Nicolson method. For better stability for non-linear terms, Adams-Bashforth, and 3 steps-Runge-Kutta is also implemented. 
+
+拡散項に対して陰解法を実装しました。対流項へアダムス・バッシュフォースを使用したもの、3段階のルンゲクッタ法の２つの時間発展を実装しています。
+
+
 ## Next to come
 
-Plan is to implement the implicit treatment for viscous terms to obtain better stability.
+The plan is to allow arbitrary boundary conditions for more fun simulations.
 
 
 
@@ -33,7 +44,7 @@ Plan is to implement the implicit treatment for viscous terms to obtain better s
 # ToDo
 
 1. Implement implicit treatment of viscous terms
-2. Imlement crank-Nicolson for the non-linear terms
+2. Implement crank-Nicolson for the non-linear terms
 3. Allow obstacles within the domain
 4. Allow inflow from the wall
 5. Make it to 3D

functions/getIntermediateU_xyCNAB.m

@@ -0,0 +1,24 @@
+% Copyright 2020 The MathWorks, Inc.
+function xu = getIntermediateU_xyCNAB(u, b, dt, Re, nx, ny, dx, dy)
+
+persistent A4u dtRe_old L4u
+
+if isempty(dtRe_old)
+    dtRe_old = dt/Re;
+end
+
+if isempty(A4u) || any(size(L4u) ~= [(nx-1)*ny,(nx-1)*ny]) || dtRe_old ~= dt/Re
+    dtRe_old = dt/Re;
+    maskU = false(nx+1,ny+2);
+    maskU(2:end-1,2:end-1) = true;
+    L4u = getL4u(nx,ny,dx,dy,maskU);
+    A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u;
+    A4u = decomposition(A4u); 
+end
+
+% x0 = u(2:end-1,2:end-1);
+% [xu,flag] = cgs(A4u,b(:),[],[],[],[],x0(:));
+xu = A4u\b(:);
+xu = reshape(xu,[nx-1,ny]);
+
+end

functions/getIntermediateU_xyRK3.m

@@ -0,0 +1,37 @@
+% Copyright 2020 The MathWorks, Inc.
+function xu = getIntermediateU_xyRK3(u, b, dt, Re, nx, ny, dx, dy, id)
+
+persistent A4u1 A4u2 A4u3 dtRe_old L4u
+
+if isempty(dtRe_old)
+    dtRe_old = dt/Re;
+end
+
+if isempty(A4u1) || any(size(L4u) ~= [(nx-1)*ny,(nx-1)*ny]) || dtRe_old ~= dt/Re
+    dtRe_old = dt/Re;
+    maskU = false(nx+1,ny+2);
+    maskU(2:end-1,2:end-1) = true;
+    L4u = getL4u(nx,ny,dx,dy,maskU);
+    A4u1 = eye(size(L4u),'like',L4u)-(8/15)*dt/(Re)*L4u;
+    A4u2 = eye(size(L4u),'like',L4u)-(2/15)*dt/(Re)*L4u;
+    A4u3 = eye(size(L4u),'like',L4u)-(2/6)*dt/(Re)*L4u;
+    A4u1 = decomposition(A4u1);
+    A4u2 = decomposition(A4u2);
+    A4u3 = decomposition(A4u3);
+
+end
+
+% x0 = u(2:end-1,2:end-1);
+% [xu,flag] = cgs(A4u,b(:),[],[],[],[],x0(:));
+switch id
+    case 1
+        xu = A4u1\b(:);
+    case 2
+        xu = A4u2\b(:);
+    case 3
+        xu = A4u3\b(:);
+end
+
+xu = reshape(xu,[nx-1,ny]);
+
+end

functions/getIntermediateV_xyCNAB.m

@@ -0,0 +1,25 @@
+% Copyright 2020 The MathWorks, Inc.
+function xv = getIntermediateV_xyCNAB(v, b, dt, Re, nx, ny, dx, dy)
+
+persistent A4v dtRe_old L4v
+
+if isempty(dtRe_old)
+    dtRe_old = dt/Re;
+end
+
+if isempty(A4v) || any(size(L4v) ~= [nx*(ny-1),nx*(ny-1)]) || dtRe_old ~= dt/Re
+    dtRe_old = dt/Re;
+    maskV = false(nx+2,ny+1);
+    maskV(2:end-1,2:end-1) = true;
+    L4v = getL4v(nx,ny,dx,dy,maskV);
+    A4v = eye(size(L4v),'like',L4v)-dt/(2*Re)*L4v;
+    A4v = decomposition(A4v);  
+end
+
+% x0 = v(2:end-1,2:end-1);
+% [xv,flag] = cgs(A4v,b(:),[],[],[],[],x0(:));
+xv = A4v\b(:);
+xv = reshape(xv,[nx,ny-1]);
+
+end
+

functions/getIntermediateV_xyRK3.m

@@ -0,0 +1,35 @@
+% Copyright 2020 The MathWorks, Inc.
+function xv = getIntermediateV_xyRK3(v, b, dt, Re, nx, ny, dx, dy,id)
+
+persistent A4v1 A4v2 A4v3 dtRe_old L4v
+
+if isempty(dtRe_old)
+    dtRe_old = dt/Re;
+end
+
+if isempty(A4v1) || any(size(L4v) ~= [nx*(ny-1),nx*(ny-1)]) || dtRe_old ~= dt/Re
+    dtRe_old = dt/Re;
+    maskV = false(nx+2,ny+1);
+    maskV(2:end-1,2:end-1) = true;
+    L4v = getL4v(nx,ny,dx,dy,maskV);
+    A4v1 = eye(size(L4v),'like',L4v)-(8/15)*dt/(Re)*L4v;
+    A4v2 = eye(size(L4v),'like',L4v)-(2/15)*dt/(Re)*L4v;
+    A4v3 = eye(size(L4v),'like',L4v)-(2/6)*dt/(Re)*L4v;
+    A4v1 = decomposition(A4v1);
+    A4v2 = decomposition(A4v2);
+    A4v3 = decomposition(A4v3);
+
+end
+
+switch id
+    case 1
+        xv = A4v1\b(:);
+    case 2
+        xv = A4v2\b(:);
+    case 3
+        xv = A4v3\b(:);
+end
+xv = reshape(xv,[nx,ny-1]);
+
+end
+

functions/getL4u.m

@@ -0,0 +1,64 @@
+% Copyright 2020 The MathWorks, Inc.
+function L4u = getL4u(nx,ny,dx,dy,mask)
+
+% Note: Boundary condition effects y derivative only
+% y-derivatives
+% ub = (u0+u1)/2 => u0 = 2*ub - u1;
+% Lu_1 = (u2 - 2*u1 + u0)/dy^2
+%      = (u2 - 2*u1 + 2*ub - u1)/dy^2
+%      = (u2 - 3*u1)/dx^2 + 2*ub/dy^2
+% 
+% x-derivatives
+% Lu_1 = (u2 - 2*u1 + ub)/dx^2
+%      = (u2 - 2*u1)/dx^2 + ub/dx^2
+
+matrixSize = [nx-1,ny];
+
+% 3
+%152
+% 4
+coefx = 1/dx^2;
+coefy = 1/dy^2;
+
+i = zeros(nx*ny*5,1);
+j = zeros(nx*ny*5,1);
+v = zeros(nx*ny*5,1);
+index1 = 1;
+index = 1;
+for jj=1:ny
+    for ii=1:nx-1
+        idx = [mask(ii+2,jj+1),mask(ii,jj+1),mask(ii+1,jj+2),mask(ii+1,jj),mask(ii+1,jj+1)];
+        stencils = [ii+1,jj
+            ii-1,jj
+            ii,jj+1
+            ii,jj-1
+            ii,jj];
+        stencils = stencils(idx,:);
+
+        tmpx = 2;
+        tmpy = 2 + sum(~idx(3:4));
+
+        coeffs = [coefx
+            coefx
+            coefy
+            coefy
+            -tmpx*coefx-tmpy*coefy];
+        coeffs = coeffs(idx);
+        linearIdx = sub2ind(matrixSize, stencils(:,1), stencils(:,2));
+
+        i(index1:index1+length(linearIdx)-1) = index;
+        j(index1:index1+length(linearIdx)-1) = linearIdx;
+        v(index1:index1+length(linearIdx)-1) = coeffs;
+        index1 = index1 + length(linearIdx);
+        index = index + 1;
+    end
+end
+
+i(i==0) = [];
+j(j==0) = [];
+v(v==0) = [];
+L4u = sparse(i,j,v,(nx-1)*ny,(nx-1)*ny);
+L4u(:,~any(L4u,1)) = [];
+L4u(~any(L4u,2),:) = [];
+
+end

functions/getL4v.m

@@ -0,0 +1,64 @@
+% Copyright 2020 The MathWorks, Inc.
+function L4v = getL4v(nx,ny,dx,dy,mask)
+
+% Note: Boundary condition effects x derivative only
+% x-derivatives
+% vb = (v0+v1)/2 => v0 = 2*vb - v1;
+% Lv_1 = (v2 - 2*v1 + v0)/dx^2
+%      = (v2 - 2*v1 + 2*vb - u1)/dx^2
+%      = (v2 - 3*v1)/dx^2 + 2*vb/dx^2
+% 
+% y-derivatives
+% Lv_1 = (v2 - 2*v1 + vb)/dy^2
+%      = (v2 - 2*v1)/dy^2 + vb/dy^2
+
+matrixSize = [nx,ny-1];
+
+% 3
+%152
+% 4
+coefx = 1/dx^2;
+coefy = 1/dy^2;
+
+i = zeros(nx*ny*5,1);
+j = zeros(nx*ny*5,1);
+v = zeros(nx*ny*5,1);
+index1 = 1;
+index = 1;
+for jj=1:ny-1
+    for ii=1:nx
+        idx = [mask(ii+2,jj+1),mask(ii,jj+1),mask(ii+1,jj+2),mask(ii+1,jj),mask(ii+1,jj+1)];
+        stencils = [ii+1,jj
+            ii-1,jj
+            ii,jj+1
+            ii,jj-1
+            ii,jj];
+        stencils = stencils(idx,:);
+
+        tmpx = 2 + sum(~idx(1:2));
+        tmpy = 2;
+
+        coeffs = [coefx
+            coefx
+            coefy
+            coefy
+            -tmpx*coefx-tmpy*coefy];
+        coeffs = coeffs(idx);
+        linearIdx = sub2ind(matrixSize, stencils(:,1), stencils(:,2));
+
+        i(index1:index1+length(linearIdx)-1) = index;
+        j(index1:index1+length(linearIdx)-1) = linearIdx;
+        v(index1:index1+length(linearIdx)-1) = coeffs;
+        index1 = index1 + length(linearIdx);
+        index = index + 1;
+    end
+end
+
+i(i==0) = [];
+j(j==0) = [];
+v(v==0) = [];
+L4v = sparse(i,j,v,nx*(ny-1),nx*(ny-1));
+L4v(:,~any(L4v,1)) = [];
+L4v(~any(L4v,2),:) = [];
+
+end

functions/operatorAu_CNAB.m

@@ -0,0 +1,19 @@
+% Copyright 2020 The MathWorks, Inc.
+function Ax = operatorAu_CNAB(x,dt,Re,nx,ny,dx,dy)
+% Getting A for u velocity (Crank-Nicolson)
+% A = [I - dt/(2*Re)*Lxy];
+
+x = reshape(x,nx-1,ny);
+xbig = zeros(nx+1,ny+2);
+xbig(2:end-1,2:end-1) = x;
+
+% Addting boudanry value
+% ub = (u0+u1)/2 => u0 = 2*ub - u1;
+xbig(:,1) = - xbig(:,2);
+xbig(:,end) = - xbig(:,end-1);
+
+Ax = x - dt/(2*Re)*(xbig(1:end-2,2:end-1)-2*xbig(2:end-1,2:end-1)+xbig(3:end,2:end-1))/dx^2; % nx-1 * ny
+Ax = Ax - dt/(2*Re)*(xbig(2:end-1,1:end-2)-2*xbig(2:end-1,2:end-1)+xbig(2:end-1,3:end))/dy^2; % nx-1 * ny
+Ax = Ax(:);
+
+end

functions/operatorAv_CNAB.m

@@ -0,0 +1,19 @@
+% Copyright 2020 The MathWorks, Inc.
+function Ax = operatorAv_CNAB(x,dt,Re,nx,ny,dx,dy)
+% Getting A for v velocity (Crank-Nicolson)
+% A = [I - dt/(2*Re)*Lxy];
+
+x = reshape(x,nx,ny-1);
+xbig = zeros(nx+2,ny+1);
+xbig(2:end-1,2:end-1) = x;
+
+% Addting boudanry value
+% ub = (u0+u1)/2 => u0 = 2*ub - u1;
+xbig(1,:) = - xbig(2,:);
+xbig(end,:) = - xbig(end-1,:);
+
+Ax = x - dt/(2*Re)*(xbig(1:end-2,2:end-1)-2*xbig(2:end-1,2:end-1)+xbig(3:end,2:end-1))/dx^2; % nx-1 * ny
+Ax = Ax - dt/(2*Re)*(xbig(2:end-1,1:end-2)-2*xbig(2:end-1,2:end-1)+xbig(2:end-1,3:end))/dy^2; % nx-1 * ny
+Ax = Ax(:);
+end
+