1.m

五点差分型多重网格方法：各种插值算子的比较）

%==============================================
%      MAIN SCRIPT OF MULTIGRID   
%==============================================
%  zhiyong2008 finish
%  solve poisson etc. equation using DIFF.
%  input can be replaced and you can do many works ...
%  ==  output pde_error(L∞) in increasing level :7-15-31-63-127-...
%  our mesh ordered by from south to north ,from left to right
%                 ---③--⑥--⑨--
%                    |   |   |
%                 ---②--⑤--⑧--
%                    |   |   |
%                  --①--④--⑦--
%  main elements:
%                 linear systems
%                 mg algrithms: highlevel-->  midlevel-->  lowlevel
%                 solver
%  [e-mail]

function demo3_run

demo_globals
include_globals
include_flags
DEMO_VAR1 == 3; 

if DEMO_VAR1 == 1  
   NX1 = [7 15];
elseif DEMO_VAR1 == 2
   NX1 = [7 15 31];
elseif DEMO_VAR1 == 3
   NX1 = [7 15 31 63];
else
   disp('*****Starting my Multigrid:Zhiyong 2008 work*****')
end
for nx1 = [  63  ];
    ny1 =nx1;                    % Set fine mesh size


    % Problem parameter defaults

    prob_args = [0 0];

    % Mesh parameter defaults

    coarse_level = 3;            % Set current maximum number of levels

    % Smoother parameter defaults

    nu1 = 2; nu2 = 1;            % Set number of pre- and post-smoothings
    wt = 0.95;                   % Set smoother weighting

    % Selection flag defaults

    problem_flag = POISSON;      % Select continuous problem
    rhs_flag = 1;                % Select right-hand side
    %------------------------------------------------------------------
   %smooth_flag = WEIGHTED_JACOBI;% Select smoother
    smooth_flag = GAUSS_SEIDEL;  
   % smooth_flag = RB_GAUSS_SEIDEL;
    %----------------------------------------------------------------
  %  restrict_flag = INJECTION; % Select restriction operator
    restrict_flag = FULL_WEIGHTING;
  %  restrict_flag = HALF_WEIGHTING;
    
    %restrict_flag = ZHIJIE; %This is INJECTION designed in my work.
   % restrict_flag = IVER_FULL_LOCAL_EQUATION;
    %restrict_flag = IVER_SEMI_LOCAL_EQUATION;
    %restrict_flag = RES_LOCAL;
    %restrict_flag = RES_LOCAL_SECOND;
    %----------------------------------------------------------------
     interp_flag = LINEAR;          % Select interpolation operator
   % interp_flag = FULL_LOCAL_EQUATION;
    
    %interp_flag = SEMI_LOCAL_EQUATION;
    %interp_flag = EIGHT_POINTS;
    %interp_flag = ONE_DIM_LOCAL_EQUATION;
    %interp_flag = EAST_SOUTH_AVERAGE;
    %interp_flag = WEST_SOUTH_AVERAGE;
    %interp_flag = TEST_FIVEPOINT ;
    %---------------------------------------------------------------
    solver_flag = VMG_SOLVE;           % Select solver
    precon_flag = NONE;          % Select preconditioner
    cycle_flag = V_CYCLE;        % Select cycle flag
    coarsening_flag = 0;         % Select coarsening
    coarse_solver_flag = DIRECT_SOLVE; % Select coarse-grid solve 
    cycle_flag = V_CYCLE;        % Select MG cycle

    % Method parameter defaults

    rtol = 1e-16;           % Set stopping tolerance on weighted residual
    prtol = 1e-16;          % Set stopping tol. on weighted pseudo-resid.
    max_it = 12;                 % Set maximum number of iterations
    max_time = 0;                % Set maximum number of seconds
    max_mflop = 0;               % Set maximum number of mflops
    num_runs = 1;                % Set number of experiments to run

    % Linear System defaults

    generate_matrix = 1;         % Must generate the matrix initially
    generate_rhs = 1;        % Must generate the right-hand side initially
    matrix_type = [];            % Matrix properties initially unknown


    if (generate_matrix)
       [A1,N1] = get_matrix(nx1,ny1);
       generate_matrix = 0;
       multigrid_setup;
    end
    
    if (generate_rhs) 
       h = 1 / (nx1+1);
       [XX,YY] = meshgrid([h:h:(1-h)]);
       %UTRUE = sin(2*pi*XX) .* sin(3*pi*YY);   % solution of equation 
       %utrue = reshape(UTRUE,N1,1);  
       %RHS = (h*h) * 13*pi*pi*UTRUE;   % right of equation
    %   w2=1000;w1=0.001;
    %   UTRUE = sin(2*pi*w1*XX) .* sin(3*pi*w2*YY);   % solution of equation 
    %   utrue = reshape(UTRUE,N1,1);  
    %   RHS = (h*h) * (4*w1*w1+9*w2*w2)*pi*pi*UTRUE;   % right of equation
       
       UTRUE = (XX.*XX-XX.*XX.*XX.*XX).*(YY.*YY.*YY.*YY-YY.*YY);
       utrue = reshape(UTRUE,N1,1);
       RHS = (h*h)*2*((1-6*XX.*XX).*YY.*YY.*(1-YY.*YY)+...
                      (1-6*YY.*YY).*XX.*XX.*(1-XX.*XX)); 
       b = reshape(RHS,N1,1);
       generate_rhs = 0;
    end

    x0 = zeros(N1,1);
    fprintf('------------------------------------------------------------------------------------\n');
    fprintf('level=%d×%d,,',nx1+1,ny1+1);
    
    demo3_vmg(A1,b,x0,utrue,rtol,prtol,max_it,max_time,max_mflop);
    
        
end