standalone/matlab/fuel_test.m

   1 function fuel_frac=fuel_test\r
   2 for i=1:4,figure(i),end\r
   3 %input('Posdslkjflkjdfition 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 display('test1')\r
  11 lfn0=[-1    -1; 1    1] % to debug fortran\r
  12 tnow=2\r
  13 %lfn0=   [0.4192792     -1.9766893E-02;   6.310914       5.983462 ]\r
  14 fuel_time= 8.235294 ;\r
  15 f_debug=[];\r
  16 n=2^(9+1);\r
  17 \r
  18     tign = [1.000000      1.000000    ;\r
  19            2.000000       2.000000    ]\r
  20     f1=fuel_burnt(lfn0,tign,tnow,fd,fuel_time) \r
  21     display('fuel_burnt result')\r
  22     f2=fuel_burnt_fd(lfn0,tign,tnow,fd,fuel_time) \r
  23     display('fuel_burnt_fd result')\r
  24     fq=fuel_burnt_quad(lfn0,tign,tnow,fd,fuel_time,n)\r
  25     display('fuel_burnt_quad result')\r
  26     test1_ls=f1-fq       \r
  27     display('fuel_burnt error')\r
  28     test1_fd=f2-fq       \r
  29     display('fuel_burnt error')\r
  30 \r
  31     display('  ')\r
  32     \r
  33     display('test2')\r
  34 \r
  35 \r
  36 lfn0=[-2    -2; -1    -1] % to debug fortran\r
  37 tnow=3\r
  38 %lfn0=   [0.4192792     -1.9766893E-02;   6.310914       5.983462 ]\r
  39 fuel_time= 8.235294 ;\r
  40 f_debug=[];\r
  41 n=2^(9+1);\r
  42 \r
  43     tign = [1.000000      1.000000    ;\r
  44            2.000000       2.000000    ]\r
  45     f1=fuel_burnt(lfn0,tign,tnow,fd,fuel_time) \r
  46     display('fuel_burnt result')\r
  47     f2=fuel_burnt_fd(lfn0,tign,tnow,fd,fuel_time) \r
  48     display('fuel_burnt_fd result')\r
  49     fq=fuel_burnt_quad(lfn0,tign,tnow,fd,fuel_time,n)\r
  50     display('fuel_burnt_quad result')\r
  51     test2_ls=f1-fq       \r
  52     display('fuel_burnt error')\r
  53     test2_fd=f2-fq       \r
  54     display('fuel_burnt error')\r
  55 \r
  56     \r
  57 display('test 3')\r
  58 \r
  59 lfn0=[-2    -2; -4    -4] % to debug fortran\r
  60 tnow=10\r
  61 %lfn0=   [0.4192792     -1.9766893E-02;   6.310914       5.983462 ]\r
  62 fuel_time= 8.235294 ;\r
  63 f_debug=[];\r
  64 n=2^(9+1);\r
  65 \r
  66     tign = [6.000000      6.000000    ;\r
  67            4.000000       4.000000    ]\r
  68     f1=fuel_burnt(lfn0,tign,tnow,fd,fuel_time) \r
  69     display('fuel_burnt result')\r
  70     f2=fuel_burnt_fd(lfn0,tign,tnow,fd,fuel_time) \r
  71     display('fuel_burnt_fd result')\r
  72     fq=fuel_burnt_quad(lfn0,tign,tnow,fd,fuel_time,n)\r
  73     display('fuel_burnt_quad result')\r
  74     test3_ls=f1-fq       \r
  75     display('fuel_burnt error')\r
  76     test3_fd=f2-fq       \r
  77     display('fuel_burnt error')\r
  78 \r
  79     \r
  80     \r
  81     \r
  82 lfn0=[-3    -1; -1    1] % to debug fortran\r
  83 tnow=10\r
  84 %lfn0=   [0.4192792     -1.9766893E-02;   6.310914       5.983462 ]\r
  85 fuel_time= 8.235294 ;\r
  86 f_debug=[];\r
  87 n=2^(9+1);\r
  88 \r
  89     tign = [4.000000      8.000000    ;\r
  90            8.000000       10.000000    ]\r
  91     f1=fuel_burnt(lfn0,tign,tnow,fd,fuel_time) \r
  92     display('fuel_burnt result')\r
  93     f2=fuel_burnt_fd(lfn0,tign,tnow,fd,fuel_time) \r
  94     display('fuel_burnt_fd result')\r
  95     fq=fuel_burnt_quad(lfn0,tign,tnow,fd,fuel_time,n)\r
  96     display('fuel_burnt_quad result')\r
  97     test4_ls=f1-fq       \r
  98     display('fuel_burnt error')\r
  99     test4_fd=f2-fq       \r
 100     display('fuel_burnt error')\r
 101 \r
 102 lfn0=[-1    1; 1    3] % to debug fortran\r
 103 tnow=2\r
 104 %lfn0=   [0.4192792     -1.9766893E-02;   6.310914       5.983462 ]\r
 105 fuel_time= 8.235294 ;\r
 106 f_debug=[];\r
 107 n=2^(9+1);\r
 108 \r
 109     tign = [1.000000      2.000000    ;\r
 110            2.000000       2.000000    ]\r
 111     f1=fuel_burnt(lfn0,tign,tnow,fd,fuel_time) \r
 112     display('fuel_burnt result')\r
 113     f2=fuel_burnt_fd(lfn0,tign,tnow,fd,fuel_time) \r
 114     display('fuel_burnt_fd result')\r
 115     fq=fuel_burnt_quad(lfn0,tign,tnow,fd,fuel_time,n)\r
 116     display('fuel_burnt_quad result')\r
 117     test5_ls=f1-fq       \r
 118     display('fuel_burnt error')\r
 119     test5_fd=f2-fq       \r
 120     display('fuel_burnt error')\r
 121     test1_ls\r
 122     test1_fd\r
 123    \r
 124     test2_ls\r
 125     test2_fd\r
 126     \r
 127     test3_ls\r
 128     test3_fd\r
 129     \r
 130     test4_ls\r
 131     test4_fd\r
 132     \r
 133     test5_ls\r
 134     test5_fd\r
 135    \r
 136     \r
 137     \r
 138     \r
 139 end\r
 140 \r
 141 \r
 142 function fuel_frac=fuel_burnt_quad(lfn,tign,tnow,fd,fuel_time,n)\r
 143 % compute the same numerically for comparison\r
 144 % if fireline passes through the cell then\r
 145 % lfn and tign must be linear\r
 146 %otherwise the result will differ\r
 147 % evaluate by quadrature on n by n points\r
 148 % the mesh\r
 149 g=(0.5+[0:n-1])/n;\r
 150 [q1,q2]=meshgrid(fd(1)*g,fd(2)*g); % quad nodes\r
 151 [c1,c2]=meshgrid([0,fd(1)],[0,fd(2)]); % corners of the mesh cell\r
 152 % T and L on quad nodes\r
 153 T = interp2(c1,c2,tign,q1,q2)-tnow;\r
 154 L = interp2(c1,c2,lfn,q1,q2);\r
 155 % the integrand\r
 156 f = (1 - exp((L<0).*T./fuel_time)).*(L<0);\r
 157 % integrate\r
 158 fuel_frac = sum(f(:))/(n*n);\r
 159 end\r
 160 \r