% Solution of Stokes and continuity equations in 3D with Multigrid
% based on V-cycle
% by using external function Stokes_smoother3D()
% Density distribution in the model corresponds to falling block test
% viscosity in the model is constant
% 
% Staggered Grid for Multigrid
% 
%     vx       vx       vx    
%
% vy  +---vy---+---vy---+   vy
%     |        |        |
%     vx   P   vx   P   vx    
%     |        |        |
% vy  +---vy---+---vy---+   vy
%     |        |        |
%     vx   P   vx   P   vx    
%     |        |        |
% vy  +---vy---+---vy---+   vy
%
%     vx       vx       vx    
% 
% Lines show basic grid
% Basic (density) nodes are shown with +
% Ghost nodes shown outside the basic grid
% are used for boundary conditions

% Clearing all variables and arrays
clear all;
% Clearing figures
clf;
% 1=inclusion, 2=solcx, 3=solkz 
% 4=linvisc2D, 5=expvisc2D 
% 6=linvisc3D, 7=expvisc3D 
gridtest=5;

% Resolution cycle
for r=1:1:6
% Number of resolution levels
levelnum=fix(r/2)+1;

% Iteration Parameters
% Total number of iteration cycles
inum=50;
% Relaxation coefficients for Gauss-Seidel iterations
% Stokes equations
krelaxs=1.25;
vkoef=1.0;
% Continuity equation
krelaxc=0.75;
pkoef=1.0;
% Acceleration of Gravity, m/s^2
gx=10; % Gravity, m/s^2
gy=10; % Gravity, m/s^2
gz=0; % Gravity, m/s^2
% Model size, m
xsize=1;
ysize=1;
zsize=1;
if(gridtest==1)
    gx=0;
    gy=0; % inclusion
    gz=0;
    % Model size, m
    xsize=2;
    ysize=2;
    zsize=2;
elseif(gridtest==2)
    gx=0;
    gy=1; % solcx
    gz=0;
    % Relaxation coefficients for Gauss-Seidel iterations
    % Stokes equations
    krelaxs=1.25;
    vkoef=1.0;
    % Continuity equation
    krelaxc=0.75;
    pkoef=1.0;
    % Number of MG levels
    levelnum=r+2;
elseif (gridtest==3)
    gx=0;
    gy=-1; % solkz
    gz=0;
    % Stokes equations
    krelaxs=1.25;
    vkoef=1.0;
    % Continuity equation
    krelaxc=0.75;
    pkoef=1.0;
    % Number of MG levels
    levelnum=r;
end


% Scales for residuals
rhoblock=3000;
etamedium=1;
g=1;


% Multigrid parameters
% Number of smoothing iterations
iternum(1)=5;
iternum(2)=5*2;
iternum(3)=5*4;
iternum(4)=5*8;
iternum(5)=5*16;
iternum(6)=5*32;
iternum(7)=5*64;
iternum(8)=5*128;


if (r==1)
% Defining resolutions on all levels
xnum(1)=13;
ynum(1)=13;
znum(1)=13;
xnum(2)=7;
ynum(2)=7;
znum(2)=7;
xnum(3)=4;
ynum(3)=4;
znum(3)=4;

elseif (r==2)
% Defining resolutions on all levels
xnum(1)=17;
ynum(1)=17;
znum(1)=17;
xnum(2)=11;
ynum(2)=11;
znum(2)=11;
xnum(3)=7;
ynum(3)=7;
znum(3)=7;
xnum(4)=4;
ynum(4)=4;
znum(4)=4;

elseif (r==3)
% Defining resolutions on all levels
xnum(1)=25;
ynum(1)=25;
znum(1)=25;
xnum(2)=13;
ynum(2)=13;
znum(2)=13;
xnum(3)=7;
ynum(3)=7;
znum(3)=7;
xnum(4)=4;
ynum(4)=4;
znum(4)=4;

elseif (r==4)
% Defining resolutions on all levels
xnum(1)=33;
ynum(1)=33;
znum(1)=33;
xnum(2)=17;
ynum(2)=17;
znum(2)=17;
xnum(3)=13;
ynum(3)=13;
znum(3)=13;
xnum(4)=7;
ynum(4)=7;
znum(4)=7;

