standalone/matlab/prop_ls_cir.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 istep=0;\r
  22 msteps=1000;\r
  23 diffLx=zeros(m,n);\r
  24 diffLy=zeros(m,n);\r
  25 diffRx=zeros(m,n);\r
  26 diffRy=zeros(m,n);\r
  27 \r
  28 method='godunov'\r
  29 split='spread';\r
  30 % split='wind';\r
  31 % split = 'orig';\r
  32 \r
  33 % time loop\r
  34 while t<time1-tol & istep < msteps\r
  35     istep=istep+1;\r
  36     % one-sided differences\r
  37     [diffLx,diffRx,diffLy,diffRy,diffCx,diffCy]=get_diff(phi,dx,dy);         \r
  38     % gradient of level function is the normal direction\r
  39     scale=sqrt(diffCx.*diffCx+diffCy.*diffCy+eps); \r
  40     nvx=diffCx./scale;\r
  41     nvy=diffCy./scale;\r
  42     % propagation speed \r
  43     %speed=spread_rate(r,vx,vy,nvx,nvy,scale);\r
  44     %\r
  45     speed=vx.*nvx + vy.*nvy;\r
  46     speed=max(speed,0);\r
  47     % to recover advection vv and spread r, transition between:\r
  48     % r = 0 & (normal,v) > const*speed => vv=speed*v/(normal,v) & rr=0\r
  49     % speed >> (normal,v) => vv=0 & rr=speed    \r
  50     nvv = nvx .* vx + nvy .* vy;\r
  51     % the CIR method\r
  52     switch method\r
  53         case 'cir'\r
  54             dt=min(min(dx./max(abs(speed.*nvx),eps)))+min(min(dy./max(abs(speed.*nvy),eps)));\r
  55             dt=min(time1-t,0.5*dt)\r
  56             % locations back on characteristics, in mesh index coordinates\r
  57             [xii,yii]=meshgrid(1:m,1:n);\r
  58             xinc=dt.*speed.*nvx./dx;\r
  59             xi = xii - xinc;\r
  60             yinc=dt.*speed.*nvy./dy;\r
  61             yi = yii - yinc; \r
  62             phi_new=interp2(xii,yii,phi,xi,yi);\r
  63             phi=min(phi_new,phi);\r
  64         case 'godunov'\r
  65             switch split\r
  66             case 'wind'\r
  67                 a=2*nvv>speed;\r
  68                 nvv(a==0)=1;\r
  69                 ! a=zeros(size(a));\r
  70                 %no_wind_cases=sum(sum(1-a))\r
  71                 rr=speed.*(1-a);\r
  72                 corr = a .* speed ./ nvv\r
  73                 vvx = vx .* corr;\r
  74                 vvy = vy .* corr;\r
  75             case 'orig'\r
  76                 vvx = vx;vvy= vy;rr=r;  % all original, if r=0 pure upwinding\r
  77             case 'spread'\r
  78                 rr = speed; vvx=0*vx; vvy=0*vy; % all in normal speed, Godunov meth.\r
  79                 % vvx=vx;vvy=vy;\r
  80             otherwise\r
  81                 error('unknown split')\r
  82             end\r
  83             % Godunov scheme for normal motion\r
  84             flowLx=(diffLx>=0 & diffCx>=0);\r
  85             flowRx=(diffRx<=0 & diffCx<0);\r
  86             diff2x=diffLx.*flowLx + diffRx.*flowRx; %\r
  87             flowLy=(diffLy>=0 & diffCy>=0);\r
  88             flowRy=(diffRy<=0 & diffCy<0);\r
  89             diff2y=diffLy.*flowLy + diffRy.*flowRy; %\r
  90             grad=sqrt(diff2x.*diff2x + diff2y.*diff2y);\r
  91             nz=find(grad);\r
  92             %tbound_n(1)=max(abs(r.*sqrt(diff2x(nz)))./grad(nz))/dx; % time step bnd\r
  93             %tbound_n(2)=max(abs(r.*sqrt(diff2y(nz)))./grad(nz))/dy;\r
  94             tbound_np=rr.*(abs(diff2x)./dx+abs(diff2y)./dy); % pointwise time step bnd\r
  95             tbound_n=max(tbound_np(nz)./abs(grad(nz))); % worst case\r
  96             tend_n =-rr.*grad;   % the result, tendency\r
  97             % standard upwinding for advection\r
  98             tend_a=-(diffLx.*max(vvx,0)+diffRx.*min(vvx,0)+...\r
  99             diffLy.*max(vvy,0)+diffRy.*min(vvy,0));\r
 100             tbound_ax=max(abs(vvx(:)))/dx;\r
 101             tbound_ay=max(abs(vvy(:)))/dy;\r
 102             % complete right-hand side\r
 103             tend=tend_n + tend_a;\r
 104             % complete time step bound\r
 105             tbound = 1/(tbound_n+tbound_ax+tbound_ay+eps); \r
 106             % decide on timestep\r
 107             dt=min(time1-t,0.5*tbound);\r
 108             % trailing edge correction - do not allow fireline to go backwards\r
 109             tend=min(tend,0);\r
 110             % advance\r
 111             phi=phi+dt*tend;\r
 112     end\r
 113     t=t+dt;\r
 114 end\r
 115 fprintf('prop_ls: %g steps from %g to %g, split %s\n',istep,time0,t,split)\r
 116 end % of the function prop_ls\r
 117 \r