%
% Solves Laplace's equation on a square region
% (see page 167 of Tannehill, Anderson, & Pletcher)
% Using Multi-grid & Gauss-Seidel
%
% The problem domain looks like :
%
%       
%         Y   / \
%              |
%              |
%              |          u=100
%              |---------------------------|
%              |                           |
%              |                           |
%              |                           |
%              |                           |
%              |                           |
%              |        __2                |
%        u=50  |        \/ u = 0           |  u=150
%              |                           |
%              |                           |
%              |                           |
%              |                           |
%              |                           |
%              |                           |
%              |---------------------------|---> X
%                         u=200


% The process is :
%
%    1. Do 3-4 iterations on the fine grid
%    2. Interpolate residual onto coarse grid (restriction operation)
%    3. Do iterations on coarse mesh
%    4. Interpolate corrections onto fine grid (prolongation operation)
%
    clear;
    echo off;
    iter = 100;

%   Choosing this will use simply Gauss-Seidel:
%      nfine   =  1;
%   Choosing this will use multigrid:
      nfine   =  4;


    res=zeros(iter,1);

    print=iter/2;

% fine grid:    
    nx = 25;

    r=zeros(nx,nx);
    du=zeros(nx,nx);
    u=zeros(nx,nx);

% coarse grid:
    nxc = (nx+1)/2;
    
    rc =zeros(nxc,nxc);
    duc=zeros(nxc,nxc);
    ducold=zeros(nxc,nxc);
    uc =zeros(nxc,nxc);


% fine grid
    x=zeros(nx,nx);
    y=zeros(nx,nx);
    for i=1:nx
      for j=1:nx
         x(i,j) = i;
         y(i,j) = j;
      end
    end

    dx = 1.0;
    dxc = 2.0 * dx ;



%  ===== fine grid bc's =====

% u=50 at x=0

    for j=1:nx
      u(1,j) = 50.0;
    end

% u=150 at x=xmax    
    for j=1:nx
      u(nx,j) = 150.0;
    end
    
% u=200 at y=0
    for i=1:nx
      u(i,1) = 200.0;
    end
    
% u=100 at y=ymax
    for i=1:nx
      u(i,nx) = 100.0;
    end


% ===== coarse grid bc's =====

% u=50 at x=0

    for j=1:nxc
      uc(1,j) = 50.0;
    end

% u=150 at x=xmax    
    for j=1:nxc
      uc(nxc,j) = 150.0;
    end
    
% u=200 at y=0
    for j=1:nxc
      uc(i,1) = 200.0;
    end
    
% u=100 at y=ymax
    for i=1:nxc
      uc(i,nxc) = 100.0;
    end


  for n=1:iter

%-------------
%  fine grid
%-------------
    for g1=1:nfine

%
% update solution
    for i=2:nx-1
      for j=2:nx-1
         u(i,j) = (u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))*0.25;
      end
    end


      end


% find residual on fine grid
    for i=2:nx-1
      for j=2:nx-1
         r(i,j) = (-4.0*u(i,j)+u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));
      end
    end

%------------
%  Coarse grid
%------------

 if  nfine==1, 

% inject fine residual onto coarse residual
    for i=2:nxc-1
      for j=2:nxc-1
         rc(i,j) = r( (2*i-1), (2*j-1) );
      end
    end

    res3=1.0;
    while res3>0.005,
      
      ducold = duc;

%      update correction
      for i=2:nxc-1
        for j=2:nxc-1
           duc(i,j) = (duc(i+1,j)+duc(i-1,j)+duc(i,j+1)+duc(i,j-1)+...
                         dxc^2 * rc(i,j)) * 0.25;
        end
      end

      res4 = duc - ducold;
      res3 = sqrt(sum(sum(res4.*res4 )/nx/nx));       
      
    end


% copy correction to fine grid

    for i=2:nxc-1
      for j=2:nxc-1
        du( (2*i-1), (2*j-1) ) =  du(i,j) ;
      end
    end

% interpolate du to other fine grid points   

     for i=2:nx:2
       for j=3:nx:2
          du(i,j) = 0.5 * ( du(i+1,j) + du(i-1,j) );
       end
      end

     for i=3:nx:2
       for j=2:nx:2
          du(i,j) = 0.5 * ( du(i+1,j) + du(i-1,j) );
       end
      end

     for i=2:nx:2
       for j=2:nx:2
          du(i,j) = 0.25 * ( du(i+1,j+1) + du(i-1,j-1) +...
                             du(i-1,j+1) + du(i+1,j-1) );
       end
      end

% find new u on fine grid
     for i=2:nx-1
       for j=2:nx-1
          u(i,j) = u(i,j) + du(i,j);
       end
      end

    end  

    n,res2 = sqrt(sum(sum(r.*r )/nx/nx)) 
    res(n) = log10( res2 );
    
     if ( rem(n,print)==0 )


% this is just for plotting purposes:
       u(nx,nx) = (u(nx-1,nx  ) + u(nx,nx-1))*0.5;
       u( 1,nx) = (u(   1,nx-1) + u( 2,nx  ))*0.5;
       u(nx, 1) = (u(nx-1,1   ) + u(nx, 2  ))*0.5;
       u( 1, 1) = (u(   1,2   ) + u( 2, 1  ))*0.5;

  %  pcolor(x,y,u)
  %   contour(u)   

%    subplot(2,1,1), contour(u)
    subplot(2,1,1), surf(x,y,u)
    xlabel('X')
    ylabel('Y')
    zlabel('PHI')
    if  nfine==1,
       title('Gauss-Seidel Applied to Laplace Equation')
    else
       title('Multigrid (& Gauss-Seidel) Applied to Laplace Equation')
     end

     pause(1)

    subplot(2,1,2),     plot(res)
    xlabel('Iteration')
    ylabel(' Log10 ( Residual) ')

    end


   end
     

% the exact solution should have u equal to the 
% average of the wall values (at steady state):

    error = 0.25*(50+100+150+200)- u( (nx+1)/2, (nx+1)/2 )