| blob | 392901623a4be0377e2068e3c762ebfee490c1f4 |
1 function core_test(exe)\r
2 \r
3 \r
4 % defaults\r
5 if ~exist('exe','var'), exe='core_test_prog.exe';end\r
6 exe\r
7 % grid\r
8 xs=10; % domain size\r
9 ys=10;\r
10 %ys=12;\r
11 %m=111 % domain mesh\r
12 m=101\r
13 n=101\r
14 time0=0\r
15 time1=10\r
16 r=0\r
17 example='w'\r
18 \r
19 switch example\r
20 case {'circlespin','c'}\r
21 wind='c'\r
22 ignition='c'\r
23 case {'gostraight','g'}\r
24 wind='s'\r
25 ignition='c'\r
26 case {'windmill','w'}\r
27 wind='c'\r
28 ignition='w'\r
29 time1=5;\r
30 case {'spinstardemo','s'},\r
31 % same as spinStarDemo('low',1)\r
32 xs=2; % domain size\r
33 ys=2;\r
34 m=101 % domain mesh\r
35 n=101\r
36 % normal propagation speed\r
37 r=0.2;\r
38 % max wind speed\r
39 w=1;\r
40 wind='c'\r
41 ignition='s'\r
42 time0=0\r
43 time1=1\r
44 end\r
45 \r
46 phi=zeros(m,n); % allocate level function\r
47 \r
48 dx=xs/m;\r
49 dy=ys/n;\r
50 x=dx*[0:m-1];\r
51 y=dy*[0:n-1];\r
52 \r
53 switch wind\r
54 case {'straight','s'}\r
55 alpha=-1.3*pi\r
56 speed=2\r
57 vx=speed*cos(alpha)*ones(m,n);\r
58 vy=speed*sin(alpha)*ones(m,n);\r
59 case {'circular','c'}\r
60 [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
61 rotation=-0.75 * pi;\r
62 rotation=0.3 * pi;\r
63 vx=yy*rotation;\r
64 vy=-xx*rotation;\r
65 otherwise\r
66 error('no wind')\r
67 end\r
68 \r
69 switch ignition\r
70 case {'star','s'}\r
71 points = 7;\r
72 shift = 2.5;\r
73 scale = 0.20;\r
74 [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
75 [ theta, rad ] = cart2pol(xx, yy);\r
76 phi = rad - scale * (cos(points * theta) + shift);\r
77 case {'circle','c'}\r
78 c=[0.3*xs,0.5*ys]; % center\r
79 d=0.5; % radius\r
80 for i=1:m\r
81 for j=1:n\r
82 % phi(i,j)=sign(norm([i*dx,j*dy]-c)-d);\r
83 phi(i,j)=(norm([i*dx,j*dy]-c)-d);\r
84 end\r
85 end\r
86 case {'widmill','w'}\r
87 phi=ones(m,n);\r
88 phi(floor(m*0.20):ceil(m*0.45),floor(n*0.49):ceil(n*0.51))=-1;\r
89 phi(floor(m*0.55):ceil(m*0.80),floor(n*0.49):ceil(n*0.51))=-1;\r
90 phi(floor(m*0.49):ceil(m*0.51),floor(n*0.20):ceil(n*0.45))=-1;\r
91 phi(floor(m*0.49):ceil(m*0.51),floor(n*0.55):ceil(n*0.80))=-1;\r
92 otherwise\r
93 error('no ignition')\r
94 end\r
95 \r
96 plotSteps = 10; % How many intermediate plots to produce?\r
97 t0 = time0; % Start time.\r
98 \r
99 % Period at which intermediate plots should be produced.\r
100 tPlot = (time1 - time0) / (plotSteps - 1);\r
101 \r
102 %r=ones(m,n).*r + 0.04*sqrt(vx.*vx+vy.*vy); % shrinkage correction\r
103 r=ones(m,n).*r; % make sure it is array\r
104 tign=-1e10*ones(m,n);\r
105 fuel_time=20*ones(m-1,n-1);\r
106 tign(phi<=0)=time0;\r
107 frac_end=ones(m-1,n-1);\r
108 frac_lost=zeros(m-1,n-1);\r
109 \r
110 tNow=time0;\r
111 vis(phi,frac_lost,vx,vy,dx,dy,tNow);\r
112 \r
113 for i=1:plotSteps\r
114 tNext=min(time1, tNow + tPlot);\r
115 \r
116 % write the input file for the fortran code\r
117 fmt='%25.15e \n';\r
118 fid = fopen('core_test_in.txt','wt');\r
119 if(fid<0), error('cannot open file core_test_in.txt'),end\r
120 fprintf(fid,fmt,1,m,1,n,tNow,tNext,dx,dy);\r
121 fprintf(fid,fmt,phi,tign,fuel_time,r,vx,vy);\r
122 fclose(fid);\r
123 \r
124 % call the fortran implementation\r
125 eval(['! ',exe])\r
126 \r
127 % read the output file from the fortran code\r
128 fid = fopen('core_test_out.txt','rt');\r
129 if(fid<0), error('cannot open file core_test_out.txt'),end\r
130 [in,c]=fscanf(fid,'%g');\r
131 fclose(fid);\r
132 \r
133 % parse the file read into a single array\r
134 count=2*m*n+2*(m-1)*(n-1);\r
135 if(c~=count), error(sprintf('read %g terms should be %g\n',c,count)),end\r
136 last=0;\r
137 terms=m*n;phi=reshape(in(1+last:terms+last),size(phi));last=terms+last;\r
138 terms=m*n;tign=reshape(in(1+last:terms+last),size(tign));last=terms+last;\r
139 tign_d=tign;tign_d(tign<0)=NaN;\r
140 terms=(m-1)*(n-1);frac_lost=reshape(in(1+last:terms+last),size(frac_lost));last=terms+last;\r
141 terms=(m-1)*(n-1);frac_end=reshape(in(1+last:terms+last),size(frac_end));last=terms+last;\r
142 if(last~=count),error('bad count'),end\r
143 frac_tend=frac_lost/(tNext-tNow);\r
144 \r
145 % display\r
146 clf;figure(1)\r
147 vis(phi,-frac_tend,vx,vy,dx,dy,tNow)\r
148 %subplot(1,3,1);vis(phi,frac_tend,vx,vy,dx,dy,tNow)\r
149 %subplot(1,3,2);mesh(tign_d);title('tign')\r
150 %subplot(1,3,3);mesh(phi);title('phi')\r
151 drawnow\r
152 pause(0.3); % for ctrl-c\r
153 tNow = tNext;\r
154 end\r
155 \r
156 end\r
157 \r
158 %-------------------------------------------\r
159 \r
160 function vis(u,f,vx,vy,dx,dy,tNow)\r
161 [m,n]=size(u);\r
162 x=[0:m-1]*dx;\r
163 y=[0:n-1]*dy;\r
164 hold off\r
165 vis_type= '3d';\r
166 drawn=false;\r
167 switch vis_type\r
168 case '3d'\r
169 figure(1);clf\r
170 h=surf(y,x,u);\r
171 set(h,'edgecolor','none')\r
172 alpha(0.3)\r
173 hold on\r
174 contour3(y,x,u,[0 0],'k')\r
175 drawn=true;\r
176 figure(2);clf\r
177 xh=[1/2:m-3/2]*dx;\r
178 yh=[1/2:n-3/2]*dy;\r
179 h=pcolor(xh,yh,-f');\r
180 set(h,'edgecolor','none')\r
181 colorbar\r
182 \r
183 case '2d'\r
184 [m,n]=size(u);\r
185 axis equal\r
186 xh=[1/2:m-3/2]*dx;\r
187 yh=[1/2:n-3/2]*dy;\r
188 h=pcolor(xh,yh,-f');\r
189 set(h,'edgecolor','none')\r
190 colorbar\r
191 hold on\r
192 contour(y,x,u',[0 0],'k')\r
193 s=20;\r
194 ix=1:ceil(m/s):m;\r
195 iy=1:ceil(n/s):n;\r
196 [xx,yy]=ndgrid(x,y);\r
197 quiver(xx(ix,iy),yy(ix,iy),vx(ix,iy),vy(ix,iy))\r
198 drawn=true;\r
199 otherwise\r
200 disp('no visualization selected, set vis=type to 2d or 3d')\r
201 end\r
202 hold off\r
203 if drawn,\r
204 title(sprintf('t=%g',tNow));\r
205 end\r
206 end\r
207 \r