function FVlinearadvectionFOU2D % File name: FVlinearadvectionFOU2D.m % Author: Clive Mingham % Date: 24 Feb 2011 % Description: Solves the 2D PDE, dU/dt + div(H) = 0 % using the FOU finite volume scheme. % This equation corresponds to the finite volume form: % doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0. % % In this case the equation is the linear advection equation so: % % flux density: H = vU % velocity: v = vx i + vy j = <vx,vy> where vx and vy are constant. % % The program is written in a structured way using subfunctions so that it % can easily be modified for other equations, schemes and meshes. % subfunctions: freadmesh, fcellarea, fsidevectors, finitialu, % fIflux, fghostcells, fcalcdt, fplotresults, fdisplay, % fdrawmesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Verification: (on uniform Cartesian mesh) % for NI=61,NJ=81,dx=2,dy=4, initial profile centre is at (61,82) % contour plots of numerical agree with exact for: % vx=10,vy=0, t=1.04, new centre (71, 82) is correct % vx=0, vy=10, t=1.12, new centre (61, 94) is correct % vx=10,vy=10, t=1.04, new centre (71, 94) is correct % % Comments: Although program appears to work the numerical scheme is highly % dissipative. %% clc; clf; clear all; runtime=1.0; % target runtime [s] t=0.0; % time [s] tlevel=0; % time level vx=10; % velocity component in x direction [m/s] vy=10; % velocity component in y direction [m/s] %vx=0;vy=10; v=[vx,vy]; % velocity vector <vx,vy> [x,y,xcen,ycen,solid]=freadmesh; % gets coordinates of the 2D computational mesh % cell centres and solid flags (excluding ghost cells) disp('read in mesh') [m,n]=size(x); NI=m-1; % number of computational cells in i direction. NJ=n-1; % number of computational cells in j direction. A=fcellarea(x,y); % compute and store cell areas in an NIxNJ array disp('calculated cell areas') S=fsidevectors(x,y); % compute and store cell side vectors in an % (NI+2)x(NJ+2)x4x2 array which contains ghost cells. disp('calculated cell side vectors') u=finitialu(xcen,ycen,0,v); % put initial cell centre values of u in NIxNJ array uinitial=u; % store for plotting initial conditions disp('inserted initial conditions') % % Append extra cells around arrays to store any ghost values. u=finsertghostcells(u); % u is now (NI+2)x(NJ+2) unext=zeros(size(u)); % (NI+2)x(NJ+2) array for u values at next time level A=finsertghostcells(A); % A is now (NI+2)x(NJ+2) %% disp('time marching starts') while(t<runtime) u=fbcs(u); % Implement boundary conditions using ghost cells. % u is (NI+2)x(NJ+2) and each computational cell is % indexed i=2 to NI+1, j=2 to NJ+1. dt=fcalcdt(A,S,v); % Finds stable time step for each iteration. t=t+dt % update time tlevel=tlevel+1; for j=2:NJ+1 for i=2:NI+1 IH=fIflux(v,u,i,j); % gets interface fluxes for cell (i,j) IH1=[IH(1,1),IH(1,2)]; % interface flux vector for side 1 IH2=[IH(2,1),IH(2,2)]; % interface flux vector for side 2 IH3=[IH(3,1),IH(3,2)]; % interface flux vector for side 3 IH4=[IH(4,1),IH(4,2)]; % interface flux vector for side 4 % s1=[S(i,j,1,1),S(i,j,1,2)]; % side vector for side 1 s2=[S(i,j,2,1),S(i,j,2,2)]; % side vector for side 2 s3=[S(i,j,3,1),S(i,j,3,2)]; % side vector for side 3 s4=[S(i,j,4,1),S(i,j,4,2)]; % side vector for side 4 % FV scheme % compute total flux out of cell (i,j) totalfluxout=dot(IH1,s1)+dot(IH2,s2)+dot(IH3,s3)+dot(IH4,s4); unext(i,j)=u(i,j)-(dt/A(i,j))*totalfluxout; end % of i loop end % of j loop u(2:NI+1,2:NJ+1)=unext(2:NI+1,2:NJ+1); % update u values % end % of while loop disp('time marching ends') uexact=finitialu(xcen,ycen,t,v); % exact solution at time t u=u(2:NI+1,2:NJ+1); % extract final computational values fdrawmesh(x,y,solid) % draws mesh if fplotresults is commented out %fdisplay(u); % print results as a matrix for a small mesh if needed fplotresults(xcen,ycen,uinitial,u,uexact) % plot results: comment in or out disp('finished main program, see figures') end % FVlinearadvectionFOU2D %%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x,y,xcen,ycen,solid]=freadmesh % The mesh is structured and has NI cells in the i direction and % NJ cells in the j direction (so each cell has 4 sides). % The x and y coordinates of the lower left hand corner of cell (i,j) are % held in arrays x and y respectively which are both (NI+1)x(NJ+1). In % this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), % (x(i+1),y(j+1)) and (x(i),y(j+1)). xcen and ycen are NI by NJ arrays % which hold cell centres. solid is an NI by NJ array which % flags solid cells. If cell (i,j) is solid then solid(i,j)=1 otherwise % solid(i,j)=0. % % % test using diamond % x(1,1)=0; y(1,1)=-1; % x(2,1)=1; y(2,1)=0; % x(2,2)=0; y(2,2)=1; % x(1,2)=-1; y(1,2)=0; % % Test using regular Cartesian mesh starting at (0,0). i direction is % x direction. % % NI=61,NJ=41,dx=2,dy=4 makes mesh centre (61,82) NI=61; NJ=41; dx=2; dy=4; x=zeros(NI+1,NJ+1); % allocate correct sized array y=zeros(NI+1,NJ+1); % allocate correct sized array xcen=zeros(NI,NJ); % allocate correct sized array ycen=zeros(NI,NJ); % allocate correct sized array solid=zeros(NI,NJ); % allocate correct sized array % create cell vertices for i=1:NI+1 for j=1:NJ+1 x(i,j)=(i-1)*dx; y(i,j)=(j-1)*dy; end % of j loop end % of i loop nonCartesian=1; % 1 for non-Cartesian mesh if (nonCartesian==1) disp('mesh is non-Cartesian') % Add 'random' numbers to mesh coordinates to produce a non-Cartesian mesh % ensure that the abs of these numbers are less than min(dx/2,dy/2). % Create non-Cartesian mesh for i=1:NI+1 for j=1:NJ+1 x(i,j)=x(i,j)+(dx/6)*cos(1.727*(i+j)); y(i,j)=y(i,j)+(dy/5)*sin(2.46723*i*j); end % of j loop end % of i loop else disp('mesh is Cartesian') end % end of if statement create mesh % % find cell centres by averaging coordinates of vertices for i=1:NI for j=1:NJ xcen(i,j)=(x(i,j)+x(i+1,j)+x(i+1,j+1)+x(i,j+1))/4; ycen(i,j)=(y(i,j)+y(i+1,j)+y(i+1,j+1)+y(i,j+1))/4; end % of j loop end % of i loop end % freadmesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A=fcellarea(x,y) % Verification: correct for single diamond % cell area = half the modulus of cross product of cell diagonal vectors. % These vectors are 2D so append 0 k component to make 3D for cross product [m,n]=size(x); NI=m-1; NJ=n-1; A=zeros(NI,NJ); for i=1:NI for j=1:NJ d1=[x(i+1,j+1)-x(i,j),y(i+1,j+1)-y(i,j),0]; d2=[x(i,j+1)-x(i+1,j),y(i,j+1)-y(i+1,j),0]; d=cross(d1,d2); A(i,j)=0.5*sqrt(dot(d,d)); end % of j loop end % of i loop end % fcellarea %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function S=fsidevectors(x,y) % Verification: correct for single diamond shaped cell. % % Finds the 4 side (S) vectors for each cell. % Cell sides are numbered 1 to 4 anti-clockwise. For cell (i,j) side 1 % is the vector s1 from (x(i,j), y(i,j)) to (x(i+1,j),y(i+1,j)). % if s1=<a,b> then the corresponding side vector S1=<b,-a>. % The side vector array has ghost cells around it so that side vector % indices (i+1,j+1) refer to cell (i,j). % % Note that vectors are 1 by 2 arrays so care must be taken because side % vectors are stored in higher dimensional arrays. % [m,n]=size(x); NI=m-1; NJ=n-1; S=zeros(NI+2,NJ+2,4,2); for i=1:NI for j=1:NJ s1=[x(i+1,j)-x(i,j),y(i+1,j)-y(i,j)]; S(i+1,j+1,1,1)=s1(2); S(i+1,j+1,1,2)=-s1(1); % s2=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)]; S(i+1,j+1,2,1)=s2(2); S(i+1,j+1,2,2)=-s2(1); % s3=[x(i,j+1)-x(i+1,j+1),y(i,j+1)-y(i+1,j+1)]; S(i+1,j+1,3,1)=s3(2); S(i+1,j+1,3,2)=-s3(1); % s4=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)]; S(i+1,j+1,4,1)=s4(2); S(i+1,j+1,4,2)=-s4(1); end % of j loop end % of i loop % [S(2,2,1,1),S(2,2,1,2)] % [S(2,2,2,1),S(2,2,2,2)] they are correct for a diamond % [S(2,2,3,1),S(2,2,3,2)] % [S(2,2,4,1),S(2,2,4,2)] end % fsidevectors %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function u=finitialu(xcen,ycen,t,v) % Inserts exact u values (and initial conditions for t=0) [NI,NJ]=size(xcen); u=zeros(NI,NJ); % create correct sized array (without ghost values) vx=v(1); vy=v(2); % % Initial 2D Gaussian function based on the centre (xc,yc) of the % computational domain. xmax=max(max(xcen)); % max x value in mesh ymax=max(max(ycen)); % max y value in mesh xmin=min(min(xcen)); % min x value in mesh ymin=min(min(ycen)); % min y value in mesh xc=(xmax+xmin)/2; % approx x coord of mesh centre yc=(ymax+ymin)/2; % approx y coord of mesh centre % cutoff=min(xmax-xc,ymax-yc)/2; % cut off radius for Gaussian for i=1:NI for j=1:NJ x=xcen(i,j)-vx*t; y=ycen(i,j)-vy*t; d=sqrt((x-xc)^2+(y-yc)^2); % distance from centre if (d<cutoff) u(i,j)=exp(-0.01*d^2); else u(i,j)=0; end end % of j loop end % of i loop end % finitialu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function IH=fIflux(v,u,i,j) % Calculates the 4 interface fluxes on the 4 sides of cell (i,j). % Cell sides are numbered 1 to 4 anti-clockwise. For cell (i,j) side 1 % is the vector s1 from (x(i,j), y(i,j)) and to (x(i+1,j),y(i+1,j)). % Interface fluxes are denoted IH1, IH2, IH3 and IH4 and are 2D vectors. % Flux H depends on U, i.e. H = H(U). % For the 2D linear advection equation, H(U) = <vx U, vy U>. % % There is considerable choice for interface flux estimation: different % choices give different FV schemes. % %% The following scheme simply takes: % IH1=H(u(i,j-1)) = v u(i,j-1) % IH2=H(u(i,j)) = v u(i,j) % IH3=H(u(i,j)) = v u(i,j) % IH4=H(u(i-1,j)) = v u(i-1,j) % This is an upwind scheme for vx > 0 and vy >0 and corresponds to the % FOU finite difference scheme on a Cartesian mesh. % IH=zeros(4,2); % array to store all 4 interface flux vectors for a cell % vx=v(1); vy=v(2); % IH(1,1)=vx*u(i,j-1); % component of IH1 in the x direction IH(1,2)=vy*u(i,j-1); % component of IH1 in the y direction IH(2,1)=vx*u(i,j); % component of IH2 in the x direction IH(2,2)=vy*u(i,j); % component of IH2 in the y direction IH(3,1)=vx*u(i,j); % component of IH3 in the x direction IH(3,2)=vy*u(i,j); % component of IH3 in the y direction IH(4,1)=vx*u(i-1,j); % component of IH4 in the x direction IH(4,2)=vy*u(i-1,j); % component of IH4 in the y direction end % fIflux %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function outarray=finsertghostcells(inarray) % Assuming that ghost cells are needed around the whole domain an mxn % array is embedded into (m+2)x(n+2) array of zeros. Hence computational % indices go from i=2 to m+1 and j=2 to n+1. [m,n]=size(inarray); dummy=zeros(m+2,n+2); dummy(2:m+1,2:n+1)=inarray; outarray=dummy; end % finsertghostcells %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function outarray=fremoveghostcells(inarray) % % Removes ghost cells so that a mxn array becomes (m-2)x(n-2) % [m,n]=size(inarray); % outarray=inarray(2:m-1,2:n-1); % end % fremoveghostcells %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function u=fbcs(u) % Implements boundary conditions using ghost cells index by i=1, i=NI+2 % j=1, j=NJ+2. [m,n]=size(u); NI=m-2; NJ=n-2; % % Dirichlet: u=0 everywhere. % Don't need to do anything as u values in ghost cells are already zero. end % bcs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dt=fcalcdt(A,S,v) % Verification: % Finds allowable time step dt using heuristic formula (which makes sense % 0<F<1 is a tunable safety factor. [m,n]=size(A); % A and S now enlarged to include ghost cells NI=m-2; NJ=n-2; dt=9999999; % large value of dt to start F=0.4; for i=2:NI+1 for j=2:NJ+1 % put side vectors into 1x2 arrays s2=[S(i,j,2,1),S(i,j,2,2)]; s3=[S(i,j,3,1),S(i,j,3,2)]; % v2=abs(dot(v,s2)); % side 2 v3=abs(dot(v,s3)); % side 3 dum=A(i,j)/max(v2,v3); % take smallest value of dt if (dum<dt) dt=dum; end end % j loop end % i loop dt=F*dt; end % fcalcdt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function fdrawmesh(x,y,solid) % Description: Draws a structured 2D finite volume mesh and fills in any % solid cells. % Date structures: % The mesh has NI cells in the i direction and NJ cells in the j direction. % The x and y coordinates of the lower left hand corner of cell (i,j) are % held in arrays x and y respectively which are both NI+1 by NJ+1. In % this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), % (x(i+1),y(j+1)) and (x(i),y(j+1)). solid is an NI by NJ array which % flags solid cells. If cell (i,j) is solid then solid(i,j)=1 otherwise % solid(i,j)=0. % [m,n]=size(x); NI=m-1; % number of cells in i direction NJ=n-1; % number of cells in j direction % % Plot the mesh hold on % keeps all plots on the same axes % draw lines in the i direction for i=1:NI+1 plot(x(i,:),y(i,:)) end % draw lines in the j direction for j=1:NJ+1 plot(x(:,j),y(:,j)) end title('computational mesh') xlabel('x') ylabel('y') % Fill in solid cells for i=1:NI for j=1:NJ if (solid(i,j)==1) solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)]; solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)]; fill(solidx,solidy,'k') end end % of j loop end % of i loop end % fdrawmesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function fplotresults(x,y,uinitial,u,uexact) % Plots results on a structured mesh - even works for non-Cartesian! % subplot(3,1,1), contour(x,y,uinitial) title('Contour plot of initial conditions') xlabel('x [m]') ylabel('y [m]') % subplot(3,1,2), contour(x,y,u) title('Contour plot of numerical solution') xlabel('x [m]') ylabel('y [m]') % subplot(3,1,3), contour(x,y,uexact) title('Contour plot of exact solution') xlabel('x [m]') ylabel('y [m]') % end % fplotresults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function fdisplay(in) % displays matrices 'correctly' - WORKS % Using index i for x which are columns and j for y which are rows in % reverse order so this routine displays them in this fashion [NX,NY]=size(in); temp=transpose(in); out=temp(NY:-1:1,:); disp('printing as a matrix') out end % fdisplay %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%