-
Notifications
You must be signed in to change notification settings - Fork 0
/
Vcycle.m
66 lines (59 loc) · 2.75 KB
/
Vcycle.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
%% Recursive V-Cycles to improve approximation of solutions
function [Psih,Wh] = Vcycle(Psih,Wh,Sph,Swh,BC,h,Re,a)
%% Condition current grid
% Relax
[Psih,Wh] = Relax(Psih,Wh,Sph,Swh,BC,h,Re,6,a);
% Compute Residual Errors
[Presh,Wresh] = residual(Psih,Wh,Sph,Swh,h,Re);
%% Restrict to Coarse grid
Psi2h = restrict(Psih); W2h = restrict(Wh);
Pres2h = restrict(Presh); Wres2h = restrict(Wresh);
% Update BCs
N2h = size(Psi2h); N2x = N2h(2)-2; N2y = N2h(1)-2; hii = 1/(2*h); [Pxxt2h,~] = BC(2*h,Re);
W2h(N2y+1,2:N2x+1) = 2*(Psi2h(N2y+1,2:N2x+1)-Psi2h(N2y,2:N2x+1))*hii^2 - Pxxt2h(:)';
W2h(2:N2y+1,2) = -2*Psi2h(2:N2y+1,3)*hii^2;
W2h(3:N2y,N2x+1) = (2*Psi2h(3:N2y,N2x+1)-Psi2h(2:N2y-1,N2x)-Psi2h(4:N2y+1,N2x))*hii^2;
% Extrapolated points
Psi2h(2:N2y+1,1) = Psi2h(2:N2y+1,3);
Psi2h(2:N2y+1,N2x+2) = Psi2h(2:N2y+1,N2x);
Psi2h(1,2:N2x+1) = 3*Psi2h(2,2:N2x+1) - 3*Psi2h(3,2:N2x+1) + Psi2h(4,2:N2x+1);
Psi2h(N2y+2,2:N2x+1) = Psi2h(N2y,2:N2x+1);
W2h(2:N2y+1,1) = 3*W2h(2:N2y+1,2) - 3*W2h(2:N2y+1,3) + W2h(2:N2y+1,4);
W2h(2:N2y+1,N2x+2) = 3*W2h(2:N2y+1,N2x+1) - 3*W2h(2:N2y+1,N2x) + W2h(2:N2y+1,N2x-1);
W2h(1,2:N2x+1) = -3*W2h(3,2:N2x+1) + W2h(4,2:N2x+1);
W2h(N2y+2,2:N2x+1) = 3*W2h(N2y+1,2:N2x+1) - 3*W2h(N2y,2:N2x+1) + W2h(N2y-1,2:N2x+1);
% form new RHS
[Sp2h,Sw2h] = RHS(Psi2h,W2h,Pres2h,Wres2h,2*h,Re);
% stop recursion at smallest grid size, otherwise continue recursion
if h>0.8
[Ap2h,Aw2h] = Relax(Psi2h,W2h,Sp2h,Sw2h,BC,2*h,Re,15,0.95);
else
[Ap2h,Aw2h] = Vcycle(Psi2h,W2h,Sp2h,Sw2h,BC,2*h,Re,a);
end
%% Interpolate upto fine grid
% determine error
Ep2h = Ap2h - Psi2h;
Ew2h = Aw2h - W2h;
% interpolate error
Eph = interp(Ep2h);
Ewh = interp(Ew2h);
% correct
Psih = Psih + Eph;
Wh = Wh + Ewh;
% Update BCs
N = size(Psih); Nx = N(2)-2; Ny = N(1)-2; hi = 1/h; [Pxxth,~] = BC(h,Re);
Wh(Ny+1,2:Nx+1) = 2*(Psih(Ny+1,2:Nx+1)-Psih(Ny,2:Nx+1))*hi^2 - Pxxth(:)';
Wh(2:Ny+1,2) = -2*Psih(2:Ny+1,3)*hi^2;
Wh(3:Ny,Nx+1) = (2*Psih(3:Ny,Nx+1)-Psih(2:Ny-1,Nx)-Psih(4:Ny+1,Nx))*hi^2;
% Extrapolated points
Psih(2:Ny+1,1) = Psih(2:Ny+1,3);
Psih(2:Ny+1,Nx+2) = Psih(2:Ny+1,Nx);
Psih(1,2:Nx+1) = 3*Psih(2,2:Nx+1) - 3*Psih(3,2:Nx+1) + Psih(4,2:Nx+1);
Psih(Ny+2,2:Nx+1) = Psih(Ny,2:Nx+1);
Wh(2:Ny+1,1) = 3*Wh(2:Ny+1,2) - 3*Wh(2:Ny+1,3) + Wh(2:Ny+1,4);
Wh(2:Ny+1,Nx+2) = 3*Wh(2:Ny+1,Nx+1) - 3*Wh(2:Ny+1,Nx) + Wh(2:Ny+1,Nx-1);
Wh(1,2:Nx+1) = -3*Wh(3,2:Nx+1) + Wh(4,2:Nx+1);
Wh(Ny+2,2:Nx+1) = 3*Wh(Ny+1,2:Nx+1) - 3*Wh(Ny,2:Nx+1) + Wh(Ny-1,2:Nx+1);
% Post-Smoothing
[Psih,Wh] = Relax(Psih,Wh,Sph,Swh,BC,h,Re,3,a);
end