quicwind/multigrid_2d_test.m

   1 function relres=multigrid_2d_test
   2 N=10;
   3 nc=[4,3];   % the coarsest grid points
   4 a=[1,3];    % domain size
   5 % solve Dirichlet problem with vertex nodes
   6 % homogeneous boundary conditions
   7 % mesh is 0:n1+1  times 0:n2+1 but the boundary is not stored
   8 % refinement:
   9 %    0       1     ...      n1          n1+1        level l+1
  10 %    0   1   2     ...     2*n1  2*n1+1 2*(n1+1)    level l
  11
  12 for l=N:-1:1
  13     n{l} = (nc+1)*2^(N-l)-1;                % mesh points
  14     h{l}=[a./(n{l}+1)];                     % mesh step
  15     A{l}=@(x)flat(@mlap,x,n{l},h{l});       % matrix vector multiply
  16     M{l}=@(x)x;                             % identity preconditioner
  17     if l<N
  18           P{l}=@(x)flat(@prolongation_2d,x,n{l+1});
  19           R{l}=@(x)flat(@restriction_2d,x,n{l});
  20     end
  21     fprintf('level %i size %i %i\n',l,n{l})
  22 end
  23 f=zeros(n{1});
  24 c=round(n{1}/2);
  25 f(c(1),c(2))=1;
  26 % set params
  27 p.tol=1e-9;       % relative residual tolerance
  28 p.maxit=7;        % maximum number of multigrid iterations
  29 p.tolsm1=p.tol;    % tolerance in pre-smoothing,
  30 p.tolsm2=p.tol;    % tolerance in pre-smoothing
  31 p.tolcr=p.tol;     % tolerance for coarse solver
  32 p.maxsm1=7;       % max pre-smoothing iterations
  33 p.maxsm2=7;       % max post-smoothing iterations
  34 p.maxcr=20;       % max coarse solve iterations
  35 p.ncoarse=1;      % number of coarse solves 1=V-cycle, 2=W-cycle
  36
  37 [x,relres]=multigrid(A,f(:),[],M,P,R,p);
  38 x = reshape(x,n{1});
  39 mesh(x)
  40
  41 end