elseif (r==5)
% Defining resolutions on all levels
xnum(1)=41;
ynum(1)=41;
znum(1)=41;
xnum(2)=31;
ynum(2)=31;
znum(2)=31;
xnum(3)=21;
ynum(3)=21;
znum(3)=21;
xnum(4)=13;
ynum(4)=13;
znum(4)=13;
xnum(5)=7;
ynum(5)=7;
znum(5)=7;

elseif (r==6)
% Defining resolutions on all levels
xnum(1)=49;
ynum(1)=49;
znum(1)=49;
xnum(2)=31;
ynum(2)=31;
znum(2)=31;
xnum(3)=21;
ynum(3)=21;
znum(3)=21;
xnum(4)=13;
ynum(4)=13;
znum(4)=13;
xnum(5)=7;
ynum(5)=7;
znum(5)=7;

elseif (r==7)
% Defining resolutions on all levels
xnum(1)=97;
ynum(1)=97;
znum(1)=97;
xnum(2)=49;
ynum(2)=49;
znum(2)=49;
xnum(3)=25;
ynum(3)=25;
znum(3)=25;
xnum(4)=13;
ynum(4)=13;
znum(4)=13;
xnum(5)=7;
ynum(5)=7;
znum(5)=7;
xnum(6)=4;
ynum(6)=4;
znum(6)=4;
end

% Defining gridsteps on all levels
% Number 
for n=levelnum:-1:1
    % Grid steps
    xstp(n)=xsize./(xnum(n)-1);
    ystp(n)=ysize./(ynum(n)-1);
    zstp(n)=zsize./(znum(n)-1);
end

% Coordinates
x=0:xstp(1):xsize;
y=0:ystp(1):ysize;
z=0:zstp(1):zsize;
xp=xstp(1)/2:xstp(1):xsize-xstp(1)/2;
yp=ystp(1)/2:ystp(1):ysize-ystp(1)/2;
zp=zstp(1)/2:zstp(1):zsize-zstp(1)/2;
xp1=-xstp(1)/2:xstp(1):xsize+xstp(1)/2;
yp1=-ystp(1)/2:ystp(1):ysize+ystp(1)/2;
zp1=-zstp(1)/2:zstp(1):zsize+zstp(1)/2;


% FINEST (PRINCIPAL) GRID
% Defining density structure rho()
% Defining initial guesses for velocity vx() vy() and pressure pr()
% Computing right part of Stokes (RX, RY) and Continuity (RC) equation
% vx, vy, P
vx1=zeros(ynum(1)+1,xnum(1),znum(1)+1);
vy1=zeros(ynum(1),xnum(1)+1,znum(1)+1);
vz1=zeros(ynum(1)+1,xnum(1)+1,znum(1));
pr1=zeros(ynum(1)-1,xnum(1)-1,znum(1)-1);
vxa=vx1;
vya=vy1;
vza=vz1;
pra=pr1;
% Right parts of equations
RX1=zeros(ynum(1)+1,xnum(1),znum(1)+1);
RY1=zeros(ynum(1),xnum(1)+1,znum(1)+1);
RZ1=zeros(ynum(1)+1,xnum(1)+1,znum(1));
RC1=zeros(ynum(1)-1,xnum(1)-1,znum(1)-1);
% Grid points cycle
rhoa=zeros(ynum(1),xnum(1),znum(1));
etaa=zeros(ynum(1),xnum(1),znum(1));
etan1=zeros(ynum(1),xnum(1),znum(1));
etaxy1=zeros(ynum(1),xnum(1),znum(1));
etaxz1=zeros(ynum(1),xnum(1),znum(1));
etayz1=zeros(ynum(1),xnum(1),znum(1));
for i=1:1:ynum(1)+1;
    for j=1:1:xnum(1)+1;
        for k=1:1:znum(1)+1;
            
            % Density Structure
            if (i<ynum(1)+1 && j<xnum(1)+1 && k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x(j),y(i),z(k),gx,gy,gz,gridtest);
                rhoa(i,j,k)=rho;
                etaa(i,j,k)=eta;
            end
            
            % Viscosity for normal stress (etan) for cells (in pressure nodes)
            if (i<ynum(1) && j<xnum(1) && k<znum(1))
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp(j),yp(i),zp(k),gx,gy,gz,gridtest);
                etan1(i,j,k)=eta;
                pra(i,j,k)=pr;
