Streamline function keeps overwriting plot in a for loop

6 Dec 2023
1 Answer
Updated 11 Dec 2023
8 Views (30 days)

0 votes

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
Show None Hide None

VBBV
VBBV on 6 Dec 2023
can you share your code ?
Shaan Stephen
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

Sign in to comment.

Sign in to answer this question.

Answers (1)

Pratyush Swain
Pratyush Swain on 11 Dec 2023
Edited: Pratyush Swain on 11 Dec 2023

0 votes

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
1.https://www.mathworks.com/help/matlab/ref/streamline.html
2.https://www.mathworks.com/help/matlab/ref/cla.html
Hope this helps.

Categories

Find more on Heat Transfer in Help Center and File Exchange

Tags

Asked:

Shaan Stephen
Shaan Stephen
on 6 Dec 2023

Edited:

Pratyush Swain
Pratyush Swain
on 11 Dec 2023

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!