/standalone/matlab/fuel_burnt.m
Objective C | 254 lines | 251 code | 3 blank | 0 comment | 36 complexity | d9b13a2d14f8d264180b16d0c42a1466 MD5 | raw file
Possible License(s): AGPL-1.0
- function fuel_fraction=fuel_burnt(lfn,tign,tnow,fd,fuel_time)
- % lfn - level set function, size (2,2)
- % tign - time of ignition, size (2,2)
- % tnow - time now
- % fd - mesh cell size (2)
- % fuel_time - time constant of fuel
- %
- % output: approximation of
- % /\
- % 1 | T(x,y)
- % fuel_burnt = ---------- | ( 1 - exp( - ----------- ) ) dxdy
- % fd(1)*fd(2) | fuel_time
- % \/
- % (x,y) in R;
- % L(x,y)<=0
- %
- % where R is the rectangle [0,fd(1)] * [0,fd(2)]
- % T is given by tign and L is given by lfn at the 4 vertices of R
- %
- % note that fuel_left = 1 - fuel_burnt
- %
- % Requirements:
- % exact in the case when T and L are linear
- % varies continuously with input
- % value of tign-tnow where lfn>0 ignored
- % assume T=0 when lfn=0
- figs=1:4;
- for i=figs,figure(i),clf,hold off,end
- % tign
- % lfn
- if all(lfn(:)>=0),
- % nothing burning, nothing to do - most common case, put it first
- fuel_fraction=0;
- elseif all(lfn(:)<=0),
- % all burning
- % T=u(1)*x+u(2)*y+u(3)
- % set up least squares(A*u-v)*inv(C)*(A*u-v)->min
- A=[0, 0, 1;
- fd(1), 0, 1;
- 0, fd(2), 1;
- fd(1),fd(2), 1];
- v=tnow-tign(:) ; % time from ignition
- rw=2*ones(4,1);
- u = lscov(A,v,rw); % solve least squares to get coeffs of T
- residual=norm(A*u-v); % zero if T is linear
- u
- display('fuel_burnt coefficients of u')
- % integrate
- uu = -u / fuel_time;
- fuel_fraction = 1 - exp(uu(3)) * intexp(uu(1)*fd(1)) * intexp(uu(2)*fd(2));
- if(fuel_fraction<0 | fuel_fraction>1),
- warning('fuel_fraction should be between 0 and 1')
- stop
- end
- else
- % part of cell is burning - the interesting case
- % set up a list of points with the corresponding values of T and L
- xylist=zeros(8,2); % allocate list of points
- Tlist=zeros(8,1);
- Llist=zeros(8,1);
- xx=[-fd(1) fd(1) fd(1) -fd(1) -fd(1)]/2; % make center (0,0)
- yy=[-fd(2) -fd(2) fd(2) fd(2) -fd(2)]/2; % cyclic, counterclockwise
- ii=[1 2 2 1 1]; % indices of corners, cyclic, counterclockwise
- jj=[1 1 2 2 1];
- points=0;
- for k=1:4
- lfn0=lfn(ii(k),jj(k));
- lfn1=lfn(ii(k+1),jj(k+1));
- if(lfn0<=0),
- points=points+1;
- xylist(points,:)=[xx(k),yy(k)]; % add corner to list
- Tlist(points)=tnow-tign(ii(k),jj(k)); % time since ignition
- Llist(points)=lfn0;
- end
- if(lfn0*lfn1<0),
- points=points+1;
- % coordinates of intersection of fire line with segment k k+1
- %lfn(t)=lfn0 + t*(lfn1-lfn0)=0
- t=lfn0/(lfn0-lfn1);
- x0=xx(k)+(xx(k+1)-xx(k))*t;
- y0=yy(k)+(yy(k+1)-yy(k))*t;
- xylist(points,:)=[x0,y0];
- Tlist(points)=0; % now at ignition
- Llist(points)=0; % at fireline
- end
- end
- % make the lists circular and trim to size
- Tlist(points+1)=Tlist(1);
- Tlist=Tlist(1:points+1);
- Llist(points+1)=Llist(1);
- Llist=Llist(1:points+1);
- xylist(points+1,:)=xylist(1,:);
- xylist=xylist(1:points+1,:);
- for k=1:5,lfnk(k)=lfn(ii(k),jj(k));end
- figure(1)
- plot3(xx,yy,lfnk,'k')
- hold on
- plot3(xylist(:,1),xylist(:,2),Tlist,'m--o')
- plot3(xylist(:,1),xylist(:,2),Llist,'g--o')
- plot(xx,yy,'b')
- plot(xylist(:,1),xylist(:,2),'-.r+')
- legend('lfn on whole cell','tign-tnow on burned area',...
- 'lfn on burned area', 'cell boundary','burned area boundary')
- patch(xylist(:,1),xylist(:,2),zeros(points+1,1),250,'FaceAlpha',0.3)
- patch(xylist(:,1),xylist(:,2),Tlist,100,'FaceAlpha',0.3)
- hold off,grid on,drawnow,pause(0.5)
- % set up least squares
- A=[xylist(1:points,1:2),ones(points,1)];
- v=Tlist(1:points);
- for i=1:points
- for j=1:points
- if(j~=i),
- dist(j)=norm(xylist(i,:)-xylist(j,:));
- else
- dist(j)=max(fd); % large
- end
- end
- rw(i)=1+min(dist); % weight = 1+min dist from other nodes
- end
- u = lscov(A,v,rw); % solve least squares to get coeffs of T
- u
- display('fuel_burnt coefficients of u')
- residual=norm(A*u-v); % should be zero if T and L are linear
- nr=sqrt(u(1)*u(1)+u(2)*u(2));
- c=u(1)/nr;
- s=u(2)/nr;
- Q=[c, s; -s, c]; % rotation such that Q*u(1:2)=[something;0]
- ut=Q*u(1:2);
- errQ=ut(2); % should be zero
- ae=-ut(1)/fuel_time;
- ce=-u(3)/fuel_time; % -T(xt,yt)/fuel_time=ae*xt+ce
- cet=ce; %keep ce for later
- xytlist=xylist*Q'; % rotate the points in the list
- xt=sort(xytlist(1:points,1)); % sort ascending in x
- fuel_fraction=0;
- for k=1:points-1
- % integrate the vertical slice from xt1 to xt2
- figure(2)
- plot(xytlist(:,1),xytlist(:,2),'-o'),grid on,hold on
- xt1=xt(k);
- xt2=xt(k+1);
- plot([xt1,xt1],sum(fd)*[-1,1]/3,'y')
- plot([xt2,xt2],sum(fd)*[-1,1]/3,'y')
- if(xt2-xt1>100*eps*max(fd)), % slice of nonzero width
- % find slice height as h=ah*x+ch
- upper=0;lower=0;
- ah=0;ch=0;
- for s=1:points % pass counterclockwise
- xts=xytlist(s,1); % start point of the line
- yts=xytlist(s,2);
- xte=xytlist(s+1,1); % end point of the line
- yte=xytlist(s+1,2);
- if (xts>xt1 & xte > xt1) | (xts<xt2 & xte < xt2)
- plot([xts,xte],[yts,yte],'--k')
- % that is not the one
- else % line y=a*x+c through (xts,yts), (xte,yte)
- a=(yts-yte)/(xts-xte);
- c=(xts*yte-xte*yts)/(xts-xte);
- if xte<xts % upper boundary
- aupp=a;
- cupp=c;
- plot([xts,xte],[yts,yte],'g')
- ah=ah+a;
- ch=ch+c;
- upper=upper+1;
- else % lower boundary
- alow=a;
- clow=c;
- plot([xts,xte],[yts,yte],'m')
- lower=lower+1;
- end
- end
- end
- if(lower~=1|upper~=1),
- error('slice does not have one upper and one lower line')
- end
- ce=cet;% use kept ce
- % shift samll amounts to avoid negative fuel burnt
- if ae*xt1+ce > 0;%disp('ae*xt1+ce');disp(ae*xt1+ce);pause;
- ce=ce-(ae*xt1+ce);end;
- if ae*xt2+ce > 0;%disp('ae*xt2+ce');disp(ae*xt2+ce);pause;
- ce=ce-(ae*xt2+ce);end;
-
- ah=aupp-alow;
- ch=cupp-clow;
- % debug only
- patch([xt1,xt2,xt2,xt1],...
- [alow*[xt1,xt2]+clow,aupp*[xt2,xt1]+cupp,],k*10)
- axis equal,drawnow, pause(0.5)
- figure(3)
- x=[xt1:(xt2-xt1)/100:xt2];
- plot(x,1-exp(ae*x+ce),x,ah*x+ch,x,...
- (ah*x+ch).*(1-exp(ae*x+ce)),[xt1,xt2],[0,0],'+k')
- xlabel x,legend('burned frac','slice height','height*burned','dividers')
- hold on,grid on,drawnow,pause(0.5)
- % integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
- % numerically sound for ae->0, ae -> infty
- % this can be important for different model scales
- % esp. if someone runs the model in single precision!!
- % s1=int((ah*x+ch),x,xt1,xt2)
- s1=(xt2-xt1)*((1/2)*ah*(xt2+xt1)+ch);
- % s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
- ceae=ce/ae;
- s2=-ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1));
- % s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
- % s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
- % expand in Taylor series around ae=0
- % collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
- % =(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
- % + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
- % + 1/2*xt1^2-1/2*xt2^2
- %
- % coefficient at ae^2 in the expansion, after some algebra
- 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));
- d=ae^4*a2;
- if abs(d)>eps
- % since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
- % for ae, ce -> 0 rounding error approx eps/ae^2
- s3=(exp(ae*(xt1+ceae))*(ae*xt1-1)-exp(ae*(xt2+ceae))*(ae*xt2-1))/ae^2;
- % we do not worry about rounding as xt1 -> xt2, then s3 -> 0
- else
- % coefficient at ae^1 in the expansion
- a1=(xt1-xt2)*((1/2)*ceae*(xt1+xt2)+(1/3)*(xt1^2+xt1*xt2+xt2^2));
- % coefficient at ae^0 in the expansion for ae->0
- a0=(1/2)*(xt1-xt2)*(xt1+xt2);
- s3=a0+a1*ae+a2*ae^2; % approximate the integral
- end
- s3=ah*s3;
- fuel_fraction=fuel_fraction+s1+s2+s3;
- if(fuel_fraction<0 | fuel_fraction>(fd(1)*fd(2))),
- fuel_fraction,(fd(1)*fd(2)),s1,s2,s3
- warning('fuel_fraction should be between 0 and 1')
- stop
- end
- end
- end
- fuel_fraction=fuel_fraction/(fd(1)*fd(2));
- display('fuel_burnt2')
- %fuel_fraction=1-fuel_fraction;
- end
- for i=figs,figure(i),hold off,end
- end % function
-
- function s=intexp(ab)
- % function s=intexp(ab)
- % s = (1/b)*int(exp(a*x),x,0,b) ab = a*b
- % = 1 if a==0
- if eps < abs(ab)^3/6,
- s = (exp(ab)-1)/(ab); % rounding error approx eps/(a*b) for small a*b
- else
- s = 1 + ab/2+ab^2/6; % last term (a*b)^2/6 for small a*b
- end
- end
-