wrfv2_fire/phys/prop_ls.m

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