femwind/multigrid_solve.m

   1 function [x,it,rate,X_coarse]=multigrid_solve(K,F,X,params)
   2 % x=multigrid_solve(K,F,X)
   3
   4 level=params.level;  % debug
   5 n = size(X{1});
   6 fprintf('multigrid level %i grid size %i %i %i rhs %g\n',params.levels,n,big(F))
   7
   8 disp(['coarsening method ',params.coarsening])
   9 disp(['smoothing method ',params.smoothing])
  10
  11 nn = size(F,1);
  12 if nn~=prod(n)
  13     error('rb_line_gs_solve: inconsistent sizes')
  14 end
  15 rate = 0;
  16
  17 % settings
  18 tol = params.restol*norm(F);
  19 if params.debug_level >= params.level
  20     stop_here=true
  21 end
  22
  23
  24 % constant arrays
  25 x = zeros(nn,1);
  26 colz=1:n(3);
  27 onez=ones(1,n(3));
  28
  29 disp('building prolongation')
  30 switch params.coarsening
  31     case 'vertical lines'
  32         nnc = n(1)*n(2);  % number of coarse
  33         P = spalloc(nn,nnc,nn);
  34         for i1=1:n(1)
  35             for i2=1:n(2)
  36                 colw=X{3}(i1,i2,end)-X{3}(i1,i2,:);
  37                 ix = sub2ind(n,i1*onez,i2*onez,colz);
  38                 ixc = sub2ind(n(1:2),i1,i2);
  39                 P(ix,ixc)=squeeze(colw);
  40             end
  41         end
  42     case '2 linear'
  43         % [P,X_coarse]=coarsening_2_linear(X,params);
  44         dx=(X{1}(2,1,1)-X{1}(1,1,1));
  45         dy=(X{2}(1,2,1)-X{1}(1,1,1));
  46         dz = squeeze(X{3}(1,1,2:end)-X{3}(1,1,1:end-1)); % ref z spacing
  47         fprintf('layers: '),disp(X{3}(1,1,:)),fprintf('dz='),disp(dz)
  48         [hzc,icl3]=coarsening_icl_fortran(dx,dy,dz,params);
  49         X_coarse=coarsening_X(hzc,icl3,X,params);
  50         nnc=prod(size(X_coarse{1}));
  51         % P=coarsening_P(icl,X,params);
  52         prol_rest_err(hzc,icl3,X,params);
  53     otherwise
  54         error(['unknown coarsening ',params.coarsening])
  55 end
  56 fprintf('coarsening level variables %i coarse %i ratio %g\n',params.levels,nn,nnc,nn/nnc);
  57 disp('computing coarse matrix')
  58 diary;diary
  59 switch params.coarse_K
  60     case {1,'variational'}
  61         disp('coarse K is variational')
  62         P=coarsening_P(icl,X,params);
  63         K_coarse = P'*K*P;
  64         check_nonzeros(params.levels-1,K_coarse,X_coarse,P,K,X);
  65     case {2,'assembly'}
  66         disp('coarse K by assembly')
  67         % assemble sparse system matrix
  68         [K_coarse,~,~] = sparse_assembly(diag(params.a),X_coarse,[],[],params);
  69         % dirichlet boundary conditions
  70         [K_coarse,~]=apply_boundary_conditions(K_coarse,[],X_coarse);
  71     otherwise
  72         disp(params.coarse_P)
  73         error('unknown coarse_P method')
  74 end
  75 if params.exact
  76     disp('exact solution, for comparison only')
  77     ex = K\F;
  78 end
  79 cycles=0;
  80 t_cycle='first cycle not complete yet';
  81 if params.levels<=1 || nn == nnc % coarsest level
  82     if params.coarsest_iter==0   % direct
  83         fprintf('multigrid coarsest level %i solving directly\n',params.level)
  84         x = K\F;
  85     else
  86         x =zeros(size(F));       % iterative
  87         fprintf('multigrid coarsest level %i solving by %i iterations\n',params.level,params.coarsest_iter)
  88         for it=1:params.coarsest_iter
  89             % fprintf('level %g iteration %g coarsest solve\n',params.level,it)
  90             xout=smoothing(K,F,X,x,params);
  91             x = xout;
  92         end
  93         res=big(K*x-F);relres=res/big(F);rate=relres^(1/params.coarsest_iter);
  94         fprintf('multigrid coarsest residual %g relative %g rate %g\n',res,relres,rate);
  95     end
  96     return
  97 end
  98 for it=1:params.maxit
  99     diary;diary
 100     params.it(params.levels)=it;  % where we are, for prints and filenames
 101     coarse = mod(it,params.nsmooth+1)==0;
 102     if coarse
 103         fprintf('level %g iteration %g coarse correction\n',params.level,it)
 104         x=coarse_correction(x,F,K,K_coarse,X_coarse,hzc,icl3,X,params);
 105         it_type=sprintf('level %g coarse correction',params.levels);
 106     else
 107         it_type='smoothing';
 108         fprintf('level %g iteration %g smoothing by %s\n',params.level,it,params.smoothing)
 109         xout=smoothing(K,F,X,x,params);
 110         x=xout;
 111     end
 112     r=F-K*x;  % residual for diagnostics only
 113     if params.exact
 114         e=x-ex;   % error
 115     else
 116         e=[];
 117     end
 118     res(it)= norm(r);
 119     if mod(it,params.nsmooth+1)==params.nsmooth
 120         cycles=cycles+1;
 121         rate = (res(it)/norm(F))^(1/cycles);
 122         t_cycle=sprintf('cycle %g level %g avg rate %g',cycles,params.levels,rate);
 123         disp(t_cycle)
 124     end
 125     % if params.debug_level >= params.level
 126     if  params.level == 1 && it == 20
 127         stop_here=true
 128     end
 129     if params.graphics > -1
 130         tstring=sprintf('it=%g %s %s %s',it,it_type,t_cycle);
 131         plot_error_slice(e,r,X,tstring,params)
 132     end
 133     if params.graphics > -2
 134         figure(params.iterations_fig)
 135         semilogy(1:it,res/norm(F),'*')
 136         legend('relative 2-norm residual')
 137         grid on
 138         title([sprintf('mesh=%g %g %g ',n),t_cycle])
 139         xlabel('iteration')
 140         drawnow,pause(0.1)
 141     end
 142     if params.exact
 143         err(it)=norm(e); % l2 error
 144         eer(it)=sqrt(e'*K*e); % energy norm error
 145         hold on, semilogy(1:it,err/err(1),'x',1:it,eer/eer(1),'+'), hold off
 146         legend('relative 2-norm residual','relative 2-norm error','relative energy norm error')
 147     end
 148     fprintf('iteration %i residual %g tolerance %g\n',it,res(it),tol)
 149     if res(it)<tol,
 150         fprintf('exiting level %i in iteration %i because residual %g < tolerance %g\n',...
 151             params.level,it,res(it),tol)
 152         if params.debug_level >= 0
 153             stop_here=true;
 154         end
 155         break
 156     end
 157 end
 158 if params.save_files > 2
 159     sfile=sprintf('%s_%s_%i.mat',params.save_file_prefix,params.id,params.levels);
 160     save(sfile,'-v7.3')
 161 end
 162 end
 163
 164