| blob | 6d9752fccea48463050d25e900c369a2382e0c73 |
1 function fuel_fraction=fuel_burnt(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 % set up least squares(A*u-v)*inv(C)*(A*u-v)->min\r
38 A=[0, 0, 1;\r
39 fd(1), 0, 1;\r
40 0, fd(2), 1;\r
41 fd(1),fd(2), 1];\r
42 v=tnow-tign(:) ; % time from ignition\r
43 rw=2*ones(4,1);\r
44 u = lscov(A,v,rw); % solve least squares to get coeffs of T\r
45 residual=norm(A*u-v); % zero if T is linear\r
46 u\r
47 display('fuel_burnt coefficients of u')\r
48 % integrate\r
49 uu = -u / fuel_time;\r
50 fuel_fraction = 1 - exp(uu(3)) * intexp(uu(1)*fd(1)) * intexp(uu(2)*fd(2));\r
51 if(fuel_fraction<0 | fuel_fraction>1),\r
52 warning('fuel_fraction should be between 0 and 1')\r
53 stop\r
54 end\r
55 else\r
56 % part of cell is burning - the interesting case\r
57 % set up a list of points with the corresponding values of T and L\r
58 xylist=zeros(8,2); % allocate list of points\r
59 Tlist=zeros(8,1);\r
60 Llist=zeros(8,1);\r
61 xx=[-fd(1) fd(1) fd(1) -fd(1) -fd(1)]/2; % make center (0,0)\r
62 yy=[-fd(2) -fd(2) fd(2) fd(2) -fd(2)]/2; % cyclic, counterclockwise\r
63 ii=[1 2 2 1 1]; % indices of corners, cyclic, counterclockwise\r
64 jj=[1 1 2 2 1];\r
65 points=0;\r
66 for k=1:4\r
67 lfn0=lfn(ii(k),jj(k));\r
68 lfn1=lfn(ii(k+1),jj(k+1));\r
69 if(lfn0<=0),\r
70 points=points+1;\r
71 xylist(points,:)=[xx(k),yy(k)]; % add corner to list\r
72 Tlist(points)=tnow-tign(ii(k),jj(k)); % time since ignition\r
73 Llist(points)=lfn0;\r
74 end\r
75 if(lfn0*lfn1<0),\r
76 points=points+1;\r
77 % coordinates of intersection of fire line with segment k k+1\r
78 %lfn(t)=lfn0 + t*(lfn1-lfn0)=0\r
79 t=lfn0/(lfn0-lfn1);\r
80 x0=xx(k)+(xx(k+1)-xx(k))*t;\r
81 y0=yy(k)+(yy(k+1)-yy(k))*t;\r
82 xylist(points,:)=[x0,y0];\r
83 Tlist(points)=0; % now at ignition\r
84 Llist(points)=0; % at fireline\r
85 end\r
86 end\r
87 % make the lists circular and trim to size\r
88 Tlist(points+1)=Tlist(1);\r
89 Tlist=Tlist(1:points+1);\r
90 Llist(points+1)=Llist(1);\r
91 Llist=Llist(1:points+1);\r
92 xylist(points+1,:)=xylist(1,:);\r
93 xylist=xylist(1:points+1,:);\r
94 for k=1:5,lfnk(k)=lfn(ii(k),jj(k));end\r
95 figure(1)\r
96 plot3(xx,yy,lfnk,'k')\r
97 hold on\r
98 plot3(xylist(:,1),xylist(:,2),Tlist,'m--o')\r
99 plot3(xylist(:,1),xylist(:,2),Llist,'g--o')\r
100 plot(xx,yy,'b')\r
101 plot(xylist(:,1),xylist(:,2),'-.r+')\r
102 legend('lfn on whole cell','tign-tnow on burned area',...\r
103 'lfn on burned area', 'cell boundary','burned area boundary')\r
104 patch(xylist(:,1),xylist(:,2),zeros(points+1,1),250,'FaceAlpha',0.3)\r
105 patch(xylist(:,1),xylist(:,2),Tlist,100,'FaceAlpha',0.3)\r
106 hold off,grid on,drawnow,pause(0.5)\r
107 % set up least squares\r
108 A=[xylist(1:points,1:2),ones(points,1)];\r
109 v=Tlist(1:points);\r
110 for i=1:points\r
111 for j=1:points\r
112 if(j~=i),\r
113 dist(j)=norm(xylist(i,:)-xylist(j,:));\r
114 else\r
115 dist(j)=max(fd); % large\r
116 end\r
117 end\r
118 rw(i)=1+min(dist); % weight = 1+min dist from other nodes\r
119 end\r
120 u = lscov(A,v,rw); % solve least squares to get coeffs of T\r
121 u\r
122 display('fuel_burnt coefficients of u')\r
123 residual=norm(A*u-v); % should be zero if T and L are linear\r
124 nr=sqrt(u(1)*u(1)+u(2)*u(2));\r
125 c=u(1)/nr;\r
126 s=u(2)/nr;\r
127 Q=[c, s; -s, c]; % rotation such that Q*u(1:2)=[something;0]\r
128 ut=Q*u(1:2);\r
129 errQ=ut(2); % should be zero\r
130 ae=-ut(1)/fuel_time;\r
131 ce=-u(3)/fuel_time; % -T(xt,yt)/fuel_time=ae*xt+ce\r
132 cet=ce; %keep ce for later\r
133 xytlist=xylist*Q'; % rotate the points in the list\r
134 xt=sort(xytlist(1:points,1)); % sort ascending in x\r
135 fuel_fraction=0;\r
136 for k=1:points-1\r
137 % integrate the vertical slice from xt1 to xt2\r
138 figure(2)\r
139 plot(xytlist(:,1),xytlist(:,2),'-o'),grid on,hold on\r
140 xt1=xt(k);\r
141 xt2=xt(k+1);\r
142 plot([xt1,xt1],sum(fd)*[-1,1]/3,'y')\r
143 plot([xt2,xt2],sum(fd)*[-1,1]/3,'y')\r
144 if(xt2-xt1>100*eps*max(fd)), % slice of nonzero width\r
145 % find slice height as h=ah*x+ch\r
146 upper=0;lower=0;\r
147 ah=0;ch=0;\r
148 for s=1:points % pass counterclockwise\r
149 xts=xytlist(s,1); % start point of the line\r
150 yts=xytlist(s,2);\r
151 xte=xytlist(s+1,1); % end point of the line\r
152 yte=xytlist(s+1,2);\r
153 if (xts>xt1 & xte > xt1) | (xts<xt2 & xte < xt2)\r
154 plot([xts,xte],[yts,yte],'--k')\r
155 % that is not the one\r
156 else % line y=a*x+c through (xts,yts), (xte,yte)\r
157 a=(yts-yte)/(xts-xte);\r
158 c=(xts*yte-xte*yts)/(xts-xte);\r
159 if xte<xts % upper boundary\r
160 aupp=a;\r
161 cupp=c;\r
162 plot([xts,xte],[yts,yte],'g')\r
163 ah=ah+a;\r
164 ch=ch+c;\r
165 upper=upper+1;\r
166 else % lower boundary\r
167 alow=a;\r
168 clow=c;\r
169 plot([xts,xte],[yts,yte],'m')\r
170 lower=lower+1;\r
171 end\r
172 end\r
173 end\r
174 if(lower~=1|upper~=1),\r
175 error('slice does not have one upper and one lower line')\r
176 end\r
177 ce=cet;% use kept ce\r
178 % shift samll amounts to avoid negative fuel burnt\r
179 if ae*xt1+ce > 0;%disp('ae*xt1+ce');disp(ae*xt1+ce);pause;\r
180 ce=ce-(ae*xt1+ce);end;\r
181 if ae*xt2+ce > 0;%disp('ae*xt2+ce');disp(ae*xt2+ce);pause;\r
182 ce=ce-(ae*xt2+ce);end;\r
183 \r
184 ah=aupp-alow;\r
185 ch=cupp-clow;\r
186 % debug only\r
187 patch([xt1,xt2,xt2,xt1],...\r
188 [alow*[xt1,xt2]+clow,aupp*[xt2,xt1]+cupp,],k*10)\r
189 axis equal,drawnow, pause(0.5)\r
190 figure(3)\r
191 x=[xt1:(xt2-xt1)/100:xt2];\r
192 plot(x,1-exp(ae*x+ce),x,ah*x+ch,x,...\r
193 (ah*x+ch).*(1-exp(ae*x+ce)),[xt1,xt2],[0,0],'+k')\r
194 xlabel x,legend('burned frac','slice height','height*burned','dividers')\r
195 hold on,grid on,drawnow,pause(0.5)\r
196 % integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2\r
197 % numerically sound for ae->0, ae -> infty\r
198 % this can be important for different model scales\r
199 % esp. if someone runs the model in single precision!!\r
200 % s1=int((ah*x+ch),x,xt1,xt2)\r
201 s1=(xt2-xt1)*((1/2)*ah*(xt2+xt1)+ch);\r
202 % s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)\r
203 ceae=ce/ae;\r
204 s2=-ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1));\r
205 % s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)\r
206 % s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)\r
207 % expand in Taylor series around ae=0\r
208 % collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)\r
209 % =(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
210 % + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae\r
211 % + 1/2*xt1^2-1/2*xt2^2\r
212 %\r
213 % coefficient at ae^2 in the expansion, after some algebra\r
214 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
215 d=ae^4*a2;\r
216 if abs(d)>eps\r
217 % since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae\r
218 % for ae, ce -> 0 rounding error approx eps/ae^2\r
219 s3=(exp(ae*(xt1+ceae))*(ae*xt1-1)-exp(ae*(xt2+ceae))*(ae*xt2-1))/ae^2;\r
220 % we do not worry about rounding as xt1 -> xt2, then s3 -> 0\r
221 else\r
222 % coefficient at ae^1 in the expansion\r
223 a1=(xt1-xt2)*((1/2)*ceae*(xt1+xt2)+(1/3)*(xt1^2+xt1*xt2+xt2^2));\r
224 % coefficient at ae^0 in the expansion for ae->0\r
225 a0=(1/2)*(xt1-xt2)*(xt1+xt2);\r
226 s3=a0+a1*ae+a2*ae^2; % approximate the integral\r
227 end\r
228 s3=ah*s3;\r
229 fuel_fraction=fuel_fraction+s1+s2+s3;\r
230 if(fuel_fraction<0 | fuel_fraction>(fd(1)*fd(2))),\r
231 fuel_fraction,(fd(1)*fd(2)),s1,s2,s3\r
232 warning('fuel_fraction should be between 0 and 1')\r
233 stop\r
234 end\r
235 end\r
236 end\r
237 fuel_fraction=fuel_fraction/(fd(1)*fd(2));\r
238 display('fuel_burnt2')\r
239 %fuel_fraction=1-fuel_fraction;\r
240 end\r
241 for i=figs,figure(i),hold off,end\r
242 end % function\r
243 \r
244 function s=intexp(ab)\r
245 % function s=intexp(ab)\r
246 % s = (1/b)*int(exp(a*x),x,0,b) ab = a*b\r
247 % = 1 if a==0\r
248 if eps < abs(ab)^3/6,\r
249 s = (exp(ab)-1)/(ab); % rounding error approx eps/(a*b) for small a*b\r
250 else\r
251 s = 1 + ab/2+ab^2/6; % last term (a*b)^2/6 for small a*b\r
252 end\r
253 end\r