1 /*
2 * $Id: test2.i,v 1.1 2005-09-18 22:06:11 dhmunro Exp $
3 * Highly vectorizing physics problem for timing Yorick.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
test2(npass)11 func test2(npass)
12 /* DOCUMENT test2
13 or test2, npass
14 Given a slab divided into zones by parallel planes, and given a
15 set of photon group boundary energies, compute the contribution of
16 each zone to the radiation flux emerging from one surface of the
17 slab. If NPASS is given, the calculation is repeated that many
18 times. The zoning, photon group structure, opacities, and source
19 functions are all computed arbitrarily, but the number of zones
20 and groups are taken to be representative of a typical 1-D
21 radiation transport calculation. */
22 {
23 if (is_void(npass) || npass<=0) npass= 1;
24
25 now= split= exponential= array(0.0, 3);
26 timer, now;
27 n= npass;
28 while (n--) esc= escout(zt, akap, srcfun);
29 timer, now, split;
30 timer_print, "Time per pass", split/npass, "Total time", split,\
31 "Computing exponentials", exponential;
32
33 if (esc(0,mxx)!=29 || !approx_eq(esc(0,max),0.001875213417) ||
34 !approx_eq(esc(0,1),0.001694423767) ||
35 !approx_eq(esc(0,0),5.434635527e-05) ||
36 esc(1,mxx)!=39 || !approx_eq(esc(1,max),0.000440460566))
37 write, "***WARNING*** values returned by escout are not correct";
38
39 if (flxout(0,mxx)!=34 || !approx_eq(flxout(0,max),0.06615472064) ||
40 !approx_eq(flxout(0,1),0.003187516911) ||
41 !approx_eq(flxout(0,0),0.005280842058) ||
42 flxout(1,mxx)!=29 || !approx_eq(flxout(1,max),0.001805157164))
43 write, "***WARNING*** values of flxout are not correct";
44
45 if (!approx_eq(min(tau(:-1,)),6.982090961e-05) ||
46 !approx_eq(max(tau),45.80160946))
47 write, "***WARNING*** values of tau are not correct";
48 }
49
approx_eq(x,y)50 func approx_eq(x, y)
51 {
52 return (abs(x-y)/(abs(x+y)+1.e-35))<1.e-6;
53 }
54
55 #if 0
56 CPU seconds per pass:
57 IDL Yorick Basis DAP FORTRAN(-O)
58 HP730 - 0.60 2.04 1.84 0.27 (0.60 -g)
59 Solbourne 2.81 1.90 6.02 5.90 1.00
60 (varies by ~10%)
61
62 Using forum (Solbourne) on 8/Dec/92:
63 forum[18] time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
64 IDL. Version 2.3.0 (sunos sparc).
65 Copyright 1989-1992, Research Systems, Inc.
66 All rights reserved. Unauthorized reproduction prohibited.
67 Site: 1491.
68 Licensed for use by: LLNL - X Division
69 IDL> .rnew /home/miggle/munro/Yorick/include/test2
70 % Compiled module: ESCOUT.
71 % Compiled module: APPROX_EQ.
72 % Compiled module: TEST2.
73 % Compiled module: BNU.
74 % Compiled module: OPACSET.
75 % Compiled module: SPAN.
76 % Compiled module: SPANL.
77 % Compiled module: $MAIN$.
78 IDL> test2,1
79 IDL> exit
80 3.2u 0.4s 0:12 29% 0+2496k 1+2io 1pf+0w
81 forum[19] time time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
82 IDL. Version 2.3.0 (sunos sparc).
83 Copyright 1989-1992, Research Systems, Inc.
84 All rights reserved. Unauthorized reproduction prohibited.
85 Site: 1491.
86 Licensed for use by: LLNL - X Division
87 IDL> .rnew /home/miggle/munro/Yorick/include/test2
88 % Compiled module: ESCOUT.
89 % Compiled module: APPROX_EQ.
90 % Compiled module: TEST2.
91 % Compiled module: BNU.
92 % Compiled module: OPACSET.
93 % Compiled module: SPAN.
94 % Compiled module: SPANL.
95 % Compiled module: $MAIN$.
96 IDL> test2,11
97 IDL> exit
98 31.3u 0.9s 0:58 54% 0+3544k 0+0io 0pf+0w
99
100
101 forum[20] time yorick
102 Yorick ready. For help type 'help'
103 > #include "/home/miggle/munro/Yorick/include/test2.i"
104 > test2,1
105 Timing Category CPU sec System sec Wall sec
106 Time per pass 1.840 0.300 2.270
107 Total time 1.840 0.300 2.270
108 Computing exponentials 0.600 0.120 0.760
109 -----Total Elapsed Times----- 2.510 0.490 27.500
110 > test2,10
111 Timing Category CPU sec System sec Wall sec
112 Time per pass 1.932 0.050 1.986
113 Total time 19.320 0.500 19.860
114 Computing exponentials 6.360 0.080 6.440
115 -----Total Elapsed Times----- 21.860 1.020 53.080
116 > quit
117 21.9u 1.0s 0:59 38% 0+3560k 31+0io 85pf+0w
118
119
120 forum[21] basis
121 Basis (basis, Version 921125)
122 Run at 13:02:49 on 12/08/92 on the forum machine, suffix 10797x
123 Initializing Basis System
124 Basis 7.0
125 Initializing EZCURVE/NCAR Graphics
126 Ezcurve/NCAR 2.0
127 Basis> echo=0
128 Basis> read /home/miggle/munro/Yorick/include/test2.bas
129 End of input file /home/miggle/munro/Yorick/include/test2.bas
130 Resuming input from TERMINAL
131 Basis> test2(1)
132
133
134 CPU (sec) SYS (sec)
135 6.017 0.300
136 Basis> test2(10)
137
138 CPU (sec) SYS (sec)
139 60.233 0.017
140 Basis> end
141
142 CPU (sec) SYS (sec)
143 68.617 0.917
144 68.6u 0.9s 2:07 54% 0+6144k 19+1io 213pf+0w
145
146
147 forum[22] dap
148 DAP> read /home/miggle/munro/Yorick/include/test2.bas
149 DAP> test2(1)
150
151 CPU (sec) SYS (sec)
152 5.883 0.400
153 DAP> test2(10)
154
155 CPU (sec) SYS (sec)
156 59.533 0.000
157 DAP> end
158 69.1u 1.5s 1:47 65% 0+6488k 20+1io 218pf+0w
159
160
161 Yorick results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
162 top showed SIZE/RES 764K/244K at prompt
163 3636K/3052K after test2,1
164 4628K/4012K after test2,20
165 tonto[9] yorick
166 Yorick ready. For help type 'help'
167 > #include "/home/miggle/munro/Yorick/include/test2.i"
168 > test2,1
169 Timing Category CPU sec System sec Wall sec
170 Time per pass 0.550 0.100 0.670
171 Total time 0.550 0.100 0.670
172 Computing exponentials 0.240 0.030 0.270
173 -----Total Elapsed Times----- 0.930 0.220 34.950
174 > test2,20
175 Timing Category CPU sec System sec Wall sec
176 Time per pass 0.599 0.001 0.906
177 Total time 11.980 0.020 18.110
178 Computing exponentials 5.030 0.000 7.270
179 -----Total Elapsed Times----- 12.920 0.240 83.120
180
181 Basis results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
182 top showed SIZE/RES 5356K/540K at prompt
183 8960K/4108K after test2(1)
184 8960K/4112K after test2(20)
185 tonto[8] basis
186 Basis (basis, Version 921125)
187 Run at 21:45:42 on 12/04/92 on the tonto machine, suffix 27713x
188 Initializing Basis System
189 Basis 7.0
190 Initializing EZCURVE/NCAR Graphics
191 Ezcurve/NCAR 2.0
192 Basis> echo=0
193 Basis> read /home/miggle/munro/Yorick/include/test2.bas
194 End of input file /home/miggle/munro/Yorick/include/test2.bas
195 Resuming input from TERMINAL
196 Basis> test2(1)
197
198 CPU (sec) SYS (sec)
199 1.990 .120
200 Basis> test2(20)
201
202 CPU (sec) SYS (sec)
203 40.800 .000
204
205 DAP results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
206 top showed SIZE/RES 6796K/844K at prompt
207 10360K/4352K after test2(1)
208 10360K/4352K after test2(20)
209 tonto[11] dap
210 DAP> echo=0
211 DAP> read /home/miggle/munro/Yorick/include/test2.bas
212 DAP> test2(1)
213
214 CPU (sec) SYS (sec)
215 1.840 .080
216 DAP> test2(20)
217
218 CPU (sec) SYS (sec)
219 36.830 .000
220 #endif
221
222 /* This routine computes the optical depth through each zone at every
223 frequency and then uses that to compute the radiation emitted in
224 each zone that escapes from the problem, assuming plane parallel
225 geometry.
226 The returned result is an nzones-by-ngroups array of
227 power per unit photon energy per unit area. */
escout(zt,akap,srcfun)228 func escout(zt, /* npoints zone boundary positions (cm) */
229 akap, /* nzones-by-ngroups opacities (1/cm) */
230 srcfun) /* nzones-by-ngroups source (specific intensity units) */
231 {
232 extern flxout; /* (output) nzones-by-ngroups outgoing fluxes */
233 extern dtau; /* (output) nzones-by-ngroups optical depths (ODs) */
234 extern tau; /* (output) npoints-by-ngroups cumulative ODs */
235
236 extern mu, wmu; /* Gauss-Legendre cos(theta) and weight arrays
237 for integration over escape angles */
238
239 /* compute tau, the optical depth to each zone along the zt-direction */
240 dtau= akap*zt(dif);
241 tau= array(0.0, dimsof(zt), dimsof(akap)(3));
242 tau(2:,)= dtau(psum,);
243 /* consider the outward going radiation, and thus use tau
244 measured from the right boundary */
245 tau= tau(0,)(-,) - tau;
246
247 /* compute exf, the fraction of the srcfun which escapes from
248 each zone in each bin at each angle mu relative to the zt direction */
249 enow= array(0.0, 3);
250 timer, enow; /* most of the actual work is computing these exponentials */
251 exf= exp(-tau(,,-)/mu(-,-,));
252 timer, enow, exponential;
253
254 /* compute the escaping flux per unit surface area in each bin
255 contributed by each zone -- units are 10^17 W/kev/cm^2 */
256 /* Note:
257 dI/dtau = -I + S (S is srcfun, which is Bnu) implies
258 I = integral of ( S * exp(tau) ) dtau from tau -infinity to tau 0
259 Hence, the contribution of a single zone to this integral is
260 S12 * (exp(tau1) - exp(tau2)), which explains the exf(dif,,) below.
261 The wmu are Gauss-Legendre integration weights, and the mu is the
262 cos(theta) to project the specific intensity onto the direction
263 normal to the surface (zt). */
264 esfun= 2*pi*exf(dif,,)*srcfun(,,-)*(mu*wmu)(-,-,);
265
266 /* Also compute what the total flux WOULD have been if each successive
267 zone were at the surface. This is the same as the "one-sided"
268 flux directed outward across each zone boundary. */
269 fuzz= 1.0e-10;
270 flxout= (esfun(psum,,)/(exf(2:,,)+fuzz))(,,sum);
271
272 return esfun(,,sum);
273 }
274
275 /* This function is used to set a dummy opacity to be used by the
276 transport calculation */
opacset(tmp,rho,gav,gb)277 func opacset(tmp, rho, gav, gb)
278 {
279 extern srcfun; /* (output) nzones-by-ngroups Bnu, LTE source function */
280 extern akap; /* (output) nzones-by-ngroups opacities (1/cm) */
281
282 /* the opacity is proportional to the density to the rhopow
283 power and the temperature to the tempow power */
284 rhopow= 2;
285 tempow= 3;
286 factr= (tmp/1.0)^tempow*(rho/1.0)^rhopow;
287
288 /* set the constants in front of the terms for the two edges */
289 con0= 1.0e+4;
290 cona= 2.5e+2;
291 conb= 5.0e+0;
292 /* set the energies of the two edges */
293 edg0= 0.1;
294 edga= 0.5;
295 edgb= 2.0;
296
297 /* set up arrays that are zero below the edges and one above them */
298 vala= double(gav>edga);
299 valb= double(gav>edgb);
300
301 /* frequency dependence is the same for all zones, so calculate it once */
302 frval= con0*(edg0/gav)^3+cona*vala*(edga/gav)^3+conb*valb*(edgb/gav)^3;
303
304 /* set the opacity */
305 akap= factr(,-)*frval(-,);
306 /* the source function is always a blackbody */
307 srcfun= bnu(tmp, gb);
308 }
309
310 /* Make a blackbody using the analytic formula --
311 the units are 10^17 W/kev/g/ster */
bnu(tmp,freqb)312 func bnu(tmp, freqb)
313 {
314 /* Note that tmp and freqb are one dimensional and the output
315 array should be 2D */
316 /* compute the derivative using the exact black body,
317 not the average over the bin */
318 exf= exp(-(freqb(zcen)(-,)/tmp));
319 return 0.0504*(freqb(zcen)^3)(-,)*exf/(1.0-exf);
320 }
321
322 /* set up default params. etc. */
323
324 /* set the number of spatial zones and photon groups */
325 npoints= 100;
326 ngroups= 64;
327 nzones= npoints-1;
328
329 /* set the zone boundary positions */
330 zt= span(0.0, 0.0050, npoints);
331
332 /* set the frequency bins */
333 gb= spanl(0.1, 4.0, ngroups+1);
334 gav= gb(zcen);
335
336 /* set the density and temperature */
337 rho= spanl(1.0, 1.0, nzones);
338 tmp= spanl(1.0, 1.0, nzones);
339
340 /* set the opacity akap and source function srcfun */
341 opacset, tmp, rho, gav, gb;
342
343 /* set up for a Gauss-Legendre angular quadrature */
344 xmu= [0.1488743389, 0.4333953941, 0.6794095682, 0.8650633666, 0.9739065285];
345 mu= -xmu(::-1);
346 grow, mu, xmu;
347 xmu= [0.2955242247, 0.2692667193, 0.2190863625, 0.1494513491, 0.0666713443];
348 wmu= xmu(::-1);
349 grow, wmu, xmu;
350 /* correct the Gauss-Legendre abcissas-- interval is only 0 to 1 */
351 mu= 0.5 + 0.5*mu;
352