You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Streamline function keeps overwriting plot in a for loop
2 views (last 30 days)
Show older comments
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 Comments
Shaan Stephen
on 6 Dec 2023
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
Answers (1)
Pratyush Swain
on 11 Dec 2023
Edited: Pratyush Swain
on 11 Dec 2023
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.
See Also
Categories
Find more on Audio I/O and Waveform Generation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- 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)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)