%                 if(i==1 || j==1 || k==1 || i==ynum(1)-1 || j==xnum(1)-1 || k==znum(1)-1)
%                     RC1(i,j,k)=pr;
%                 end
            end
            
            % Shear viscosity for SIGMAxy stress (etaxy) for z-directed edges
            if (i<ynum(1)+1 && j<xnum(1)+1 && k<znum(1))
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x(j),y(i),zp(k),gx,gy,gz,gridtest);
                etaxy1(i,j,k)=eta;
            end
            
            % Shear viscosity for SIGMAxz stress (etaxz) for y-directed edges
            if (i<ynum(1) && j<xnum(1)+1 && k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x(j),yp(i),z(k),gx,gy,gz,gridtest);
                etaxz1(i,j,k)=eta;
            end
            
            % Shear viscosity for SIGMAyz stress (etayz) for x-directed edges
            if (i<ynum(1)+1 && j<xnum(1)  && k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp(j),y(i),z(k),gx,gy,gz,gridtest);
                etayz1(i,j,k)=eta;
            end

           
            % Right part of x-Stokes Equation
            if(j>1 && i>1 && k>1 && i<ynum(1)+1 && j<xnum(1) && k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x(j),yp(i-1),zp(k-1),gx,gy,gz,gridtest);
                vxa(i,j,k)=vx;
                RX1(i,j,k)=-rho*gx;
            elseif(j<xnum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,x(j),yp1(i),zp1(k),gx,gy,gz,gridtest);
                vxa(i,j,k)=vx;
                RX1(i,j,k)=vx;
            end
            
            % Right part of y-Stokes Equation
            if(j>1 && k>1 && i>1 && i<ynum(1) && j<xnum(1)+1 && k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp(j-1),y(i),zp(k-1),gx,gy,gz,gridtest);
                vya(i,j,k)=vy;
                RY1(i,j,k)=-rho*gy;
            elseif(i<ynum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp1(j),y(i),zp1(k),gx,gy,gz,gridtest);
                vya(i,j,k)=vy;
                RY1(i,j,k)=vy;
            end
            
            % Right part of z-Stokes Equation
            if(j>1 && i>1 && k>1 && i<ynum(1)+1 && j<xnum(1)+1 && k<znum(1))
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp(j-1),yp(i-1),z(k),gx,gy,gz,gridtest);
                vza(i,j,k)=vz;
                RZ1(i,j,k)=-rho*gz;
            elseif(k<znum(1)+1)
                [eta,rho,vx,vy,vz,pr]=sol_visc3d(xsize,ysize,zsize,xp1(j),yp1(i),z(k),gx,gy,gz,gridtest);
                vza(i,j,k)=vz;
                RZ1(i,j,k)=vz;
            end
            
        end
    end
end
% Initial approximation
pr1=pra;
vx1=vxa;
vy1=vya;
vz1=vza;
% Pressure in the origin corner
prnorm=pra(1,1,1);

vxbeg=mean(mean(vxa(2:ynum(1),1,2:znum(1))))
vxend=mean(mean(vxa(2:ynum(1),xnum(1),2:znum(1))))
vybeg=mean(mean(vya(1,2:xnum(1),2:znum(1))))
vyend=mean(mean(vya(ynum(1),2:xnum(1),2:znum(1))))
vzbeg=mean(mean(vza(2:ynum(1),2:xnum(1),1)))
vzend=mean(mean(vza(2:ynum(1),2:xnum(1),znum(1))))
dvxdx=(vxend-vxbeg)/xsize
dvydy=(vyend-vybeg)/ysize
dvzdz=(vzend-vzbeg)/zsize
divv=dvxdx+dvydy+dvzdz



