From c4134d76c368dead79152fadc73971af616ba34b Mon Sep 17 00:00:00 2001 From: Volodymyr Kondratenko Date: Tue, 27 Jul 2010 23:29:31 -0600 Subject: [PATCH] fuel_burnt_fd.m is done in Matlab --- standalone/matlab/fuel_burnt_fd.m | 504 +++++++++++++++++++++----------------- standalone/matlab/fuel_test.m | 2 +- 2 files changed, 277 insertions(+), 229 deletions(-) rewrite standalone/matlab/fuel_burnt_fd.m (78%) diff --git a/standalone/matlab/fuel_burnt_fd.m b/standalone/matlab/fuel_burnt_fd.m dissimilarity index 78% index 0ca20ce1..4a9340f4 100644 --- a/standalone/matlab/fuel_burnt_fd.m +++ b/standalone/matlab/fuel_burnt_fd.m @@ -1,228 +1,276 @@ -function fuel_fraction=fuel_burnt_fd(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 -% mesh is from fd(1) to fd(2), so its length is fd(2)-fd(1) -if (fd(2)0, but fd(2)-fd(1)<0 '); -else - if all(lfn(:)>=0), - % nothing burning, nothing to do - most common case, put it first - fuel_fraction=0; - else -% t(0,0)=tign(1,1) -% t(0,msh_sz)=tign(1,2) -% t(msh_sz,0)=tign(2,1) -% t(msh_sz,msh_sz)=tign(2,2) -% msh_sz=h -% t(g/2,h/2)=sum(tign(i,i))/4 -% dt/dx=(1/2h)*(t10-t00+t11-t01) -% dt/dy=(1/2h)*(t01-t00+t11-t10) -% approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences -% t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy -% u(1)=dt/dx -% u(2)=dt/dy -% u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy) - - msh_sz=fd(2)-fd(1); - tign_middle=(tign(1,1)+tign(1,2)+tign(2,1)+tign(2,2))/4; - dt_dx=(tign(2,1)-tign(1,1)+tign(2,2)-tign(1,2))/(2*msh_sz); - dt_dy=(tign(1,2)-tign(1,1)+tign(2,2)-tign(2,1))/(2*msh_sz); - u(1)=dt_dx; - u(2)=dt_dy; - u(3)=tign_middle-msh_sz*(dt_dx+dt_dy)/2; - if all(lfn(:)<=0), - % all burning - % 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 -% 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: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 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 -% ends -% 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)) -% %fuel_fraction=1-fuel_fraction; -% end -% for i=figs,figure(i),hold off,end -% end % function - - - - end - - end -end -end -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 -% intexp=1/b*integral(0,b) of exp(ax)dx -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 - +function fuel_fraction=fuel_burnt_fd2(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 + +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) + % t(0,0)=tign(1,1) +% t(0,msh_sz)=tign(1,2) +% t(msh_sz,0)=tign(2,1) +% t(msh_sz,msh_sz)=tign(2,2) +% msh_sz=h +% t(g/2,h/2)=sum(tign(i,i))/4 +% dt/dx=(1/2h)*(t10-t00+t11-t01) +% dt/dy=(1/2h)*(t01-t00+t11-t10) +% approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences +% t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy +% u(1)=dt/dx +% u(2)=dt/dy +% u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy) + + tign_middle=(tign(1,1)+tign(1,2)+tign(2,1)+tign(2,2))/4; + dt_dx=(tign(2,1)-tign(1,1)+tign(2,2)-tign(1,2))/(2*fd(1)); + dt_dy=(tign(1,2)-tign(1,1)+tign(2,2)-tign(2,1))/(2*fd(2)); + u(1)=dt_dx; + u(2)=dt_dy; + u(3)=tign_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2; + + % 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) + %%%%% Approximation of the plane for lfn using finite differences +% approximate L(x,y)=u(1)*x+u(2)*y+u(3) + lfn_middle=(lfn(1,1)+lfn(1,2)+lfn(2,1)+lfn(2,2))/4; + dt_dx=(lfn(2,1)-lfn(1,1)+lfn(2,2)-lfn(1,2))/(2*fd(1)); + dt_dy=(lfn(1,2)-lfn(1,1)+lfn(2,2)-lfn(2,1))/(2*fd(2)); + u(1)=dt_dx; + u(2)=dt_dy; + u(3)=lfn_middle-(dt_dx*fd(1)+dt_dy*fd(2))/2; +% finding the coefficient c, reminder we work over one subcell only +% T(x,y)=c*L(x,y)+time_now + a=0; b=0; + for i=1:2 + for j=1:2 + if (lfn(i,j)<=0) + a=a+lfn(i,j)*lfn(i,j); + if (tign(i,j)>tnow) + disp('fuel_burnt_fd: tign(i1) should be less then time_now'); + else + b=b+(tign(i,j)-tnow)*lfn(i,j); + end + end + end + end + if (a==0) + display('fuel_burnt_fd: if c is on fire then one of cells should be on fire'); + end + c=b/a; + u(1)=u(1)*c; + u(2)=u(2)*c; + u(3)=u(3)*c+tnow; + u=u'; + 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 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)) + %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 + \ No newline at end of file diff --git a/standalone/matlab/fuel_test.m b/standalone/matlab/fuel_test.m index dfce9b67..2adf0694 100644 --- a/standalone/matlab/fuel_test.m +++ b/standalone/matlab/fuel_test.m @@ -26,7 +26,7 @@ for off=[0:0.05:1.1] % tign=tnow+2*lfn+0.0*randn(2,2); tign = [2.000000 1.565244 ; 2.000000 2.000000 ] - f=fuel_burnt(lfn,tign,tnow,fd,fuel_time); + f=fuel_burnt_fd(lfn,tign,tnow,fd,fuel_time); f_debug=[f_debug f]; tmp(k)=f; k=k+1; -- 2.11.4.GIT