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