comments in fuel_left_cell_2, more checking there
[wrffire.git] / standalone / matlab / core_test.m
blob392901623a4be0377e2068e3c762ebfee490c1f4
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
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
46 phi=zeros(m,n); % allocate level function\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
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
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
96 plotSteps = 10;               % How many intermediate plots to produce?\r
97 t0 = time0;                      % Start time.\r
99 % Period at which intermediate plots should be produced.\r
100 tPlot = (time1 - time0) / (plotSteps - 1);\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
110 tNow=time0;\r
111 vis(phi,frac_lost,vx,vy,dx,dy,tNow);\r
113 for i=1:plotSteps\r
114     tNext=min(time1, tNow + tPlot);\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
124     %   call the fortran implementation\r
125     eval(['! ',exe])\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
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
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
156 end\r
158 %-------------------------------------------\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
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