| blob | 9148c6899df90d601ce8ef2f8107c4860de294a0 |
1 function fuel_fraction=fuel_burnt_fd(lfn,tign,tnow,fd,fuel_time)\r
2 % lfn - level set function, size (2,2)\r
3 % tign - time of ignition, size (2,2)\r
4 % tnow - time now\r
5 % fd - mesh cell size (2)\r
6 % fuel_time - time constant of fuel\r
7 %\r
8 % output: approximation of\r
9 % /\\r
10 % 1 | T(x,y)\r
11 % fuel_burnt = ---------- | ( 1 - exp( - ----------- ) ) dxdy\r
12 % fd(1)*fd(2) | fuel_time\r
13 % \/\r
14 % (x,y) in R;\r
15 % L(x,y)<=0\r
16 %\r
17 % where R is the rectangle [0,fd(1)] * [0,fd(2)]\r
18 % T is given by tign and L is given by lfn at the 4 vertices of R\r
19 %\r
20 % note that fuel_left = 1 - fuel_burnt\r
21 %\r
22 % Requirements:\r
23 % exact in the case when T and L are linear\r
24 % varies continuously with input\r
25 % value of tign-tnow where lfn>0 ignored\r
26 % assume T=0 when lfn=0\r
27 figs=1:4;\r
28 for i=figs,figure(i),clf,hold off,end\r
29 % tign\r
30 % lfn\r
31 if all(lfn(:)>=0),\r
32 % nothing burning, nothing to do - most common case, put it first\r
33 fuel_fraction=0;\r
34 elseif all(lfn(:)<=0),\r
35 % all burning\r
36 % T=u(1)*x+u(2)*y+u(3)\r
37 % t(0,0)=tign(1,1)\r
38 % t(0,msh_sz)=tign(1,2)\r
39 % t(msh_sz,0)=tign(2,1)\r
40 % t(msh_sz,msh_sz)=tign(2,2)\r
41 % msh_sz=h\r
42 % t(g/2,h/2)=sum(tign(i,i))/4\r
43 % dt/dx=(1/2h)*(t10-t00+t11-t01)\r
44 % dt/dy=(1/2h)*(t01-t00+t11-t10)\r
45 % approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences\r
46 % t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy\r
47 % u(1)=dt/dx\r
48 % u(2)=dt/dy\r
49 % u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)\r
50 \r
51 tign_middle=(tign(1,1)+tign(1,2)+tign(2,1)+tign(2,2))/4;\r
52 dt_dx=(tign(1,1)-tign(2,1)+tign(1,2)-tign(2,2))/(2*fd(1));\r
53 dt_dy=(tign(1,1)-tign(1,2)+tign(2,1)-tign(2,2))/(2*fd(2));\r
54 % dt_dx=(tign(2,1)-tign(1,1)+tign(2,2)-tign(1,2))/(2*fd(1));\r
55 % dt_dy=(tign(1,2)-tign(1,1)+tign(2,2)-tign(2,1))/(2*fd(2));\r
56 \r
57 u(1)=dt_dx;\r
58 u(2)=dt_dy;\r
59 u(3)=tign_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2;\r
60 u1=u'\r
61 display('fuel_burnt_fd coefficients of u')\r
62 % integrate\r
63 uu = -u / fuel_time;\r
64 fuel_fraction = 1 - exp(uu(3)) * intexp(uu(1)*fd(1)) * intexp(uu(2)*fd(2));\r
65 if(fuel_fraction<0 | fuel_fraction>1),\r
66 warning('fuel_fraction should be between 0 and 1')\r
67 stop\r
68 end\r
69 else\r
70 % part of cell is burning - the interesting case\r
71 % set up a list of points with the corresponding values of T and L\r
72 xylist=zeros(8,2); % allocate list of points\r
73 Tlist=zeros(8,1);\r
74 Llist=zeros(8,1);\r
75 xx=[-fd(1) fd(1) fd(1) -fd(1) -fd(1)]/2; % make center (0,0)\r
76 yy=[-fd(2) -fd(2) fd(2) fd(2) -fd(2)]/2; % cyclic, counterclockwise\r
77 ii=[1 2 2 1 1]; % indices of corners, cyclic, counterclockwise\r
78 jj=[1 1 2 2 1];\r
79 points=0;\r
80 for k=1:4\r
81 lfn0=lfn(ii(k),jj(k));\r
82 lfn1=lfn(ii(k+1),jj(k+1));\r
83 if(lfn0<=0),\r
84 points=points+1;\r
85 xylist(points,:)=[xx(k),yy(k)]; % add corner to list\r
86 Tlist(points)=tnow-tign(ii(k),jj(k)); % time since ignition\r
87 Llist(points)=lfn0;\r
88 end\r
89 if(lfn0*lfn1<0),\r
90 points=points+1;\r
91 % coordinates of intersection of fire line with segment k k+1\r
92 %lfn(t)=lfn0 + t*(lfn1-lfn0)=0\r
93 t=lfn0/(lfn0-lfn1);\r
94 x0=xx(k)+(xx(k+1)-xx(k))*t;\r
95 y0=yy(k)+(yy(k+1)-yy(k))*t;\r
96 xylist(points,:)=[x0,y0];\r
97 Tlist(points)=0; % now at ignition\r
98 Llist(points)=0; % at fireline\r
99 end\r
100 end\r
101 % make the lists circular and trim to size\r
102 Tlist(points+1)=Tlist(1);\r
103 Tlist=Tlist(1:points+1);\r
104 Llist(points+1)=Llist(1);\r
105 Llist=Llist(1:points+1);\r
106 xylist(points+1,:)=xylist(1,:);\r
107 xylist=xylist(1:points+1,:);\r
108 for k=1:5,lfnk(k)=lfn(ii(k),jj(k));end\r
109 figure(1)\r
110 plot3(xx,yy,lfnk,'k')\r
111 hold on\r
112 plot3(xylist(:,1),xylist(:,2),Tlist,'m--o')\r
113 plot3(xylist(:,1),xylist(:,2),Llist,'g--o')\r
114 plot(xx,yy,'b')\r
115 plot(xylist(:,1),xylist(:,2),'-.r+')\r
116 legend('lfn on whole cell','tign-tnow on burned area',...\r
117 'lfn on burned area', 'cell boundary','burned area boundary')\r
118 patch(xylist(:,1),xylist(:,2),zeros(points+1,1),250,'FaceAlpha',0.3)\r
119 patch(xylist(:,1),xylist(:,2),Tlist,100,'FaceAlpha',0.3)\r
120 hold off,grid on,drawnow,pause(0.5)\r
121 %%%%% Approximation of the plane for lfn using finite differences\r
122 % approximate L(x,y)=u(1)*x+u(2)*y+u(3)\r
123 lfn_middle=(lfn(1,1)+lfn(1,2)+lfn(2,1)+lfn(2,2))/4;\r
124 dt_dx=(lfn(1,1)-lfn(2,1)+lfn(1,2)-lfn(2,2))/(2*fd(1));\r
125 dt_dy=(lfn(1,1)-lfn(1,2)+lfn(2,1)-lfn(2,2))/(2*fd(2));\r
126 % I feel that it is right but not least squares tnow-tign\r
127 u(1)=dt_dx;\r
128 u(2)=dt_dy;\r
129 u(3)=lfn_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2;\r
130 % finding the coefficient c, reminder we work over one subcell only\r
131 % T(x,y)=c*L(x,y) in the sence of least squares\r
132 a=0; b=0;\r
133 for i=1:2\r
134 for j=1:2\r
135 if (lfn(i,j)<=0) \r
136 a=a+lfn(i,j)*lfn(i,j);\r
137 if (tign(i,j)>tnow)\r
138 disp('fuel_burnt_fd: tign(i1) should be less then time_now');\r
139 else\r
140 b=b+(tign(i,j)-tnow)*lfn(i,j);\r
141 end\r
142 end\r
143 end\r
144 end\r
145 if (a==0)\r
146 display('fuel_burnt_fd: if c is on fire then one of cells should be on fire');\r
147 end\r
148 c=b/a;\r
149 u(1)=u(1)*c;\r
150 u(2)=u(2)*c;\r
151 u(3)=u(3)*c;\r
152 u=u'\r
153 display('fuel_burnt_fd coefficients of u')\r
154 nr=sqrt(u(1)*u(1)+u(2)*u(2));\r
155 c=u(1)/nr;\r
156 s=u(2)/nr;\r
157 Q=[c, s; -s, c]; % rotation such that Q*u(1:2)=[something;0]\r
158 ut=Q*u(1:2);\r
159 errQ=ut(2); % should be zero\r
160 ae=-ut(1)/fuel_time;\r
161 ce=-u(3)/fuel_time; % -T(xt,yt)/fuel_time=ae*xt+ce\r
162 cet=ce; %keep ce for later\r
163 xytlist=xylist*Q'; % rotate the points in the list\r
164 xt=sort(xytlist(1:points,1)); % sort ascending in x\r
165 fuel_fraction=0;\r
166 for k=1:points-1\r
167 % integrate the vertical slice from xt1 to xt2\r
168 figure(2)\r
169 plot(xytlist(:,1),xytlist(:,2),'-o'),grid on,hold on\r
170 xt1=xt(k);\r
171 xt2=xt(k+1);\r
172 plot([xt1,xt1],sum(fd)*[-1,1]/3,'y')\r
173 plot([xt2,xt2],sum(fd)*[-1,1]/3,'y')\r
174 if(xt2-xt1>100*eps*max(fd)), % slice of nonzero width\r
175 % find slice height as h=ah*x+ch\r
176 upper=0;lower=0;\r
177 ah=0;ch=0;\r
178 for s=1:points % pass counterclockwise\r
179 xts=xytlist(s,1); % start point of the line\r
180 yts=xytlist(s,2);\r
181 xte=xytlist(s+1,1); % end point of the line\r
182 yte=xytlist(s+1,2);\r
183 if (xts>xt1 & xte > xt1) | (xts<xt2 & xte < xt2)\r
184 plot([xts,xte],[yts,yte],'--k')\r
185 % that is not the one\r
186 else % line y=a*x+c through (xts,yts), (xte,yte)\r
187 a=(yts-yte)/(xts-xte);\r
188 c=(xts*yte-xte*yts)/(xts-xte);\r
189 if xte<xts % upper boundary\r
190 aupp=a;\r
191 cupp=c;\r
192 plot([xts,xte],[yts,yte],'g')\r
193 ah=ah+a;\r
194 ch=ch+c;\r
195 upper=upper+1;\r
196 else % lower boundary\r
197 alow=a;\r
198 clow=c;\r
199 plot([xts,xte],[yts,yte],'m')\r
200 lower=lower+1;\r
201 end\r
202 end\r
203 end\r
204 if(lower~=1|upper~=1),\r
205 error('slice does not have one upper and one lower line')\r
206 end\r
207 ce=cet;% use kept ce\r
208 % shift samll amounts to avoid negative fuel burnt\r
209 if ae*xt1+ce > 0;%disp('ae*xt1+ce');disp(ae*xt1+ce);pause;\r
210 ce=ce-(ae*xt1+ce);end;\r
211 if ae*xt2+ce > 0;%disp('ae*xt2+ce');disp(ae*xt2+ce);pause;\r
212 ce=ce-(ae*xt2+ce);end;\r
213 \r
214 ah=aupp-alow;\r
215 ch=cupp-clow;\r
216 % debug only\r
217 patch([xt1,xt2,xt2,xt1],...\r
218 [alow*[xt1,xt2]+clow,aupp*[xt2,xt1]+cupp,],k*10)\r
219 axis equal,drawnow, pause(0.5)\r
220 figure(3)\r
221 x=[xt1:(xt2-xt1)/100:xt2];\r
222 plot(x,1-exp(ae*x+ce),x,ah*x+ch,x,...\r
223 (ah*x+ch).*(1-exp(ae*x+ce)),[xt1,xt2],[0,0],'+k')\r
224 xlabel x,legend('burned frac','slice height','height*burned','dividers')\r
225 hold on,grid on,drawnow,pause(0.5)\r
226 % integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2\r
227 % numerically sound for ae->0, ae -> infty\r
228 % this can be important for different model scales\r
229 % esp. if someone runs the model in single precision!!\r
230 % s1=int((ah*x+ch),x,xt1,xt2)\r
231 s1=(xt2-xt1)*((1/2)*ah*(xt2+xt1)+ch);\r
232 % s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)\r
233 ceae=ce/ae;\r
234 s2=-ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1));\r
235 % s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)\r
236 % s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)\r
237 % expand in Taylor series around ae=0\r
238 % collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)\r
239 % =(1/8*xt1^4+1/3*xt1^3*ceae+1/4*xt1^2*ceae^2-1/8*xt2^4-1/3*xt2^3*ceae-1/4*xt2^2*ceae^2)*ae^2\r
240 % + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae\r
241 % + 1/2*xt1^2-1/2*xt2^2\r
242 %\r
243 % coefficient at ae^2 in the expansion, after some algebra\r
244 a2=(xt1-xt2)*((1/4)*(xt1+xt2)*ceae^2+(1/3)*(xt1^2+xt1*xt2+xt2^2)*ceae+(1/8)*(xt1^3+xt1*xt2^2+xt1^2*xt2+xt2^3));\r
245 d=ae^4*a2;\r
246 if abs(d)>eps\r
247 % since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae\r
248 % for ae, ce -> 0 rounding error approx eps/ae^2\r
249 s3=(exp(ae*(xt1+ceae))*(ae*xt1-1)-exp(ae*(xt2+ceae))*(ae*xt2-1))/ae^2;\r
250 % we do not worry about rounding as xt1 -> xt2, then s3 -> 0\r
251 else\r
252 % coefficient at ae^1 in the expansion\r
253 a1=(xt1-xt2)*((1/2)*ceae*(xt1+xt2)+(1/3)*(xt1^2+xt1*xt2+xt2^2));\r
254 % coefficient at ae^0 in the expansion for ae->0\r
255 a0=(1/2)*(xt1-xt2)*(xt1+xt2);\r
256 s3=a0+a1*ae+a2*ae^2; % approximate the integral\r
257 end\r
258 s3=ah*s3;\r
259 fuel_fraction=fuel_fraction+s1+s2+s3;\r
260 if(fuel_fraction<0 | fuel_fraction>(fd(1)*fd(2))),\r
261 fuel_fraction,(fd(1)*fd(2)),s1,s2,s3\r
262 warning('fuel_fraction should be between 0 and 1')\r
263 stop\r
264 end\r
265 end\r
266 end\r
267 fuel_fraction=fuel_fraction/(fd(1)*fd(2));\r
268 % display('fuel_burnt_fd')\r
269 %fuel_fraction=1-fuel_fraction;\r
270 end\r
271 for i=figs,figure(i),hold off,end\r
272 end % function\r
273 \r
274 function s=intexp(ab)\r
275 % function s=intexp(ab)\r
276 % s = (1/b)*int(exp(a*x),x,0,b) ab = a*b\r
277 % = 1 if a==0\r
278 if eps < abs(ab)^3/6,\r
279 s = (exp(ab)-1)/(ab); % rounding error approx eps/(a*b) for small a*b\r
280 else\r
281 s = 1 + ab/2+ab^2/6; % last term (a*b)^2/6 for small a*b\r
282 end\r
283 end\r