Streamline function keeps overwriting plot in a for loop
2 次查看(过去 30 天)
显示 更早的评论
I am trying to create a streamline plot of a 2D field that is evolving in time. I would like to visualize how the streamlines change over time and I am currently using a for loop to loop through each iteration in time. However, when I use the streamline function to plot it, I can't seem to clear the previous iteration of streamlines and instead they keep layering on top of the previous one. Using the close command on the figure does work, but this causes the figure to flash on the screen for a quick second before it closes again, and I am trying to get it to stay on the screen and keep updating like a normal plot would.
2 个评论
Shaan Stephen
2023-12-6
I am trying to do the lid driven cavity problem for some CFD code so it is pretty bulky. The section where I print the streamline function is in the fractial step portion
clc; close all; clear;
global np nu nx ny dx dy dt iu iv ip noo uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B
nx = 64 ;
ny = 64 ;
Re = 1000 ;
noo = 1/Re;
%BCs ...
uBC_L = zeros(1, ny);
uBC_R = zeros(1, ny);
uBC_B = 0;
uBC_T = 1;
vBC_L = 0;
vBC_R = 0;
vBC_B = 0;
vBC_T = 0;
% % w/ jet
% midpoint = ny/2;
% midjet = .25;
% uBC_L(midpoint) = midjet; uBC_L(midpoint+1) = midjet; uBC_L(midpoint-1) = midjet;
% Restart?
Restart = false ;
% pointer and grid
[ip, iu, iv] = GenPointer(nx, ny) ;
% time step
dt = .01;
dx = 1 / nx ;
dy = 1 / ny ;
np = nx * ny ;
nu = 2*nx*ny - nx - ny ;
nt = 60/dt;
uinitial = zeros(nu, 1);
u_BC_div = BC_Div(uBC_L, uBC_R, vBC_T, vBC_B);
u_BC_Lap = BC_Laplace();
pressure = cell([nt 1]);
% explicit euler
uf = uinitial + dt * (Adv(uinitial) + noo * Laplace(uinitial));
b_EE = (1/dt) * Div(uf) + (1/dt) * u_BC_div;
pressure{1} = CG_EE(b_EE);
u_new = uf - dt * Grad(cell2mat(pressure{1}));
u{1} = uinitial;
u{2} = u_new;
x = 0:dx:1;
y = 0:dy:1;
x = x(2:end-1);
y = y(2:end);
%%
P_array = nan(np, 1);
for i = 1:np
P_array(i) = i;
end
test1 = Grad(P_array);
h = 0;
xy = cell(nt, 1);
%% time stepping using fractial step
for i = 2:1:nt
% fractial step: stage 1
RHS_b = Rinverse(u{i}) + dt * (3 * Adv(u{i})- Adv(u{i-1})) * .5 + dt * noo * u_BC_Lap;
uf = CG_stage1(RHS_b) ;
% fractial step: stage 2
RHS_b = Div(cell2mat(uf)) /dt + BC_Div(uBC_L, uBC_R, vBC_T, vBC_B)/dt;
pressure{i} = CG_stage2(RHS_b) ;
% fractial step: stage 3 (assemble u_new)
u{i+1} = uf{1} - dt * Grad(cell2mat(pressure{i}));
s = [nx ny];
uval = u{end}(1:((nx-1) * ny)); % take only horizontal velocity values
ugraph = reshape(uval, [ny (nx-1)]);
%ugraph = [uBC_L ; ugraph; uBC_R];
vval = u{end}((nx-1)*ny + 1:end);
vgraph = reshape(vval, [(nx-1) (ny)]);
% plot
figure(1)
contourf(x, y, ugraph)
% N = 1000; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
% xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
% figure(i)
% xy{i} = stream2(x, y(1:end-1), ugraph(1:end-1, :), vgraph(:, 1:end-1), xstart, ystart,[0.3, 50]);
% streamlinegraph = streamline(xy{i});
%
i
end
%save('latest_solution.mat', 'u', 'p', ???)
%% testing online code
% close all
uval = u{end}(1:((nx-1) * ny));
vval = u{end}((nx-1)*ny + 1:end);
N = 1000; xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
xstart = max(x)*rand(N,1); ystart = max(y)*rand(N,1);
ugraph = reshape(uval, [ny (nx-1)]);
vgraph = reshape(vval, [(nx-1) ny]);
figure(2)
h = streamline(x, y(1:end-1), ugraph(1:end-1, :), vgraph(:, 1:end-1), xstart, ystart,[0.3, 50])
%% Req'ed functions
% ========================================================================
function [ip, iu, iv] = GenPointer(nx, ny)
% Memory allocation
ip = nan(nx, ny) ;
iu = nan(nx, ny) ;
iv = nan(nx, ny) ;
% Pointer matrix for P
id_p = 1 ; % index to be used in vector variable P
for i = 1:1:nx
for j = 1:1:ny
ip(i, j) = id_p ;
id_p = id_p + 1 ;
end
end
% Pointer matrix for ux
id_u = 1 ; % index to be used in vector variable u = [ux; uy]
for i = 2:1:nx
for j = 1:1:ny
iu(i, j) = id_u ;
id_u = id_u + 1 ;
end
end
% Pointer matrix for uy
for i = 1:1:nx
for j = 2:1:ny
iv(i, j) = id_u ;
id_u = id_u + 1 ;
end
end
end
% ========================================================================
function qo = Grad(qi)
global np nu nx ny dx dy iu iv ip
% Gradient operator: Pressure gradient at u and v points
% input: p-type (np elements)
% output: u-type (nu elements)
% Input size check:
if size(qi, 2) ~= 1
error('Grad:: input is not a column vector')
end
if size(qi, 1) ~= np
error('Grad:: input vector size mismatch!!!')
end
% Initialize output
qo = nan(nu, 1) ;
% inner domain
% x-direction gradient
for i = 2:1:nx
for j = 1:1:ny
qo(iu(i, j)) = ( -qi(ip(i-1, j)) + qi(ip(i, j)) ) / dx ;
end
end
% y-direction gradient
for i = 1:1:nx
for j = 2:1:ny
qo(iv(i, j)) = ( -qi(ip(i, j-1)) + qi(ip(i, j)) ) / dy ;
end
end
end
% ========================================================================
function qo = Div(qi)
global np nu nx ny dx dy iu iv ip
% divergence operator: To find velocity at p values
% input: u-type (nu elements)
% output: p-type (np elements)
% Input size check:
if size(qi, 2) ~= 1
error('Div:: input is not a column vector')
end
if size(qi, 1) ~= nu
error('Div:: input vector size mismatch!!!')
end
% Initialize output
qo = nan(np, 1) ;
% inner domain
for i = 2:1:(nx-1)
for j = 2:1:(ny-1)
qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
+ ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
end
end
% Edges
% bottom inner
j = 1 ;
for i = 2:1:(nx-1)
qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
+ ( qi(iv(i, j+1)) ) / dy ; % - qi(iv(i, j)) +
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
qo(ip(i, j)) = ( - qi(iu(i, j)) + qi(iu(i+1, j)) ) / dx ...
+ ( - qi(iv(i, j)) ) / dy ; % + qi(iv(i, j+1))
end
% left inner
i = 1 ;
for j = 2:1:(ny-1)
qo(ip(i, j)) = ( + qi(iu(i+1, j)) ) / dx ... % - qi(iu(i, j))
+ ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
end
% right inner
i = nx ;
for j = 2:1:(ny-1)
qo(ip(i, j)) = ( - qi(iu(i, j)) ) / dx ... % + qi(iu(i+1, j))
+ ( - qi(iv(i, j)) + qi(iv(i, j+1)) ) / dy ;
end
% Corners
% bottom left (pinning)
i = 1 ;
j = 1 ;
qo(ip(i, j)) = 0 ;
% bottom right
i = nx ;
j = 1 ;
qo(ip(i, j)) = ( - qi(iu(i, j)) ) / dx ... % + qi(iu(i+1, j))
+ ( + qi(iv(i, j+1)) ) / dy ; % - qi(iv(i, j))
% top left
i = 1 ;
j = ny ;
qo(ip(i, j)) = ( + qi(iu(i+1, j)) ) / dx ... % - qi(iu(i, j))
+ ( - qi(iv(i, j)) ) / dy ; % + qi(iv(i, j+1))
% top right
i = nx ;
j = ny ;
qo(ip(i, j)) = ( - qi(iu(i, j)) ) / dx ... % + qi(iu(i+1, j))
+ ( - qi(iv(i, j)) ) / dy ; % + qi(iv(i, j+1))
end
% ========================================================================
function bcD = BC_Div(uBC_L, uBC_R, vBC_T, vBC_B)
global np ip nx ny dx dy
% BC vector for divergence operator:
% input: BCs
% output: p-type (np elements)
% Initialize output
bcD = zeros(np, 1) ;
% Edges
% bottom inner
j = 1 ;
for i = 2:1:(nx-1)
bcD(ip(i, j)) = - vBC_B / dy ; % qi(iv(i, j+1))
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
bcD(ip(i, j)) = vBC_T / dy ; % qi(iv(i, j+1))
end
% left inner
i = 1 ;
for j = 2:1:(ny-1)
bcD(ip(i, j)) = - uBC_L(j)/ dx ;
end
% right inner
i = nx ;
for j = 2:1:(ny-1)
bcD(ip(i, j)) = uBC_R(j)/ dx ;
end
% Corners
% bottom left (need pinning)
i = 1 ;
j = 1 ;
bcD(ip(i, j)) = 0 ;
% bottom right
i = nx ;
j = 1 ;
bcD(ip(i, j)) = + uBC_R(j) / dx ...
- vBC_B / dy ;
% top left
i = 1 ;
j = ny ;
bcD(ip(i, j)) = - uBC_L(j) / dx ... % - qi(iu(i, j))
+ vBC_T / dy ; % + qi(iv(i, j+1))
% top right
i = nx ;
j = ny ;
bcD(ip(i, j)) = + uBC_R(j) / dx ...
+ vBC_T / dy ;
end
% ========================================================================
function qo = Laplace(qi)
global nu iu iv nx ny dx dy
% Laplace operator:
% input: u-type (nu elements)
% output: u-type (nu elements)
% Input size check:
if size(qi, 2) ~= 1
error('Div:: Need a colume vector input!!')
end
if size(qi, 1) ~= nu
error('Div:: input vector size mismatch!!!!!')
end
% Initialize output
qo = nan(nu, 1) ;
%% 1. ex-Component
% inner domain
for i = 3:1:(nx-1)
for j = 2:1:(ny-1)
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
end
end
% Edges
% left inner
i = 2 ;
for j = 2:1:(ny-1)
qo(iu(i, j)) = ( -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
end
% bottom inner
j = 1 ;
for i = 3:1:(nx-1)
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
+ ( -qi(iu(i, j)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
end
% right inner
i = nx ;
for j = 2:1:(ny-1)
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) ) / (dx^2) ... % + uBC_R / (dx^2)
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ;
end
% top inner
j = ny ;
for i = 3:1:(nx-1)
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ...
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j)) ) / (dy^2) ; % + 2*uBC_T / (dy^2)
end
% Corners
% bottom left
i = 2 ;
j = 1 ;
qo(iu(i, j)) = ( -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( -qi(iu(i, j)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
% bottom right
i = nx ;
j = 1 ;
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) ) / (dx^2) ... % + uBC_R / (dx^2)
+ ( -qi(iu(i, j)) -2*qi(iu(i, j)) + qi(iu(i, j+1)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
% top left
i = 2 ;
j = ny ;
qo(iu(i, j)) = ( -2*qi(iu(i, j)) + qi(iu(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j)) ) / (dy^2) ; % + 2*uBC_T / (dy^2)
% top right
i = nx ;
j = ny ;
qo(iu(i, j)) = ( +qi(iu(i-1, j)) -2*qi(iu(i, j)) ) / (dx^2) ... % + uBC_R / (dx^2)
+ ( +qi(iu(i, j-1)) -2*qi(iu(i, j)) - qi(iu(i, j)) ) / (dy^2) ; % + 2*uBC_T / (dy^2)
%% 2. ey-Component
% inner domain
for i = 2:1:(nx-1)
for j = 3:1:(ny-1)
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ...
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;
end
end
% Edges
% left inner
i = 1 ;
for j = 3:1:(ny-1)
qo(iv(i, j)) = ( -qi(iv(i, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + 2*vBC_L / (dx^2)
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ;
end
% bottom inner
j = 2 ;
for i = 2:1:(nx-1)
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ...
+ ( -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ; % + vBC_B
end
% right inner
i = nx ;
for j = 3:1:(ny-1)
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j)) ) / (dx^2) ... % + vBC_R
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ; % + vBC_B
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + vBC_R
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) ) / (dy^2) ; % + vBC_B
end
% Corners
% bottom left
i = 1 ;
j = 2 ;
qo(iv(i, j)) = ( -qi(iv(i, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
% bottom right
i = nx ;
j = 2 ;
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( -2*qi(iv(i, j)) + qi(iv(i, j+1)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
% top left
i = 1 ;
j = ny ;
qo(iv(i, j)) = ( -qi(iv(i, j)) -2*qi(iv(i, j)) + qi(iv(i+1, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
% top right
i = nx ;
j = ny ;
qo(iv(i, j)) = ( +qi(iv(i-1, j)) -2*qi(iv(i, j)) - qi(iv(i, j)) ) / (dx^2) ... % + uBC_L / (dx^2)
+ ( +qi(iv(i, j-1)) -2*qi(iv(i, j)) ) / (dy^2) ; % + 2*uBC_B / (dy^2)
end
% ========================================================================
function bcL = BC_Laplace()
global nu iu iv nx ny dx dy uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B
% BC vector for divergence operator:
% input: BCs
% output: u-type (nu elements)
% Input size check:
% input BC's are all scalars
% Initialize output
bcL = zeros(nu, 1) ;
%% 1. U-Component
% inner domain
% Edges
% left inner
i = 2 ;
for j = 2:1:(ny-1)
bcL(iu(i, j)) = + uBC_L(j) / (dx^2) ;
end
% bottom inner
j = 1 ;
for i = 3:1:(nx-1)
bcL(iu(i, j)) = + 2*uBC_B / (dy^2) ;
end
% right inner
i = nx ;
for j = 2:1:(ny-1)
bcL(iu(i, j)) = + uBC_R(j) / (dx^2) ;
end
% top inner
j = ny ;
for i = 3:1:(nx-1)
bcL(iu(i, j)) = + 2*uBC_T / (dy^2) ;
end
% Corners
% bottom left
i = 2 ;
j = 1 ;
bcL(iu(i, j)) = + uBC_L(j) / (dx^2) + 2*uBC_B / (dy^2) ;
% bottom right
i = nx ;
j = 1 ;
bcL(iu(i, j)) = + uBC_R(j) / (dx^2) + 2*uBC_B / (dy^2) ;
% top left
i = 2 ;
j = ny ;
bcL(iu(i, j)) = + uBC_L(j) / (dx^2) + 2*uBC_T / (dy^2) ;
% top right
i = nx ;
j = ny ;
bcL(iu(i, j)) = + uBC_R(j) / (dx^2) + 2*uBC_T / (dy^2) ;
%% 2. V-Component
% inner domain
% Edges
% bottom inner
j = 2 ;
for i = 2:1:(nx-1)
bcL(iv(i, j)) = + vBC_B / (dy^2) ;
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
bcL(iv(i, j)) = + vBC_T / (dy^2) ;
end
% left inner
i = 1 ;
for j = 3:1:(ny-1)
bcL(iv(i, j)) = + 2*vBC_L / (dx^2) ;
end
% right inner
i = nx ;
for j = 3:1:(ny-1)
bcL(iv(i, j)) = + 2*vBC_R / (dx^2) ;
end
% Corners
% bottom left
i = 1 ;
j = 2 ;
bcL(iv(i, j)) = + 2*vBC_L / (dx^2) + vBC_B / (dy^2) ;
% bottom right
i = nx ;
j = 2 ;
bcL(iv(i, j)) = + 2*vBC_R / (dx^2) + vBC_B / (dy^2) ;
% top left
i = 1 ;
j = ny ;
bcL(iv(i, j)) = + 2*vBC_L / (dx^2) + vBC_T / (dy^2) ;
% top right
i = nx ;
j = ny ;
bcL(iv(i, j)) = + 2*vBC_R / (dx^2) + vBC_T / (dy^2) ;
end
% ========================================================================
function qo = Adv(qi)
global nu iu iv nx ny dx dy uBC_L uBC_R uBC_B uBC_T vBC_L vBC_R vBC_T vBC_B
% advection operator (BC embedded): -\nabla \cdot (uu)
% input: u-type (nu elements)
% output: u-type (nu elements)
%
% Input size check:
if size(qi, 2) ~= 1
error('Adv:: Need a colume vector input!!')
end
if size(qi, 1) ~= nu
error('Adv:: input vector size mismatch!!')
end
% Initialize output
qo = nan(nu, 1) ;
% 1. U-Component
% inner domain
for i = 3:1:(nx-1)
for j = 2:1:(ny-1)
qo(iu(i, j)) = - (1/dx) * ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) ;
end
end
% Edges
% left inner
i = 2 ;
for j = 2:1:(ny-1)
qo(iu(i, j)) = - ( - ( uBC_L(j) + qi(iu(i ,j )) ) / 2 * ( uBC_L(j) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) / dx ...
- ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) / dy ;
end
% bottom inner
j = 1 ;
for i = 3:1:(nx-1)
qo(iu(i, j)) = - ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) / dx ...
- ( - uBC_B * vBC_B ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) / dy ;
end
% right inner
i = nx ;
for j = 2:1:(ny-1)
qo(iu(i, j)) = - ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + uBC_R(j) ) / 2 * ( qi(iu(i ,j )) + uBC_R(j) ) / 2 ) / dx ...
- ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) / dy ;
end
% top inner
j = ny ;
for i = 3:1:(nx-1)
qo(iu(i, j)) = - ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) / dx ...
- ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_T * vBC_T ) / dy ;
end
% Corners
% bottom left
i = 2 ;
j = 1 ;
qo(iu(i, j)) = - ( - ( uBC_L(j) + qi(iu(i ,j )) ) / 2 * ( uBC_L(j) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) / dx ...
- ( - uBC_B * vBC_B ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) / dy ;
% bottom right
i = nx ;
j = 1 ;
qo(iu(i, j)) = - ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + uBC_R(j) ) / 2 * ( qi(iu(i ,j )) + uBC_R(j) ) / 2 ) / dx ...
- ( - uBC_B * vBC_B ...
+ ( qi(iu(i ,j )) + qi(iu(i ,j+1)) ) / 2 * ( qi(iv(i-1,j+1)) + qi(iv(i ,j+1)) ) / 2 ) / dy ;
% top left
i = 2 ;
j = ny ;
qo(iu(i, j)) = - ( - ( uBC_L(j) + qi(iu(i ,j )) ) / 2 * ( uBC_L(j) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 * ( qi(iu(i ,j )) + qi(iu(i+1,j )) ) / 2 ) / dx ...
- ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_T * vBC_T ) / dy ;
% top right
i = nx ;
j = ny ;
qo(iu(i, j)) = - ( - ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 * ( qi(iu(i-1,j )) + qi(iu(i ,j )) ) / 2 ...
+ ( qi(iu(i ,j )) + uBC_R(j) ) / 2 * ( qi(iu(i ,j )) + uBC_R(j) ) / 2 ) / dx ...
- ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_T * vBC_T ) / dy ;
% 2. V-Component
% inner domain
for i = 2:1:(nx-1)
for j = 3:1:(ny-1)
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
end
end
% Edges
% left inner
i = 1 ;
for j = 3:1:(ny-1)
qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j) * vBC_L ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
end
% bottom inner
j = 2 ;
for i = 2:1:(nx-1)
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( vBC_B + qi(iv(i ,j )) ) / 2 * ( vBC_B + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
end
% right inner
i = nx ;
for j = 3:1:(ny-1)
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_R(j) * vBC_R ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
end
% top inner
j = ny ;
for i = 2:1:(nx-1)
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + vBC_T ) / 2 * ( qi(iv(i ,j )) + vBC_T ) / 2 ) ;
end
% Corners
% bottom left
i = 1 ;
j = 2 ;
qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j) * vBC_L ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( vBC_B + qi(iv(i ,j )) ) / 2 * ( vBC_B + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
% bottom right
i = nx ;
j = 2 ;
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_R(j) * vBC_R ) ...
- (1/dy) * ( - ( vBC_B + qi(iv(i ,j )) ) / 2 * ( vBC_B + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i ,j+1)) ) / 2 ) ;
% top left
j = ny ;
i = 1 ;
qo(iv(i, j)) = - (1/dx) * ( - uBC_L(j) * vBC_L ...
+ ( qi(iu(i+1,j-1)) + qi(iu(i+1,j )) ) / 2 * ( qi(iv(i ,j )) + qi(iv(i+1,j )) ) / 2 ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + vBC_T ) / 2 * ( qi(iv(i ,j )) + vBC_T ) / 2 ) ;
% top right
j = ny ;
i = nx ;
qo(iv(i, j)) = - (1/dx) * ( - ( qi(iu(i ,j-1)) + qi(iu(i ,j )) ) / 2 * ( qi(iv(i-1,j )) + qi(iv(i ,j )) ) / 2 ...
+ uBC_R(j) * vBC_R ) ...
- (1/dy) * ( - ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 * ( qi(iv(i ,j-1)) + qi(iv(i ,j )) ) / 2 ...
+ ( qi(iv(i ,j )) + vBC_T ) / 2 * ( qi(iv(i ,j )) + vBC_T ) / 2 ) ;
end
% ========================================================================
function qout = R(qin)
global dt noo
qout = qin - .5 * dt * noo * Laplace(qin);
end
% ========================================================================
function qout = Rinverse(qin)
global dt noo
qout = qin + .5 * dt * noo * Laplace(qin);
end
% ========================================================================
function x_sol = CG_EE(b)
global np
x(:, 1) = rand(np, 1) ;
n = 1;
r(:, n) = b - Div(Grad(x(:, 1)));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
alpha = (r(:, n)' * r(:, n)) / (d' * Div(Grad(d)));
x(:, n+1) = x(:, n) + alpha * d;
r(:, n+1) = r(:, n) - alpha * Div(Grad(d));
beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
d = r(:, n+1) + beta * d;
error = norm(r(:, n+1) - r(:, n));
n = n+1;
x_sol{n} = x(:, n);
end
x_sol = x_sol(end);
end
% ========================================================================
function x_sol = CG_stage1(b)
global nu
x(:, 1) = zeros(nu, 1) ;
n = 1;
r(:, n) = b - R(x(:, n));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
alpha = (r(:, n)' * r(:, n)) / (d' * R(d));
x(:, n+1) = x(:, n) + alpha * d;
r(:, n+1) = r(:, n) - alpha * R(d);
beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
d = r(:, n+1) + beta * d;
error = norm(r(:, n+1) - r(:, n));
n = n+1;
x_sol{n} = x(:, n);
end
x_sol = x_sol(end);
end
% ========================================================================
function x_sol = CG_stage2(b)
global nx ny
x(:, 1) = zeros(nx * ny, 1) ;
n = 1;
r(:, n) = b - Div(Grad(x(:, n)));
d = r(:, n);
error = 1;
x_sol{n} = x(:, n);
while error > .001
alpha = (r(:, n)' * r(:, n)) / (d' * Div(Grad(d)));
x(:, n+1) = x(:, n) + alpha * d;
r(:, n+1) = r(:, n) - alpha * Div(Grad(d));
beta = (r(:, n+1)' * r(:, n+1))/(r(:, n)' * r(:, n));
d = r(:, n+1) + beta * d;
error = norm(r(:, n+1) - r(:, n));
n = n+1;
x_sol{n} = x(:, n);
end
x_sol = x_sol(end);
end
回答(1 个)
Pratyush Swain
2023-12-11
编辑:Pratyush Swain
2023-12-11
Hi Shaan,
I understand you are unable to clear the the previous iteration of streamlines as it keeps layering on top of the previous ones. One workaround is to use the "cla" command to clear the axes which will delete all visible figures in each iteration. Please refer to the example implementation below:
% Load Dataset %
load wind
% Meshgrid for starting points %
[startX,startY] = meshgrid(80,20:10:50);
% Cell array to hold the streamline object handles %
lineobjs = cell(1,10);
% Looping over new set of data points in each iteration %
for i=1:10
xi = x(:,:,i);
yi = y(:,:,i);
ui = u(:,:,i);
vi = v(:,:,i);
% Compute the 2-D streamline vertex data %
verts = stream2(xi,yi,ui,vi,startX,startY);
% Visualize the 2-D matrix of vector fields by calling streamline function %
lineobjs{i} = streamline(verts);
% Add title and labels %
title(['Streamlines at time step ' num2str(i)]);
xlabel('X');
ylabel('Y');
% Add pause and then clear the figure object %
pause(0.5)
cla
end
The "cla" command along with "pause" function will help to visualize the streamlines over each iteration without any overlapping.For more information please refer to
Hope this helps.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Audio I/O and Waveform Generation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)