-
Notifications
You must be signed in to change notification settings - Fork 3
/
solveVs22.m
134 lines (102 loc) · 2.85 KB
/
solveVs22.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
function solveVs22(Vs)
global Us2;
global Vs2;
global rr;
global tt;
global delr;
global delth;
global r_n;
global t_n;
global Re;
global dt;
global k;
b1 =zeros(size(Vs2));
b =zeros(size(Vs2));
% Define the RHS
for i=2:r_n-1
for j=1:t_n
if i==r_n-1
% if (tt(i,j)<pi/2 || tt(i,j)>3*pi/2 ) % dirichlet
% b1(i,j) = -Vs(i,j);
% else
b1(i,j) = -Vs(i,j)-(dt/Re)*( (1/delr^2)*Vs(i+1,j) + (1/rr(i,j))*((1/(2*delr))*Vs(i+1,j)));
% end
else
b1(i,j) = -Vs(i,j) ;
end
end
end
for i=2:r_n-1
for j=1:t_n
if j==1
dutheta = (1/(2*delth))*(Us2(i,j+1)-Us2(i,t_n-1));
elseif j==t_n
dutheta=(1/(2*delth))*(Us2(i,2)-Us2(i,j-1));
else
dutheta = (1/(2*delth))*(Us2(i,j+1)-Us2(i,j-1)) ;
end
b(i,j) = b1(i,j) - (dt/Re)*(2/rr(i,j)^2)*dutheta;
end
end
% Biconjugate Gradient Stabilized method %%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%---------------Gonna try matlab bicgstab------------------
% Vectors for BICGSTAB
v0=zeros(size(Vs2));
d0=zeros(size(Vs2));
dd0=zeros(size(Vs2));
x0=eye(size(Vs2));%zeros(size(XX));
%x0(10,10)=1;
% Scalars
p0=v0;
alpha =1;
w0=1;
rho0=1;
iter = 0;
Ax0=ComputeAx_collocatedVs2(x0,rr,tt, delr, delth, r_n, t_n,dt,Re);
for i=2:r_n-1
for j=1:t_n
dd0(i,j)=b(i,j)-Ax0(i,j);
d0(i,j)=dd0(i,j);
end
end
%dd0=b-Ax0;
%d0=dd0;
Tol=1e-15;
Res=100;
while Res>Tol,
rho_i=innerpUV(dd0,d0,r_n,t_n);
if (rho_i==0) break, end
if (iter>0)
beta = (rho_i/rho0)*(alpha/w0);
p_i = d0 + beta.*(p0 - w0.*v0);
else
p_i=d0;
end
v_i = ComputeAx_collocatedVs2(p_i,rr,tt,delr,delth,r_n,t_n,dt,Re);
alpha = rho_i/innerpUV(dd0,v_i,r_n,t_n);
s = d0 - alpha.*v_i;
t=ComputeAx_collocatedVs2(s,rr,tt,delr,delth,r_n,t_n,dt,Re);
w_i=innerpUV(t,s,r_n,t_n)/innerpUV(t,t,r_n,t_n);
x_i=x0+alpha.*p_i + w_i.*s;
d_i=s - w_i.*t;
%P= x_i;
%%%%%%%%
d0=d_i;
Res = innerpUV((x_i - x0 ), (x_i -x0),r_n,t_n);
Res = sqrt ( Res / ((r_n-2)*(t_n)));
x0=x_i;
p0=p_i;
w0=w_i;
v0=v_i;
rho0=rho_i;
iter=iter+1;
%fprintf('\n Iter is %d Residual is %0.3g ', iter, Res);
end
for i=2:r_n-1
for j=1:t_n
Vs2(i,j) = x_i(i,j);
end
end
fprintf('\n V converged at Iter %d , Residual is %0.3g at time step %d', iter, Res, k);
end