function nma_math_228b_HW4_problem_3() %solves diffusion problem $u_t + a * u_x = 0$ using finite volume method with flux limiter functions %This function solves problem 3 HW 4 for Math 228B, UC Davis, winter 2011 % %It solves the diffusion problem u_t + a * u_x = 0 using %finite volume method with the use of flux limiter functions % %by Nasser M. Abbasi %March 18, 2011 % close all; nma_set_figure_position(); %adjust figure position for display doingRefinement = false; %flag to use if doing refinment study if doingRefinement % % Do refinement study. Use refinement_study object % % the grid space to use, each time it is half the last time h = [1/8 1/16 1/32 1/64 1/128 1/256 1/512 1/1024]; num_refinement_study_runs = length(h); %these are the methods to do refinement study on methods = {'upwind','lax','mc','superbee','minmod','vanleer','warming'}; tags = {'-^','-o','-s','-v','-+','-*'}; %for plotting for method=1:length(methods) %create a refinement study object R = refinement_study_manager(num_refinement_study_runs,1); for j = 1:num_refinement_study_runs [numerical_solution,exactSolution,k] = run(h(j),methods{method}); R.add_exact(numerical_solution, exactSolution,h(j),k); [hd,bdy] = R.format_table(); %hd %bdy if j == num_refinement_study_runs %get the table from the refinement object and plot it tbl = R.get_table(); figure(5); hold on; set(0,'defaultaxesfontsize',8) ; loglog(tbl(4:end,1),tbl(4:end,6),tags{method}); xlabel('number of points','FontSize',7); title('refinement study result, ','FontSize',7); ylabel('1-norm of error','FontSize',7); grid on; set(gca,'FontSize',7); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); drawnow(); end end end figure(5); legend(methods); else run(0.005,'vanleer'); end end %---------------------------- function u_initial=initial_data_as_cell_averages(X,data) %This function is called to initialize the data. It finds %the average value of the function for each cell %it uses the function u0() which is the initial data function u_initial = zeros(1,length(X)-1); for i=1:length(X)-1 u_initial(i) = quad(@(x) u0(x,data) ,X(i),X(i+1)) / data.h; end end %------------------------------ function y=u00(X,shift,data) %This function is also used for initialization, but for plotting %only, as it handles shifting the functions so that data appears %on the same place. %These just tell what initial data to use SIN=1; HANDOUT=2; PARTA=3; PARTB=4; PARTC=5; switch data.IC case SIN y=sin(4*pi*X); case HANDOUT y=(X > (0.1-shift)).*(X < (0.3-shift)) + ... exp(-200*(X-(0.75-shift)).^2); case PARTA y=cos(16*pi*(X)).*exp(-50*(X-(0.5)).^2); case PARTB y=sin(2*pi*X).*sin(4*pi*X); case PARTC y=(X > (.25-shift)).*(X < (0.75-shift)) ; end end %---------------------------- function y=u0(X,data) %this function is called to initialize the domain SIN=1; HANDOUT=2; PARTA=3; PARTB=4; PARTC=5; switch data.IC case SIN y=sin(4*pi*X); case HANDOUT y=(X > 0.1).*(X < 0.3) + ... exp(-200*(X-0.75).^2); case PARTA y=cos(16*pi*X).*exp((-50*(X-0.5).^2)); case PARTB y=sin(2*pi*X).*sin(4*pi*X); case PARTC y=(X > .25).*(X < 0.75) ; end end %---------------------------- function [u,u_exact,delt] = run(h,method) %This is the main function where the algorithm is run % %INPUT: % h: grid spacing % method: a string which tells which flux limiter function to use % %OUTPUT % u: numerical solution % u_exact: exact solution % delt: time delta used % %These are the initial conditions SIN=1; HANDOUT=2; PARTA=3; PARTB=4; PARTC=5; data.IC = HANDOUT; %change this to select different initial data data.method = method; data.h = h; data.a = 1; data.L = 1; data.courant_number = 0.9; X = 0:data.h:data.L; X_with_ghost_points = [X(1)-2*data.h X(1)-data.h X X(end)+data.h X(end)+2*data.h]; Xh = X_with_ghost_points(1:end-1)+data.h/2; u_initial0 = initial_data_as_cell_averages(X, data); data.y_limits = find_y_limits_to_use(u_initial0); u_initial = [u_initial0(end-2) u_initial0(end-1) u_initial0 u_initial0(2) u_initial0(3)]; u = u_initial; %solution goes here data.X = h/2:data.h:data.L-h/2; %grid with no ghost points data.delt = data.courant_number*data.h/data.a; delt = data.delt; data.max_t = 1; data.current_t = 0; set(gcf,'doublebuffer', 'on') ; data.ANIMATE = false; %flag to do animation data.current_time = 0; makePlot(Xh,u,0,data); %initial solution at t=0 nSteps = data.max_t/data.delt; %finds how many steps to run %this calculates if we need a partial step at the end completeNumberOfSteps = floor(nSteps); if (abs(nSteps-completeNumberOfSteps)>2*eps) partialStep = (nSteps - completeNumberOfSteps)*data.delt; else partialStep = 0; end % % start the main loop, check if doing animation % if data.ANIMATE f = getframe(gcf); [im,map] = rgb2ind(f.cdata,512,'nodither'); end animation_counter = 0; for n = 1:completeNumberOfSteps data.current_time = n*data.delt; [u,u_exact] = makeStep(u,Xh,n,data); % Solution step if data.ANIMATE if mod(n,10) ==0 f = getframe(gcf); animation_counter = animation_counter + 1; im(:,:,1,animation_counter) = rgb2ind(f.cdata,map,'nodither'); end end end %check if need to make partial step at end if partialStep ~=0 data.current_time = data.max_t; data.delt = partialStep; u = makeStep(u,Xh,n,data); end u=u(3:end-2); %actual solution without ghost points if data.ANIMATE imwrite(im,map,'problem_3.gif','DelayTime',0,'LoopCount',0); end end %---------------------- function u_exact = makePlot(Xh,u_new,stepNumber,data) %Called to make plot. This function also calcuates the exact %solution at this time % SIN=1; HANDOUT=2; PARTA=3; PARTB=4; PARTC=5; shift = data.a*data.current_time; %Depending on intitial data, handle the shifting part to allow the display %to appear fixed in space. For nonperiodic initial data we do something %different, hence need to check what was the initial data switch data.IC case {SIN , PARTA, PARTB} u_exact=u00(Xh(3:end-2),shift,data); case {HANDOUT , PARTC} u_exact=u00(Xh(3:end-2)-shift,shift,data); end if data.ANIMATE if mod(stepNumber,10) ~=0 return; end end %this below is a trick to handle data to make it appear fixed in space domain_length = Xh(end-2)-Xh(3); adjust_X = mod(Xh(3:end-2)-shift,domain_length); A = [adjust_X' u_new(3:end-2)']; A = sortrows(A); plot(A(:,1),A(:,2),'ro-'); xlim([Xh(1) Xh(end)]); hold on; plot( Xh(3:end-2),u_exact,'-','LineWidth',1); title({sprintf('method=%s, time=%f, step =%d, time step=%f',... data.method,data.current_time,stepNumber,data.delt);... sprintf('space step=%f, speed=%f, Courant number=%2.2f',... data.h,data.a,data.courant_number)},... 'FontSize',10,'interpreter','latex'); ylim(data.y_limits); xlabel('x (meter)'); ylabel('solution'); legend('numerical','exact'); grid; %gridxy(Xh(1:8:end),[]) ; hold off; set(gca,'FontSize',8); drawnow(); end %----------------------- function [u,u_exact] = makeStep(u,Xh,stepNumber,data) % Called to make solution step % u_new = zeros(size(u)); %this is lax-wendroff in finite volume using flux limiters for j=3:length(Xh)-2 u_new(j) = u(j) - (data.delt/data.h)* ( F('right',u,j,data) - ... F('left' ,u,j,data) ); end u_exact = makePlot(Xh,u_new,stepNumber,data); u = u_new; %update boundary points, periodic u(2) = u_new(end-3); u(1) = u_new(end-4); u(end-1) = u_new(4); u(end) = u_new(5); end %----------------------- function u_new = make_step_lax(u,data) %Not used. this function find Lax-Wendoff solution no finite volume %method. Used for verification. % i = 2:length(data.X)-1; u_new = zeros(size(u)); A = data.a*data.delt/(2*data.h); B = data.a^2*data.delt^2/(2*data.h^2); u_new(i) = u(i)-(A)*(u(i+1)-u(i-1)) + B*(u(i-1)-2*u(i)+u(i+1)); u_new(1) = u(1)-A*(u(2)-u(end-1))+B*(u(end-1)-2*u(1)+u(2)); u_new(end) = u(end)-A*(u(2)-u(end-1))+B*(u(end-1)-2*u(end)+u(2)); end %-------------------- function r = minmod(b) r=max(0,min(b,1)); end %----------------------------------- function r = PHI(val,method) %The PHI function, the flux limiter function switch method case 'lax' r = 1; case 'superbee' r = max([0,min(2*val,1),min(2,val)]); case 'minmod' r = minmod(val); case 'upwind' r = 0; case 'mc' r = max(0,min( [(1+val)/2,2,2*val])); case 'vanleer' r = (val+abs(val))./(1+abs(val)); case 'warming' r = val; end end %-------------------------- function r = F(edge,u,j,data) %Find the F function depending on the edge switch edge case 'left' du = u(j) - u(j-1); case 'right' du = u(j+1) - u(j) ; end if abs(du)<2*eps limited_diff = 0; else limited_diff = PHI ( theta(edge,j,u,data.a,du), data.method ) * du; end r = F_upwind(edge,j,data.a,u) + abs(data.a)/2 * ... (1- abs(data.a*data.delt/data.h))*limited_diff; end %---------------------------- function r = theta(edge,j,u,a,du) %evaluates the THETA function switch edge case 'left' if a>0 r = u(j-1)-u(j-2); else r = u(j+1)-u(j); end case 'right' if a>0 r = u(j) - u(j-1) ; else r = u(j+2)-u(j+1); end end r = r/du; end %----------------------------- function r = F_upwind(edge,j,a,u) %Evaluates F upwind switch edge case 'left' if a>0 r = a*u(j-1); else r = a*u(j); end case 'right' if a>0 r = a*u(j); else r = a*u(j+1); end end end %---------------------------- function y_limits=find_y_limits_to_use(u_initial0) %THis function finds good y-axis limits to use for plotting %depending on the initial data % ymin = min(u_initial0); ymax = max(u_initial0); if abs(ymin)<2*eps ymin= -0.2; end if abs(ymax)<2*eps ymax= 0.1; end y_limits = [1.2*ymin 1.3*ymax]; end