惯性矩阵中惯量项的离散化矩阵与右手边的修正。

%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
%  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

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
% 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