function nma_math228b_HW2_prob2() % implements the refinement study for HW2, ath 228B UC Davis % % and also does the part(d) verification of discrete conservation % property % % problem 2, Math 228B, Winter 2011, UC Davis % % By Nasser M. Abbasi 2/14/2011 %refinment_study(); verification_of_problem2_part_d(); end %-------------------------------------- function verification_of_problem2_part_d() % verification_of_problem2_part_d() This function verifies % that the solution to the discrete equation of the diffusion problem % in 2D satisfies the discete conservation property. % % This solves part (d) problem 2, HW 2 for Math 228B % UC Davis, Winter 2011 % % Algorithm: % Select run parameters % Call ADI solver with option to display sum(sum(U)) % over the whole 2D space at each time step. % Verify by looking at the result that the value above does not change % as run time is changing. % % Initialize variables before start grid_size = 3^4; %initial grid size, select one D = 0.1; %diffusion constant time_to_run = 1; %simulation time h = 1/grid_size; delt = h; %make time step the same as space step. ADI is second order in %in space and time. SHOW_SUM_AT_EACH_STEP = false; %flag to tell ADI solver to print values %make initial data ic = @(X,Y) exp( -10*((X-0.4).^2 + (Y-0.4).^2 )); [X,Y] = meshgrid(h/2:h:1-h/2,h/2:h:1-h/2); % Generate the A,B Matrices to use for ADI solver. THis is done once % and not changed, since the same grid size is used for the whole % simulation time [A,A_rhs]= ... nma_generate_A_and_ARHS_for_2D_diffusion_Neumman(grid_size,D,delt,h); %Call solver. The table will be printed on the screen since the flag below %is set to true nma_solve_2D_diffusion_ADI(A,A_rhs,h,delt,D,time_to_run,ic(X,Y),... SHOW_SUM_AT_EACH_STEP); % For future changes, make solver uses Lx,Ly on 1D only instead of the way % currently implemented and compare which is faster. Need more time. % A=nma_make_LX_Neumann_BC_1D(grid_size); % I=speye(grid_size); % diagA=-2*I; % solve_2D_diffusion_ADI(A,diagA,I,h,delt,D,time_to_run,ic(X,Y),SHOW_SUM_AT_EACH_STEP); end %-------------------------------------- function refinment_study() % refinment_study() is function which does refinement study for % problem 2, HW 2. % use 5 runs, and allocate the table to store the error and ratios % change as needed, but 5 seems enough to converge to ratio expected number_of_runs = 5; % Initialize variables before start grid_size = 1; %initial grid size h2 = figure(); %for log plot D = 0.1; %diffusion constant time_to_run = 1; %How long to run solver SKIP = 3; %factor to increase grid size by see inside loop R = refinement_study_manager(number_of_runs,SKIP); ic = @(X,Y) exp( -10*((X-0.4).^2 + (Y-0.4).^2 )); % initial data %ic = @(X,Y) cos(2*pi*X).*cos(2*pi*Y); %ic = @(X,Y) sin(X).*cos(Y); %ic= @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); SHOW_SUM_AT_EACH_STEP = false; %flag for ADI. for current_run_number = 1:number_of_runs % Simulatiously divide space step and time step. grid_size = grid_size * SKIP; % will go like 3,9,27,81,.... % change h and delt at the same time h = 1/grid_size; delt = h; u = solve_ADI(grid_size,h,delt,D,time_to_run,ic,SHOW_SUM_AT_EACH_STEP,h2); R.add(u, h, delt); end % The refinment study is completed. Generate plots and error table h2 = figure(); set(0,'CurrentFigure',h2); set(0,'defaultaxesfontsize',8) ; [hdr,body] = R.format_table(); fprintf('%s\n',hdr); fprintf('%s\n',body); table = R.get_table(); loglog(table(2:end,4),table(2:end,5),'-d'); xlabel('log(h)','FontSize',8); title({'refinement study result, ';'log vs successive errors difference'},... 'FontSize',8); ylabel('log(error norm)','FontSize',8); grid on; end %--------------------------------- function u = solve_ADI(grid_size,h,delt,D,time_to_run,ic,SHOW_SUM_AT_EACH_STEP,h2) % Generate Mesh and operators to use by ADI solver [X,Y]=meshgrid(h/2:h:1-h/2,h/2:h:1-h/2); [A,A_rhs]= ... nma_generate_A_and_ARHS_for_2D_diffusion_Neumman(grid_size,D,delt,h); %call the solver, it will return the numerical solution in u [u,X,Y,~] = nma_solve_2D_diffusion_ADI(A,A_rhs,h,... delt,D,time_to_run,ic(X,Y),SHOW_SUM_AT_EACH_STEP); % plot current solution set(0,'CurrentFigure',h2); surf(X,Y,u); colormap cool; title('Solution at end of run'); xlabel('x'); ylabel('y'); zlabel('u(x,y,t)'); drawnow(); end