Skip to content

Commit

Permalink
first commit
Browse files Browse the repository at this point in the history
  • Loading branch information
cfd1 committed Sep 3, 2015
0 parents commit 70bb9b9
Show file tree
Hide file tree
Showing 12 changed files with 1,158 additions and 0 deletions.
250 changes: 250 additions & 0 deletions BICGSTABIter_p.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,250 @@
function [p_int] = BICGSTABIter_p(u,r,Nr,Np,hr,hp,f)

%ip = Nr-1 ; jp = Np/2 ;

STOP = false ;

BICGEPS = 1.0e-16 ;

Temp=ComputeAX_p(u,r,Nr,Np,hr,hp) ;

% Initial vector r_0 = b - Ax_0, and r0* = r_0

for i = 2:(Nr-1)
for j =1:Np

Uj(i,j) = u(i,j) ;
rj(i,j) = f(i,j) - Temp(i,j) ;
r0_star(i,j) = rj(i,j) ;

%{
if (i == ip) && (j == jp)
Uj(i,j) = 0.0 ;
rj(i,j) = 0.0 ;
r0_star(i,j) = 0.0 ;
end
%}

end
end

BICG_ITER = 0 ; norm = 0.0 ;

while ((BICG_ITER < 10000 ) && (~STOP == 1))

% compute rhoj = (r0, r0*)

rhoj = 0.0 ;
for i = 2:(Nr-1)
for j =1:Np
rhoj = rhoj + rj(i,j)*r0_star(i,j) ;
end
end

if ( sqrt(rhoj/((Nr-1)*(Np+1))) < BICGEPS )

STOP = true ;

else

if ( BICG_ITER == 0 )

for i = 2:(Nr-1)
for j = 1:Np
pj(i,j) = rj(i,j); % p0 = r0
end
end

else
betaj = (rhoj/rhoj_Minus)*(alphaj/omegaj) ;

for i = 1:(Nr-1)
for j = 1:Np
pj(i,j) = rj(i,j) + betaj*(pj(i,j) - omegaj*Var(i,j));
end
end
end

%Solve for Upstar, Vpstar from Ku* = u...., where K is the preconditioning matrix

for i = 2:(Nr-1)
for j = 1:Np
% No preconditioning
pstar(i,j) = pj(i,j) ;
end
end

% compute vj = A*pstar

for i = 2:(Nr-1)
for j = 1:Np

%{
if (i == ip) && (j == jp)
u(i,j) = 0.0 ;
else
u(i,j) = pstar(i,j) ;
end
%}


u(i,j) = pstar(i,j) ;

end
end

Temp=ComputeAX_p(u,r,Nr,Np,hr,hp) ;

for i = 2:(Nr-1)
for j = 1:Np
Var(i,j) = Temp(i,j) ;
end
end
H1 = 0.0 ;
for i = 2:(Nr-1)
for j = 1:Np
H1 = H1 + Var(i,j)*r0_star(i,j) ;
end
end
alphaj = rhoj/H1 ;

% find sj

for i = 2:(Nr-1)
for j = 1:Np
sj(i,j) = rj(i,j) - alphaj*Var(i,j) ;
end
end

% Solve for Upstar, Vpstar from Ku* = u...., where K is the preconditioning matrix

for i = 2:(Nr-1)
for j = 1:Np
% No preconditioning
sstar(i,j) = sj(i,j) ;
end
end
norm = 0.0 ;
for i = 2:(Nr-1)
for j = 1:Np
norm = norm + sstar(i,j)*sstar(i,j) ;
end
end
norm = sqrt(norm/((Nr-1)*(Np+1))) ;
if( norm < BICGEPS)
STOP = true ; % if ||s||_2 is small x_i = x_{i-1} + alphai*p_i
for i = 2:(Nr-1)
for j = 1:Np

%{
if (i == ip) && (j == jp)
Uj(i,j) = 0.0 ;
else
Uj(i,j) = Uj(i,j)+ alphaj*pstar(i,j) ;
end
%}
Uj(i,j) = Uj(i,j)+ alphaj*pstar(i,j) ;
end
end
else

% compute t = As

for i = 2:(Nr-1)
for j = 1:Np

%{
if (i == ip) && (j == jp)
u(i,j) = 0.0 ;
else
u(i,j) = sstar(i,j) ;
end
%}
u(i,j) = sstar(i,j) ;
end
end

Temp=ComputeAX_p(u,r,Nr,Np,hr,hp) ;

H1 = 0.0 ; H2 = 0.0 ;
for i = 2:(Nr-1)
for j = 1:Np
H1 = H1 + Temp(i,j)*sj(i,j) ;
H2 = H2 + Temp(i,j)*Temp(i,j) ;
end
end
omegaj = H1/H2;

