%Who  : Abdallah Sayyed Ahmad.
%What : 2D Poission equation solver;DelV=f and bc1 at the bottom, bc2 at the right, bc3 at top and bc4 at left of a rectangular box
%Where: Computational Physics 338/Department of Physics/Birzeit University
%When : December 2020
clear;
clc
solver=2;
f = inline('(6*x*y*(1-y)-2*x^3)', 'x', 'y');
f = inline('0.0', 'x', 'y');
bc1 = inline('0.0', 'x');
bc2 = inline('y*(1-y)', 'y');
bc3 = inline('0.0', 'x');
bc4 = inline('0.0', 'y');

N=100;
x0=0;
x1=1;
y0=x0;
y1=x1;
dx=(x1-x0)/(N-1);
dy=dx;
NN=N-1;

% NU is the number of unknows.
NU=(N-2)*(N-2);
Row=zeros(5*NU,1);
Col=zeros(5*NU,1);
Val=zeros(5*NU,1);
b=zeros(NU,1);
unknown=1;
k=1;
row=0;
for j=2:(N-1)
  yy=y0+(j-1)*dy;      
    for i=2:(N-1)
      xx=x0+(i-1)*dx;      
        row=row+1;       
        Row(k)=row;
        Col(k)=row;
        Val(k)=-4;
        k=k+1;
        %right hand side
        b(unknown)=dx*dx*f(xx,yy);
        
        %left part.
        if(i>2)
            Row(k)=row;
            Col(k)=row-1;
            Val(k)=1;
            k=k+1;
        else
            %left boundary.
            b(unknown)=b(unknown)-(bc4(yy));
        end
        %right part.
        if(i<(N-2))
            Row(k)=row;
            Col(k)=row+1;
            Val(k)=1;
            k=k+1;
        else
            %right boundary.
            b(unknown)=b(unknown)-(bc2(yy));
        end
              
        
        %lower part.
        if(j>2)
            Row(k)=row;
            Col(k)=row-N+2;
            Val(k)=1;
            k=k+1;
        else
            %bottom boundary.
            b(unknown)=b(unknown)-(bc1(xx));
        end
        %upper part.
        if(j<(N-2))
            Row(k)=row;
            Col(k)=row+N-2;
            Val(k)=1;
            k=k+1;
        else
            %upper boundary.
            b(unknown)=  b(unknown)-(bc3(xx));
        end
        
        unknown=unknown+1;
        
    end
end
Row=Row(1:k-1);
Col=Col(1:k-1);
Val=Val(1:k-1);
A=sparse(Row,Col,Val);

tic
tol=1e-12;
maxits =600;
switch solver
case 1
    %matlab build in direct solver.
    sol=A\b;
case 2
    %Biconjugate gradient solver without preconditioer.
    sol= bicgstab(A,b,tol,maxits); 
 
case 3
    %calculate the incomplete LU factorization of A. Given droptol.
    opts.type = "crout";
    opts.droptol = 1e-2;
    [L,U] = ilu(A,opts); 
    %biconjugate gradient solver with the incomplete LU factors as a preconditioner.
   sol= bicgstab(A,b,tol,maxits,L,U); 
 end
toc


x = linspace(x0,x1,N);
y = linspace(y0,y1,N);
[xx,yy] = meshgrid(x,y);
ZZ = reshape(sol(:), N-2, []).';

Z=zeros(N);
Z(2:N-1,2:N-1)=ZZ;
vectorize(bc1);
for i=1:N
 x=x0+(i-1)*dx; 
Z(i,1)=bc1(x);
Z(i,N)=bc3(x);
end

for j=1:N
y=y0+(j-1)*dy; 
Z(1,j)=bc2(y);
Z(N,j)=bc4(y);
end

%draw solution
contourf(xx,yy,Z)
% compute exact solution; this need to be changed when the source term or the boundary conditions changed
exact=yy.*(1-yy).*xx.^3;
% compute error in solution per grid point
error=norm(Z-exact)/N^2