function [k,u] = nma_solver_Vcycle(u,f,mu1,mu2,METHOD,tol,PLOT_OPTION)
% Solve poisson 2D pde on unit square zero BC using multigrid V cycle method

%------------------------
% Function to solve poisson 2D pde on unit square with zero boundary
% conditions using multigrid V cycle method. This function sets up cycle,
% then calls nma_V_cycle() to actually perform the cycle computation in
% a loop until convegence. file name: nma_solver_Vcycle.m
% INPUT:
%    u: 2D grid contains the initial guess of solution
%    f: 2D grid contains the force field
%    mu1: number of pre-smoothers
%    mu2: number of post-smoothers
%    METHOD: the smoother to use, supported ones are
%         [1] Gauss-Seidel red/black method, [2] means use Jacobi
%         [3] means use Gauss-Seidel LEX     [4] means use SOR
%    tol: numerical value for tolerance. Such as 0.01*h^2 where h is
%         the mesh spacing.
%    PLOT_OPTION: flag, either true or false. If true, then the solution
%         is plotted after each iteration on the screen.
%  OUTPUT:
%    k: number of iterations to reach convergence
%    u: the solution reached, 2D grid
%
% Nasser M. Abbasi
% Math 228a, UC Davis, Fall 2010


[n,~] = size(f);
h     = 1/(n-1);
[X,Y] = meshgrid(0:h:1,0:h:1); % coordinates
done  = false;
k     = 0;  % loop counter

normf = nma_find_norm(f);  % find ||f|| once.

while not(done)   %keep looping until convergence, multigrid solver
   
   k = k+1;
   u = nma_V_cycle(u,f,mu1,mu2,METHOD);  %--- CALL V CYCLE
   the_residue      = nma_find_residue(u,f);
   relative_residue = nma_find_norm(the_residue)/normf;
   
   if relative_residue < tol
      done = true;
   end
   
   if PLOT_OPTION
      mesh(X,Y, u);
      drawnow; title(sprintf('k=%d, tol=%f, n=%d, h=%f\n',k,tol,n,h));
   end
end
end