% Interpolation (restriction) of viscosity to coarser levels
for n=1:1:levelnum-1;
    % Action depends on the level of resolution
    switch n
            
        % Level 1 (principal grid)
        case 1
        if (levelnum>n)
            [etaxy2 etaxz2 etayz2 etan2]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy1,etaxz1,etayz1,etan1);
        end
        % Level 2 
        case 2
        if (levelnum>n)
            [etaxy3 etaxz3 etayz3 etan3]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy2,etaxz2,etayz2,etan2);
        end
        % Level 3 
        case 3
        if (levelnum>n)
            [etaxy4 etaxz4 etayz4 etan4]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy3,etaxz3,etayz3,etan3);
        end
        % Level 4 
        case 4
        if (levelnum>n)
            [etaxy5 etaxz5 etayz5 etan5]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy4,etaxz4,etayz4,etan4);
        end
        % Level 5 
        case 5
        if (levelnum>n)
            [etaxy6 etaxz6 etayz6 etan6]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy5,etaxz5,etayz5,etan5);
        end
        % Level 6 
        case 6
        if (levelnum>n)
            [etaxy7 etaxz7 etayz7 etan7]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy6,etaxz6,etayz6,etan6);
        end
        % Level 7 
        case 7
        if (levelnum>n)
            [etaxy8 etaxz8 etayz8 etan8]=Viscosity_restriction3D(n,xnum,ynum,znum,xstp,ystp,zstp,etaxy7,etaxz7,etayz7,etan7);
        end
        % No further restriction from last Level 8
    end
end




