standalone/matlab/prop_ls_ros.m

   1 function [t,phi]=prop_ls(phi0,time0,time1,vvx,vvy,r,dx,dy)\r
   2 % this version always propagates in the normal direction\r
   3 % the normal direction is approximated by gradient of level function\r
   4 % wind is projected into the rate of spread\r
   5 % but this does not work well for the windmill..\r
   6 \r
   7 % a jm/mkim version of propagate\r
   8 % follows Osher-Fedkiw and Mitchell toolbox\r
   9 \r
  10 % phi       level function out\r
  11 % phi0      level function in\r
  12 % time0     starting time\r
  13 % time1     end time\r
  14 % vx        component x of velocity field\r
  15 % vy        component y of velocity field\r
  16 % r         propagation speed in normal direction\r
  17 % dx        mesh step in x direction\r
  18 % dy        mesh step in y direction\r
  19 \r
  20 % initialize\r
  21 phi=phi0; % allocate level function\r
  22 [m,n]=size(phi);\r
  23 t=time0;  % current time\r
  24 tol=300*eps;\r
  25 istep=0;\r
  26 msteps=1000;\r
  27 diffLx=zeros(m,n);\r
  28 diffLy=zeros(m,n);\r
  29 diffRx=zeros(m,n);\r
  30 diffRy=zeros(m,n);\r
  31 \r
  32 % time loop\r
  33 while t<time1-tol & istep < msteps\r
  34     istep=istep+1;\r
  35     % one-sided differences\r
  36     [diffLx,diffRx,diffLy,diffRy,diffCx,diffCy]=get_diff(phi,dx,dy);     \r
  37     % Godunov scheme for normal motion\r
  38     flowLx=(diffLx>=0 & diffCx>=0);\r
  39     flowRx=(diffRx<=0 & diffCx<0);\r
  40     diff2x=diffLx.*flowLx + diffRx.*flowRx; %\r
  41     flowLy=(diffLy>=0 & diffCy>=0);\r
  42     flowRy=(diffRy<=0 & diffCy<0);\r
  43     diff2y=diffLy.*flowLy + diffRy.*flowRy; %\r
  44     grad=sqrt(diff2x.*diff2x + diff2y.*diff2y);\r
  45     nz=find(grad);\r
  46     %tbound_n(1)=max(abs(r.*sqrt(diff2x(nz)))./grad(nz))/dx; % time step bnd\r
  47     %tbound_n(2)=max(abs(r.*sqrt(diff2y(nz)))./grad(nz))/dy;\r
  48     %tbound_np=r.*(abs(diff2x)./dx+abs(diff2y)./dy); % pointwise time step bnd\r
  49     %tbound_n=max(tbound_np(nz)./abs(grad(nz))); % worst case\r
  50     %tend_n =-r.*grad;   % the result, tendency\r
  51     tbound_n=0; tend_n=0;\r
  52     % replace wind by normal advection\r
  53     \r
  54     % scale [diffCx,diffCy] to normal vector\r
  55     scale=sqrt(realmin+diffCx.*diffCx+diffCy.*diffCy);\r
  56     nvx = diffCx ./ scale;\r
  57     nvy = diffCy ./ scale;\r
  58 \r
  59     % given spread rate plus size of velocity field projected on the normal direction\r
  60     spread_rate=ros(r,vvx,vvy,nvx,nvy);\r
  61     % empirical correction\r
  62     spread_rate=max(0,spread_rate - 0.5*max(dx,dy)*scale);\r
  63 \r
  64     % get the advection size in the normal direction: \r
  65     % project [vx,vy] on the normal vector\r
  66     vx = nvx .* spread_rate;\r
  67     vy = nvy .* spread_rate;\r
  68     % standard upwinding for advection\r
  69     tend_a=-(diffLx.*max(vx,0)+diffRx.*min(vx,0)+...\r
  70         diffLy.*max(vy,0)+diffRy.*min(vy,0));\r
  71     tbound_a=max(abs(vx(:)))/dx+max(abs(vy(:)))/dy;\r
  72     % complete right-hand side\r
  73     tend=tend_n + tend_a;\r
  74     % complete time step bound\r
  75     tbound = 1/(tbound_n+tbound_a+realmin); \r
  76     % decide on timestep\r
  77     dt=min(time1-t,0.5*tbound);\r
  78     % trailing edge correction - do not allow fireline to go backwards\r
  79     %tt=max(dx,dy);\r
  80     %ins=find(phi<=0);\r
  81     %tend(ins)=min(tend(ins),-0.5);\r
  82     % tend=min(tend,0);\r
  83     % advance\r
  84     phi=phi+dt*tend;\r
  85     t=t+dt;\r
  86 end\r
  87 fprintf('prop_ls: %g steps from %g to %g\n',istep,time0,t)\r
  88 end % of the function prop_ls\r
  89 \r
  90 function spread_rate=ros(r,vx,vy,nvx,nvy)\r
  91 % get the spread rate from wind vector and\r
  92 % fireline normal vector\r
  93 spread_rate=max(r+vx.*nvx+vy.*nvy,0);\r
  94 end\r