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