% Main Multigrid cycle
% We put 8 levels (for the case if we want to increase resolution)
% We use different arrays for different levels 
% since grid resolutions is different
% Multigrid schedule = V-cycle 
ynreset=0
for niter=1:1:inum;
    
    
    % Smoothing+restriction cycle
    for n=1:1:levelnum;
        % Action depends on the level of resolution
        switch n
            
            % Leve 1 (principal grid)
            case 1
            % Smoothing operation: solving of Stokes and Continuity equations on nodes
            % and computing residuals
            % by calling function Stokes_smoother()
            [vx1,resx1,vy1,resy1,vz1,resz1,pr1,resc1]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy1,etaxz1,etayz1,etan1,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX1,vx1,RY1,vy1,RZ1,vz1,RC1,pr1);
            % Restriction operation: 
            % Interpolating residuals to coarcer level (k+1) 
            % to produce right parts for this coarcer level
            if (levelnum>n)
                [RX2,RY2,RZ2,RC2]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx1,resy1,resz1,resc1,etan1,etan2);
            end
        
            % Level 2
            case 2
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx2=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy2=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz2=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr2=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx2,resx2,vy2,resy2,vz2,resz2,pr2,resc2]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy2,etaxz2,etayz2,etan2,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX2,vx2,RY2,vy2,RZ2,vz2,RC2,pr2);
            % Restriction operation: 
            if (levelnum>n)
                [RX3,RY3,RZ3,RC3]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx2,resy2,resz2,resc2,etan2,etan3);
            end
        
            % Level 3
            case 3
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx3=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy3=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz3=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr3=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx3,resx3,vy3,resy3,vz3,resz3,pr3,resc3]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy3,etaxz3,etayz3,etan3,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX3,vx3,RY3,vy3,RZ3,vz3,RC3,pr3);
            % Restriction operation: 
            if (levelnum>n)
                [RX4,RY4,RZ4,RC4]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx3,resy3,resz3,resc3,etan3,etan4);
            end
            
            % Level 4
            case 4
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx4=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy4=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz4=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr4=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx4,resx4,vy4,resy4,vz4,resz4,pr4,resc4]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy4,etaxz4,etayz4,etan4,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX4,vx4,RY4,vy4,RZ4,vz4,RC4,pr4);
            % Restriction operation: 
            if (levelnum>n)
                [RX5,RY5,RZ5,RC5]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx4,resy4,resz4,resc4,etan4,etan5);
            end
            
            % Level 5
            case 5
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx5=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy5=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz5=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr5=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx5,resx5,vy5,resy5,vz5,resz5,pr5,resc5]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy5,etaxz5,etayz5,etan5,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX5,vx5,RY5,vy5,RZ5,vz5,RC5,pr5);
            % Restriction operation: 
            if (levelnum>n)
                [RX6,RY6,RZ6,RC6]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx5,resy5,resz5,resc5,etan5,etan6);
            end
             
            % Level 6
            case 6
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx6=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy6=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz6=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr6=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx6,resx6,vy6,resy6,vz6,resz6,pr6,resc6]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy6,etaxz6,etayz6,etan6,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX6,vx6,RY6,vy6,RZ6,vz6,RC6,pr6);
            % Restriction operation: 
            if (levelnum>n)
                [RX7,RY7,RZ7,RC7]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx6,resy6,resz6,resc6,etan6,etan7);
            end
             
            % Level 7
            case 7
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx7=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy7=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz7=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr7=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx7,resx7,vy7,resy7,vz7,resz7,pr7,resc7]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy7,etaxz7,etayz7,etan7,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX7,vx7,RY7,vy7,RZ7,vz7,RC7,pr7);
            % Restriction operation: 
            if (levelnum>n)
                [RX8,RY8,RZ8,RC8]=Stokes_Continuity3D_viscous_restriction(n,xnum,ynum,znum,xstp,ystp,zstp,resx7,resy7,resz7,resc7,etan7,etan8);
            end
             
            % Level 8
            case 8
            %Making initial approximation for vx,vy,vz, P corrections (zeros)
            vx8=zeros(ynum(n)+1,xnum(n),znum(n)+1);
            vy8=zeros(ynum(n),xnum(n)+1,znum(n)+1);
            vz8=zeros(ynum(n)+1,xnum(n)+1,znum(n));
            pr8=zeros(ynum(n)-1,xnum(n)-1,znum(n)-1);
            % Smoothing operation: 
            [vx8,resx8,vy8,resy8,vz8,resz8,pr8,resc8]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy8,etaxz8,etayz8,etan8,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX8,vx8,RY8,vy8,RZ8,vz8,RC8,pr8);
            % No Restriction operation on this last level 

        end
    end
   
    


    
    
    % Smoothing+prolongation cycle
    for n=levelnum:-1:1;
        % Action depends on the level of resolution
        switch n

            % Leve 1 (principal grid)
            case 1
            % Smoothing operation: solving of Stokes and Continuity equations on nodes
            % and computing residuals
            % by calling function Stokes_smoother()
            [vx1,resx1,vy1,resy1,vz1,resz1,pr1,resc1]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy1,etaxz1,etayz1,etan1,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX1,vx1,RY1,vy1,RZ1,vz1,RC1,pr1);
            % No prolongation operation on finest level
        
            % Level 2
            case 2
            % Smoothing operation: 
            [vx2,resx2,vy2,resy2,vz2,resz2,pr2,resc2]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy2,etaxz2,etayz2,etan2,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX2,vx2,RY2,vy2,RZ2,vz2,RC2,pr2);
            % Prolongation operation: 
            % Interpolating dvx,dvy,dpr corrections to finer level (k+1) 
            % and computing solution update for this finer level
            [dvx1,dvy1,dvz1,dpr1]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx2,vy2,vz2,pr2);
            vx1=vx1+dvx1*vkoef;
            vy1=vy1+dvy1*vkoef;
            vz1=vz1+dvz1*vkoef;
            pr1=pr1+dpr1*pkoef;
        
            % Level 3
            case 3
            % Smoothing operation: 
            [vx3,resx3,vy3,resy3,vz3,resz3,pr3,resc3]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy3,etaxz3,etayz3,etan3,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX3,vx3,RY3,vy3,RZ3,vz3,RC3,pr3);
            % Prolongation operation: 
            [dvx2,dvy2,dvz2,dpr2]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx3,vy3,vz3,pr3);
            vx2=vx2+dvx2*vkoef;
            vy2=vy2+dvy2*vkoef;
            vz2=vz2+dvz2*vkoef;
            pr2=pr2+dpr2*pkoef;
            
            % Level 4
            case 4
            % Smoothing operation: 
            [vx4,resx4,vy4,resy4,vz4,resz4,pr4,resc4]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy4,etaxz4,etayz4,etan4,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX4,vx4,RY4,vy4,RZ4,vz4,RC4,pr4);
            % Prolongation operation: 
            [dvx3,dvy3,dvz3,dpr3]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx4,vy4,vz4,pr4);
            vx3=vx3+dvx3*vkoef;
            vy3=vy3+dvy3*vkoef;
            vz3=vz3+dvz3*vkoef;
            pr3=pr3+dpr3*pkoef;
            
             % Level 5
            case 5
            % Smoothing operation: 
            [vx5,resx5,vy5,resy5,vz5,resz5,pr5,resc5]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy5,etaxz5,etayz5,etan5,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX5,vx5,RY5,vy5,RZ5,vz5,RC5,pr5);
            % Prolongation operation: 
            [dvx4,dvy4,dvz4,dpr4]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx5,vy5,vz5,pr5);
            vx4=vx4+dvx4*vkoef;
            vy4=vy4+dvy4*vkoef;
            vz4=vz4+dvz4*vkoef;
            pr4=pr4+dpr4*pkoef;
            
            
            % Level 6
            case 6
            % Smoothing operation: 
            [vx6,resx6,vy6,resy6,vz6,resz6,pr6,resc6]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy6,etaxz6,etayz6,etan6,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX6,vx6,RY6,vy6,RZ6,vz6,RC6,pr6);
            % Prolongation operation: 
            [dvx5,dvy5,dvz5,dpr5]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx6,vy6,vz6,pr6);
            vx5=vx5+dvx5*vkoef;
            vy5=vy5+dvy5*vkoef;
            vz5=vz5+dvz5*vkoef;
            pr5=pr5+dpr5*pkoef;
            
            % Level 7
            case 7
            % Smoothing operation: 
            [vx7,resx7,vy7,resy7,vz7,resz7,pr7,resc7]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy7,etaxz7,etayz7,etan7,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX7,vx7,RY7,vy7,RZ7,vz7,RC7,pr7);
            % Prolongation operation: 
            [dvx6,dvy6,dvz6,dpr6]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx7,vy7,vz7,pr7);
            vx6=vx6+dvx6*vkoef;
            vy6=vy6+dvy6*vkoef;
            vz6=vz6+dvz6*vkoef;
            pr6=pr6+dpr6*pkoef;
            
            % Level 8
            case 8
            % Smoothing operation: 
            [vx8,resx8,vy8,resy8,vz8,resz8,pr8,resc8]=Stokes_Continuity3D_viscous_smoother(prnorm,etaxy8,etaxz8,etayz8,etan8,iternum(n),krelaxs,krelaxc,xnum(n),ynum(n),znum(n),xstp(n),ystp(n),zstp(n),RX8,vx8,RY8,vy8,RZ8,vz8,RC8,pr8);
            % Prolongation operation: 
            [dvx7,dvy7,dvz7,dpr7]=Stokes_Continuity3D_prolongation(n,xnum,ynum,znum,xstp,ystp,zstp,vx8,vy8,vz8,pr8);
            vx7=vx7+dvx7*vkoef;
            vy7=vy7+dvy7*vkoef;
            vz7=vz7+dvz7*vkoef;
            pr7=pr7+dpr7*pkoef;
        end
    end 
 
    % Defining scale for Stokes residuals from y-Stokes equation
    % dSIGMAij/dj-dP/di=-RHO*gi=0  => Stokes scale=abs(RHO*gi)
    stokesscale= rhoblock*g;
    % Defining scale for Continuity residuals from y-Stokes equation
    % dvx/dx+dvy/dy=0 can be transformed to 2ETA(dvx/dx+dvy/dy)/dx=0 
    % which is similar to dSIGMAij/dj and has scale given above
    % therefore continuity scale = scale=abs(RHO*gi/ETA*dx)
    continscale= rhoblock*g/etamedium*xstp(1);


    % Errors
    vxe=vx1-vxa;
    vye=vy1-vya;
    vze=vz1-vza;
    pre=pr1-pra;
    
    
     % Plotting Properties
    figure(4);clf
    % Etaa
    subplot(1,2,1)
    pcolor(log10(etaa(:,:,(znum(1)-1)/2+1)));
    shading interp; colorbar
    title('Etaa at z=0.5')
    % Rhoaa
    subplot(1,2,2)
    pcolor(log10(rhoa(:,:,(znum(1)-1)/2+1)));
    shading interp; colorbar
    title('Rhoa at z=0.5')
    
    
    
   % Plotting Results
    figure(1);clf
    % Numerical
    % Vxn
    subplot(3,4,1)
    pcolor(vx1(2:ynum(1),2:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vxn at z=0.5')
    % Vyn
    subplot(3,4,2)
    pcolor(vy1(2:ynum(1)-1,2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vyn at z=0.5')
    % Vzn
    subplot(3,4,3)
    pcolor(vz1(2:ynum(1),2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vzn at z=0.5')
    % Prn
    subplot(3,4,4)
    pcolor(pr1(1:ynum(1)-1,1:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Prn at z=0.5')
    
    % Analytical
    % Vxa
    subplot(3,4,5)
    pcolor(vxa(2:ynum(1),2:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vxa at z=0.5')
    % Vya
    subplot(3,4,6)
    pcolor(vya(2:ynum(1)-1,2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vya at z=0.5')
    % Vza
    subplot(3,4,7)
    pcolor(vza(2:ynum(1),2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vza at z=0.5')
    % Pra
    subplot(3,4,8)
    pcolor(pra(1:ynum(1)-1,1:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Pra at z=0.5')
    
    % Error
    % Vxe
    subplot(3,4,9)
    pcolor(vxe(2:ynum(1),2:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vxe at z=0.5')
    % Vye
    subplot(3,4,10)
    pcolor(vye(2:ynum(1)-1,2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vye at z=0.5')
    % Vze
    subplot(3,4,11)
    pcolor(vze(2:ynum(1),2:xnum(1),(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Vze at z=0.5')
    % Pre
    subplot(3,4,12)
    pcolor(pre(1:ynum(1)-1,1:xnum(1)-1,(znum(1)-1)/2+1));
    shading interp; colorbar
    title('Pre at z=0.5')
    
    
    
    %Residuals 
    resx01=resx1/stokesscale;
    resy01=resy1/stokesscale;
    resz01=resz1/stokesscale;
    resc01=resc1/continscale;
    
    % Computing mean square residuals
    resx001(niter)=0;
    resy001(niter)=0;
    resz001(niter)=0;
    resc001(niter)=0;
    for i=1:1:ynum(1)
        for j=1:1:xnum(1)
            for k=1:1:znum(1)
                % x-Stokes
                if (i>1 && j>1 && k>1 && j<xnum(1))
                    resx001(niter)=resx001(niter)+resx01(i,j,k)^2;
                end
                % y-Stokes
                if (i>1 && j>1 && k>1 && i<ynum(1))
                    resy001(niter)=resy001(niter)+resy01(i,j,k)^2;
                end
                % z-Stokes
                if (i>1 && j>1 && k>1 && k<znum(1))
                    resz001(niter)=resz001(niter)+resz01(i,j,k)^2;
                end
                % Continuity
                if (j<xnum(1) && i<ynum(1) && k<znum(1))
                    resc001(niter)=resc001(niter)+resc01(i,j,k)^2;
                end
            end
        end
    end
    resx001(niter)=log10((resx001(niter)/((ynum(1)-1)*(xnum(1)-2)*(znum(1)-1)))^0.5);
    resy001(niter)=log10((resy001(niter)/((ynum(1)-2)*(xnum(1)-1)*(znum(1)-1)))^0.5);
    resz001(niter)=log10((resz001(niter)/((ynum(1)-1)*(xnum(1)-1)*(znum(1)-2)))^0.5);
    resc001(niter)=log10((resc001(niter)/((ynum(1)-1)*(xnum(1)-1)*(znum(1)-1)))^0.5);
   

    % Plotting Mean residuals
    figure(2);clf
    plot(resc001, 'k');
    hold on
    plot(resx001, 'b');
    hold on
    plot(resy001, 'g');
    hold on
    plot(resz001, 'r');
    xlabel('V-cycles')
    ylabel('log(Residuals)')

    zlabel(['RMS errors, cycle ',num2str(niter)])
%     title(['RMS errors, Multigrid V-cycle = ',num2str(niter)])
%     colorbar;

 

    pause(.001)
    
    
if(resc001(niter)<-12) 
    break
end

    
end

% ERROR NORMS
% P
pravr=mean(mean(mean(abs(pra(1:ynum(1)-1,1:xnum(1)-1,1:znum(1)-1)))));
vxavr=mean(mean(mean(abs(vxa(2:ynum(1),2:xnum(1)-1,2:znum(1))))));
vyavr=mean(mean(mean(abs(vya(2:ynum(1)-1,2:xnum(1),2:znum(1))))));
vzavr=mean(mean(mean(abs(vza(2:ynum(1),2:xnum(1),2:znum(1)-1)))));
vavr=(vxavr^2+vyavr^2+vzavr^2)^0.5;
hres(r)=xstp(1);
L12inf(r,1)=mean(mean(mean(abs(pre(1:ynum(1)-1,1:xnum(1)-1,1:znum(1)-1)))))/pravr; % L1
L12inf(r,2)=mean(mean(mean(pre(1:ynum(1)-1,1:xnum(1)-1,1:znum(1)-1).^2))).^0.5/pravr; % L2
L12inf(r,3)=max(max(max(abs(pre(1:ynum(1)-1,1:xnum(1)-1,1:znum(1)-1)))))/pravr; % Linf
% vx
L12inf(r,4)=mean(mean(mean(abs(vxe(2:ynum(1),2:xnum(1)-1,2:znum(1))))))/vavr; % L1
L12inf(r,5)=mean(mean(mean(vxe(2:ynum(1),2:xnum(1)-1,2:znum(1)).^2))).^0.5/vavr; % L2
L12inf(r,6)=max(max(max(abs(vxe(2:ynum(1),2:xnum(1)-1,2:znum(1))))))/vavr; % Linf
% vy
L12inf(r,7)=mean(mean(mean(abs(vye(2:ynum(1)-1,2:xnum(1),2:znum(1))))))/vavr; % L1
L12inf(r,8)=mean(mean(mean(vye(2:ynum(1)-1,2:xnum(1),2:znum(1)).^2))).^0.5/vavr; % L2
L12inf(r,9)=max(max(max(abs(vye(2:ynum(1)-1,2:xnum(1),2:znum(1))))))/vavr; % Linf
% vz
L12inf(r,10)=mean(mean(mean(abs(vze(2:ynum(1),2:xnum(1),2:znum(1)-1)))))/vavr; % L1
L12inf(r,11)=mean(mean(mean(vze(2:ynum(1),2:xnum(1),2:znum(1)-1).^2))).^0.5/vavr; % L2
L12inf(r,12)=max(max(max(abs(vze(2:ynum(1),2:xnum(1),2:znum(1)-1)))))/vavr; % Linf

figure(3); clf
hold on
plot(log10(hres),log10(L12inf(:,1)),'b -');
plot(log10(hres),log10(L12inf(:,2)),'b .-');
plot(log10(hres),log10(L12inf(:,3)),'b --');
plot(log10(hres),log10(L12inf(:,4)),'r -');
plot(log10(hres),log10(L12inf(:,5)),'r .-');
plot(log10(hres),log10(L12inf(:,6)),'r --');
plot(log10(hres),log10(L12inf(:,7)),'k -');
plot(log10(hres),log10(L12inf(:,8)),'k .-');
plot(log10(hres),log10(L12inf(:,9)),'k --');
plot(log10(hres),log10(L12inf(:,10)),'g -');
plot(log10(hres),log10(L12inf(:,11)),'g .-');
plot(log10(hres),log10(L12inf(:,12)),'g --');
hold off
title('error norms convergence');
pause(5)
end