% find xj
norm = 0.0 ;
for i = 2:(Nr-1)
for j = 1:Np
H1 = (alphaj*pstar(i,j) + omegaj*sstar(i,j)) ;
%{
if (i == ip) && (j == jp)
Uj(i,j) = 0.0 ;
else
Uj(i,j) = Uj(i,j) + H1 ;
end
%}
Uj(i,j) = Uj(i,j) + H1 ;
norm = norm + H1*H1 ;
end
end
norm = sqrt(norm/((Nr-1)*(Np+1))) ;

if(norm < BICGEPS)
STOP = true ;
end

% find rjplusone

for i = 2:(Nr-1)
for j = 1:Np
rj(i,j) = sj(i,j) - omegaj*Temp(i,j);
end
end
rhoj_Minus = rhoj ;
end
end

BICG_ITER = BICG_ITER + 1 ;

end

for i = 2:(Nr-1)
for j = 1:Np
p_int(i,j) = Uj(i,j) ;
end
end

fprintf('\n Pressure converged at Iter %d, ', BICG_ITER);

end

59 changes: 59 additions & 0 deletions ComputeAX_p.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
function [Temp] = ComputeAX_p(u,r,Nr,Np,hr,hp)

%ip = Nr-1 ; jp = Np/2 ;

for i = 2:(Nr-1)
for j =1:Np

if i==2

u_r = (u(i+1,j) - u(i,j))/(2.0*hr) ;
u_rr = (u(i+1,j) - u(i,j))/(hr*hr) ;

elseif i==(Nr-1)

u_r = (u(i,j) - u(i-1,j))/(2.0*hr) ;
u_rr = (-u(i,j) + u(i-1,j))/(hr*hr) ;

else

u_r = (u(i+1,j) - u(i-1,j))/(2.0*hr) ;
u_rr = (u(i+1,j) - 2.0*u(i,j) + u(i-1,j))/(hr*hr) ;

end

if j==1

u_pp = (u(i,j+1) - 2.0*u(i,j) + u(i,Np))/(hp*hp) ;

elseif j==Np

u_pp = (u(i,1) - 2.0*u(i,j) + u(i,j-1))/(hp*hp) ;

else

u_pp = (u(i,j+1) - 2.0*u(i,j) + u(i,j-1))/(hp*hp) ;

end


%{
if i==ip && j==jp
Temp(i,j) = 0 ;
else
Temp(i,j) = u_rr + (1.0/r(i))*u_r + (1.0/(r(i)*r(i)))*u_pp ;
end
%}

Temp(i,j) = u_rr + (1.0/r(i))*u_r + (1.0/(r(i)*r(i)))*u_pp ;

end
end

end

66 changes: 66 additions & 0 deletions ComputeAx_collocatedUs2.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@

% ComputeAx
function R=ComputeAx_collocatedUs2(X, rr, tt ,delr, delth, r_n, t_n,dt,Re)
% int ii r direction,jj theta direction; double P_rr, P_tt, P_r
%X=reshape(X1,r_n,t_n);

R=zeros(r_n,t_n);
for ii=2:r_n-1 %r direction
for jj=1:t_n % theta direction
% for theta derivative
if(jj==1)

Uss_tt = (X(ii,jj+1)+X(ii,t_n)-2*X(ii,jj))/(delth^2);

elseif (jj==t_n)
Uss_tt = (X(ii,1)+X(ii,jj-1)-2*X(ii,jj))/(delth^2);
else
Uss_tt = (X(ii,jj+1)+X(ii,jj-1)-2*X(ii,jj))/(delth^2);
end

% for first r derivatives
if (ii==2)
Uss_r=(X(ii+1,jj)-0)/(2*delr);% dirichlet
elseif (ii==r_n-1)
% At the farfield
% if (tt(ii,jj)<pi/2 || tt(ii,jj)>3*pi/2 ) % neumann
% Uss_r=(X(ii,jj)-X(ii-1,jj))/(2*delr);
% % downstream
% else

Uss_r=(0-X(ii-1,jj))/(2*delr);% Dirichlet
% end

else
Uss_r=(X(ii+1,jj)-X(ii-1,jj))/(2*delr); %
end

% for second r derivatives
if (ii==2)
Uss_rr=(X(ii+1,jj)-2*X(ii,jj) + 0)/(delr^2);

elseif (ii==r_n-1)

% if (tt(ii,jj)<pi/2 || tt(ii,jj)>3*pi/2 )
% Uss_rr=(-1*X(ii,jj)+X(ii-1,jj))/(delr^2);

% else


Uss_rr=(0-2*X(ii,jj)+X(ii-1,jj))/(delr^2);% Dirichlet
% end

else
Uss_rr=(X(ii+1,jj)-2*X(ii,jj)+X(ii-1,jj))/(delr^2);
end


R(ii,jj)= (dt/Re)*( Uss_rr + (1/(rr(ii,jj)^2))*Uss_tt + (3/rr(ii,jj))*Uss_r ) ...
+((dt/(Re*rr(ii,jj)^2)) -1)*X(ii,jj);

end
end
%R(r_n-1,1:2)=X(r_n-1,1:2);
%R=R1(:);

end
Loading

0 comments on commit 70bb9b9

Please sign in to comment.