wrfv2_fire/phys/fuel_test.m

   1 function fuel_frac=fuel_test \r
   2 for i=1:4,figure(i),end\r
   3 input('Position figure windows and press Enter >');\r
   4 fd=[1,1];\r
   5 [c1,c2]=ndgrid([0,fd(1)],[0,fd(2)]); % corners of the mesh cell\r
   6 maxerr=0;\r
   7 \r
   8 %tnow=1;\r
   9 %fuel_time=1;\r
  10 %lfn0=[0    0.3479;0.6521    1.0000]; % to debug fortran\r
  11 tnow=2\r
  12 lfn0=   [0.4192792     -1.9766893E-02    ;\r
  13          6.310914       5.983462 ]\r
  14 fuel_time= 8.235294 ;\r
  15 f_debug=[];\r
  16 \r
  17 %!./ifmake clean\r
  18 %!./ifmake fuel_burnt_test\r
  19 \r
  20 tmp=zeros(23,1);\r
  21 k=1;\r
  22 for off=[0:0.05:1.1]\r
  23     if off >= 1.0 \r
  24         disp(off)\r
  25     end    \r
  26     lfn=lfn0-off;\r
  27 %    tign=tnow+2*lfn+0.0*randn(2,2);\r
  28     tign = [2.000000      1.565244    ;\r
  29            2.000000       2.000000    ]\r
  30     f=fuel_burnt(lfn,tign,tnow,fd,fuel_time); \r
  31     f_debug=[f_debug f];\r
  32     tmp(k)=f;\r
  33     k=k+1;\r
  34     f_debug=[f_debug f];\r
  35     err=[];\r
  36     for i=[1:9]\r
  37        n=2^(i+1);\r
  38        fq=fuel_burnt_quad(lfn,tign,tnow,fd,fuel_time,n);\r
  39         nn(i)=n;\r
  40         err(i)=abs(f-fq);\r
  41         last=max(1,i-1);\r
  42         figure(4)\r
  43         loglog(nn(last:i),err(last:i),'-ok')\r
  44         xlabel n,ylabel difference,title('Comparison with numerical quadrature')\r
  45         axis([3 1100 1e-8 1]),grid on, hold on\r
  46     end\r
  47      hold off\r
  48 end\r
  49 \r
  50 fid=fopen('tmp.txt','r');\r
  51 A= fscanf(fid,'%g'); \r
  52 fclose(fid);\r
  53 err_mat_fortran=abs(A-tmp)\r
  54 save('err.mat','err_mat_fortran')\r
  55 disp(f_debug);\r
  56 end\r
  57 \r
  58 \r
  59 function fuel_frac=fuel_burnt_quad(lfn,tign,tnow,fd,fuel_time,n)\r
  60 % compute the same numerically for comparison\r
  61 % if fireline passes through the cell then \r
  62 % lfn and tign must be linear\r
  63 %otherwise the result will differ\r
  64 % evaluate by quadrature on n by n points\r
  65 % the mesh\r
  66 g=(0.5+[0:n-1])/n;\r
  67 [q1,q2]=meshgrid(fd(1)*g,fd(2)*g); % quad nodes\r
  68 [c1,c2]=meshgrid([0,fd(1)],[0,fd(2)]); % corners of the mesh cell\r
  69 % T and L on quad nodes\r
  70 T = interp2(c1,c2,tign,q1,q2)-tnow;\r
  71 L = interp2(c1,c2,lfn,q1,q2);\r
  72 % the integrand\r
  73 f = (1 - exp((L<0).*T./fuel_time)).*(L<0);\r
  74 % integrate\r
  75 fuel_frac = sum(f(:))/(n*n);\r
  76 end\r