| blob | de55250219c929dd520a3669cbf291323bd513d3 |
1 function prop_test\r
2 % a jm/mkim version of propagate\r
3 % follows Osher-Fedkiw and Mitchell toolbox\r
4 \r
5 % defaults\r
6 % grid\r
7 %xs=15; % domain size\r
8 xs=10;\r
9 ys=10;\r
10 %m=121 % domain mesh\r
11 m=101\r
12 n=101\r
13 time0=0\r
14 time1=10\r
15 r=0\r
16 speed=2 % wind speed\r
17 wind='c'\r
18 alpha=0.75*pi\r
19 dt=0.2;\r
20 plotSteps=50\r
21 normal_spread_c=1 % speed = r + wind^e * c\r
22 normal_spread_e=1\r
23 \r
24 example='p'\r
25 \r
26 switch example\r
27 case {'circlespin','c'}\r
28 wind='c'\r
29 ignition='c'\r
30 case {'gostraight','g'}\r
31 wind='s'\r
32 ignition='c'\r
33 case {'physical','p'}\r
34 wind='s';\r
35 ignition='c';\r
36 alpha=pi;\r
37 r=0.02;\r
38 m=401\r
39 n=401\r
40 %m=151\r
41 %n=151\r
42 speed=10;\r
43 xs=6*(m-1);\r
44 ys=6*(n-1);\r
45 time1=200;\r
46 dt=1;\r
47 plotSteps=400\r
48 normal_spread_c=0.185060861\r
49 normal_spread_e=1.310758329\r
50 %normal_spread_c=1\r
51 %normal_spread_e=1\r
52 case {'windmill','w'}\r
53 time1=3;\r
54 r=0.1;\r
55 wind='c'\r
56 ignition='w' \r
57 case {'spinstardemo','s'},\r
58 % same as spinStarDemo('low',1)\r
59 xs=2; % domain size\r
60 ys=2;\r
61 m=101 % domain mesh\r
62 n=101\r
63 % normal propagation speed\r
64 r=0.2;\r
65 % max wind speed\r
66 w=1;\r
67 wind='c'\r
68 ignition='s'\r
69 time0=0\r
70 time1=1\r
71 case {'o','only growth'}\r
72 speed=0\r
73 wind='s'\r
74 ignition='c'\r
75 r=0.1\r
76 end\r
77 \r
78 phi=zeros(m,n); % allocate level function\r
79 \r
80 dx=xs/(m-1);\r
81 dy=ys/(n-1);\r
82 x=dx*[0:m-1];\r
83 y=dy*[0:n-1];\r
84 \r
85 switch wind\r
86 case {'straight','s'}\r
87 vx=speed*cos(alpha)*ones(m,n);\r
88 vy=speed*sin(alpha)*ones(m,n);\r
89 case {'circular','c'}\r
90 [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
91 rotation=-0.75 * pi;\r
92 rotation=0.3 * pi;\r
93 vx=yy*rotation;\r
94 vy=-xx*rotation;\r
95 otherwise\r
96 error('no wind')\r
97 end\r
98 \r
99 switch ignition\r
100 case {'star','s'}\r
101 points = 7;\r
102 shift = 2.5;\r
103 scale = 0.20;\r
104 [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
105 [ theta, rad ] = cart2pol(xx, yy);\r
106 phi = rad - scale * (cos(points * theta) + shift);\r
107 case {'circle','c'}\r
108 c=[0.8*xs,0.5*ys]; % center\r
109 d=3*max(dx,dy); % radius\r
110 for i=1:m\r
111 for j=1:n\r
112 % phi(i,j)=sign(norm([i*dx,j*dy]-c)-d);\r
113 phi(i,j)=(norm([i*dx,j*dy]-c)-d);\r
114 end\r
115 end\r
116 case {'widmill','w'}\r
117 phi=ones(m,n);\r
118 phi(floor(m*0.20):ceil(m*0.45),floor(n*0.49):ceil(n*0.51))=-1;\r
119 phi(floor(m*0.55):ceil(m*0.80),floor(n*0.49):ceil(n*0.51))=-1;\r
120 phi(floor(m*0.49):ceil(m*0.51),floor(n*0.20):ceil(n*0.45))=-1;\r
121 phi(floor(m*0.49):ceil(m*0.51),floor(n*0.55):ceil(n*0.80))=-1;\r
122 otherwise\r
123 error('no ignition')\r
124 end\r
125 \r
126 \r
127 %r=ones(m,n).*r + 0.04*sqrt(vx.*vx+vy.*vy); % shrinkage correction\r
128 r=ones(m,n).*r; % make sure it is array\r
129 \r
130 data2=phi;\r
131 data3=phi;\r
132 tNow=time0;\r
133 dofort=1;\r
134 domat=0;\r
135 for i=1:plotSteps\r
136 tNext=tNow+dt\r
137 % % call the fortran implementation\r
138 if dofort,\r
139 [data3]=prop_test_f(data3,tNow,tNext,...\r
140 normal_spread_c,normal_spread_e,vx,vy,r,dx,dy);\r
141 end\r
142 %call own imp\r
143 %[ux,uy]=get_advection(phi,r,vx,vy,dx,dy);\r
144 %vis_wind(ux,uy,dx,dy)\r
145 if domat,\r
146 [dummy,data2]=prop_ls_cir(data2,tNow,tNext,vx,vy,r,dx,dy,@spread_rate);\r
147 end\r
148 %[tNow,data2]=prop_ls(data2,tNow,tNext,ux,uy,0,dx,dy);\r
149 if dofort & domat, \r
150 err_fort=norm(data3(:)-data2(:))\r
151 if err_fort>1e-10,\r
152 i,warning('large difference between matlab and fortran')\r
153 end\r
154 end\r
155 % display\r
156 if dofort,\r
157 vis(data3,vx,vy,dx,dy,tNow);\r
158 else\r
159 vis(data2,vx,vy,dx,dy,tNow);\r
160 end\r
161 tNow = tNext;\r
162 end\r
163 %err_fort=norm(data3(:)-data2(:))\r
164 % call the toolbox routines \r
165 % for that need to be in toolboxls/Examples/OsherFedkiw\r
166 if 0\r
167 data = normal_advection(phi,time0,time1,vx,vy,r,dx,dy,'low');\r
168 fprintf('\n');\r
169 err_toolbox=norm(data(:)-data2(:))\r
170 end \r
171 \r
172 end\r
173 \r
174 %-------------------------------------------\r
175 \r
176 function vis(u,vx,vy,dx,dy,tNow)\r
177 [m,n]=size(u);\r
178 x=[0:m-1]*dx;\r
179 y=[0:n-1]*dy;\r
180 figure(1);clf\r
181 hold off\r
182 vis_type= '2d';\r
183 drawn=false;\r
184 switch vis_type\r
185 case '3d'\r
186 h=surf(y,x,u);\r
187 set(h,'edgecolor','none')\r
188 alpha(0.3)\r
189 hold on\r
190 contour3(y,x,u,[0 0],'k')\r
191 hold off\r
192 drawn=true;\r
193 case '2d'\r
194 [m,n]=size(u);\r
195 contour(y,x,u,[0 0],'k')\r
196 hold on\r
197 vis_wind(vx,vy,dx,dy)\r
198 axis equal\r
199 drawn=true;\r
200 otherwise\r
201 disp('no visualization selected, set vis=type to 2d or 3d')\r
202 end\r
203 if drawn,\r
204 title(sprintf('t=%g',tNow));\r
205 drawnow;\r
206 pause(0.3); % for ctrl-c\r
207 end\r
208 end\r
209 \r
210 function [phi_out]=prop_test_f(phi,ts,te,c,e,vx,vy,r,dx,dy)\r
211 fmt='%25.15e \n';\r
212 if any(any(isnan(phi))), error('phi NaN'), end\r
213 [m,n]=size(phi);\r
214 fid = fopen('prop_test_in.txt','wt');\r
215 fprintf(fid,fmt,m,n,ts,te,c,e,r,dx,dy);\r
216 fprintf(fid,fmt,phi,vx,vy);\r
217 fclose(fid);\r
218 ! prop_test_prog.exe\r
219 fid = fopen('prop_test_out.txt','rt');\r
220 if fid < 0,\r
221 error('cannot open prop_test_out.txt')\r
222 end\r
223 [in,c]=fscanf(fid,'%g');\r
224 fclose(fid);\r
225 if c~=m*n,\r
226 fprintf('read %g terms should have %g\n',c,m*n+1)\r
227 error('bad number of terms in output file')\r
228 end\r
229 %t=in(1);\r
230 phi_out=reshape(in(1:c),size(phi));\r
231 if any(any(isnan(phi_out))), error('phi_out NaN'), end\r
232 end\r
233 \r
234 function speed=spread_rate(r,vx,vy,nvx,nvy,scale);\r
235 speed = r + max(vx.*nvx + vy.*nvy,0);\r
236 end\r
237 \r