draw_grid.m
代码语言：javascript复制%DRAW_GRID
% Screen plot of grid
tic
[X,Y] = meshgrid([0,cumsum(dx)],[0,cumsum(dy)]);
figure(1), clf, hold on
title('Grid','FontSize',16)
plot(X,Y,'k') % Vertical grid lines
axis('equal')
plot(X',Y','k') % Horizontal grid lines
tijd = toc; disp(['draw_grid time = ',num2str(tijd)])
errornorm.m
代码语言：javascript复制%ERRORNORM
% Maximum norm of error in outflow profile if exact solution is available
% Function called: exact_solution
if geval == 1 % Horizontal Poiseuille flow
uq = reshape(u1,size(XU')); uq = uq';
disp(['outflow_errormaxnorm = ',...
num2str(norm(uq(:,J 1) - exact_solution(0, 0, yu', 'right'),inf))])
projexsol = (exact_solution(0, 0, yu'- dy'/2, 'right') ...
4*exact_solution(0, 0, yu', 'right') ...
exact_solution(0, 0, yu' dy'/2, 'right'))/6;
disp(['projected_outflow_errormaxnorm = ',...
num2str(norm(uq(:,J 1) - projexsol,inf))])
elseif geval == 2 % Vertical Poiseuille flow
vq = reshape(v1,size(XV')); vq = vq';
disp(['outflow_errormaxnorm = ',...
num2str(norm(vq(K 1,:) - exact_solution(0, xv, 0, 'upper'),inf))])
projexsol = (exact_solution(0, xv - dx/2, 0, 'upper') ...
4*exact_solution(0, xv, 0, 'upper') ...
exact_solution(0, xv dx/2, 0, 'upper'))/6;
disp(['projected_outflow_errormaxnorm = ',...
num2str(norm(vq(K 1,:) - projexsol,inf))])
else
end
exact_solution.m
代码语言：javascript复制function ye = exact_solution(t, x, y, side)
%EXACT_SOLUTION
% Possible values for side: 'lower', 'upper', 'left', 'right'
global geval xu yu xv yv
ye = 0;
if geval == 1 % Horizontal Poiseuille flow
if strcmp(side, 'right') == 1
ye = 6*(y - yv(1)).*(yv(end) - y)/(yv(end) - yv(1))^3;
else
ye = 0;
end
elseif geval == 2 % Vertical Poiseuille flow
if strcmp(side, 'upper') == 1
ye = 6*(x - xu(1)).*(xu(end) - x)/(xu(end) - xu(1))^3;
else
ye = 0;
end
elseif geval == 7 % Horizontal Poiseuille flow
if strcmp(side, 'right') == 1
ye = -6*(y - yv(1)).*(yv(end) - y)/(yv(end) - yv(1))^3;
else
ye = 0;
end
else
error('No exact solution for this case')
end
grad_and_div.m
代码语言：javascript复制%GRAD_AND_DIV
% Generates pressure gradient matrices Pu, Pv and divergence matrix D
% Vectorized matrix generation by evaluation of stencil coefficients
% Output: Pu, Pv, Du, Dv
J= length(dx); K = length(dy);
%.................Pressure gradient matrix Pu ....................
% | 0 |
% Stencil: [Pu] = |p1 p2 0 |
% | 0 |
% Indexing convention in staggered grid:
% --k 1--
% | |
% j jk j 1
% | |
% ---k---
tic
p2 = 1./DXU(1,:)'; p1 = zeros(size(p2)); p1(1:J) = -p2(2:J 1);
Puu = spdiags([p1 p2],[-1;0], J 1, J); Pu = kron(speye(K),Puu);
%......................Boundary corrections....................................
for seg = 1:length(ybc(1,:))
if (ybc(1,seg) == 1)|(ybc(1,seg) == 2) % No-slip or inflow at left boundary
k = 1 kseg(seg):kseg(seg 1); Pu(1 (k-1)*(J 1),:) = 0;
end
if (ybc(2,seg) == 1)|(ybc(2,seg) == 2) % No-slip or inflow at right boundary
k = 1 kseg(seg):kseg(seg 1); Pu(k*(J 1),:) = 0;
end
end
tijd = toc; disp(['Breakdown of grad_and_div time'])
disp([' Pu time = ',num2str(tijd)])
%.................Pressure gradient matrix Pv ....................
% | 0 |
% Stencil: [Pv] = |0 p2 0|
% | p1 |
tic
pp1 = zeros(size(XV)); pp2 = pp1; % Diagonals of Pv
pp1(2:K 1,:) = -1./DYV(2:K 1,:); pp2 = -pp1;
pp2(1,:) = 1./DYV(1,:); pp2(K 1,:) = 0;
%......................Boundary corrections....................................
for seg = 1:length(xbc(1,:))
if (xbc(1,seg) == 1)|(xbc(1,seg) == 2) % No-slip or inflow at lower boundary
j = 1 jseg(seg):jseg(seg 1); pp1(1,j) = 0; pp2(1,j) = 0;
end
if (xbc(2,seg) == 1)|(xbc(2,seg) == 2) % No-slip or inflow at upper boundary
j = 1 jseg(seg):jseg(seg 1); pp1(K 1,j) = 0; pp2(K 1,j) = 0;
end
end
n = J*(K 1); p1 = reshape(pp1',n,1); p2 = reshape(pp2',n,1);
p1 = [p1(J 1:n); zeros(J,1)]; % Shift to accomodate spdiags
Pv = spdiags([p1 p2],[-J;0], n,n-J);
tijd = toc; disp([' Pv time = ',num2str(tijd)])
%.................u divergence matrix Du ....................
% | 0 |
% Stencil: [Du] = |0 p1 p2|
% | 0 |
tic
Duu = spdiags([-ones(J,1) ones(J,1)], [0;1], J,J 1);
Du = kron(spdiags(dy',0,K,K),Duu);
tijd = toc; disp([' Du time = ',num2str(tijd)])
%.................v divergence matrix Dv ....................
% | p2 |
% Stencil: [Dv] = |0 p1 0|
% | 0 |
tic
Dvv = spdiags([-ones(K,1) ones(K,1)], [0;1], K,K 1);
Dv = kron(Dvv,spdiags(dx',0,J,J));
tijd = toc; disp([' Dv time = ',num2str(tijd)])
clear pp1 pp2 p1 p2 Puu Duu Dvv % Save storage
