Nuclear Engineering at NC State

In order to overcome the difficulty in using current MCNP[1] weight window technique, a new weight window generator and matching weight window game were developed and incorporated into MCNP 4A through a patch file, and the patch file were later updated to match the latest MCNP version -- MCNP 4B. The new weight window method differs from MCNP weight window in that it uses a fine spatial mesh system, which is independent of the physical geometry, to describe thespatial dependence of the importance function, rather than physical cells as used in MCNP weight window. The de-coupling of importance mesh structure and the physical geometry results in user-friendliness.

It has been demonstrated that the new weight window method is more efficient than current MCNP cell-based weight window. This efficiency gain can be attributed to the better spatial resolution and less overhead in computing intersections between the particle tracks and mesh surfaces. Besides the user-friendliness and the high efficiency, the new weight window also provides an easy way to generate global distribution information, which would otherwise be very difficult, if not impossible, to obtain, in ordinary Monte Carlo calculations. With this version of the patch file, three kinds of global functions are provided; they are regular (forward) particle flux into each mesh, contributon flux into each mesh, and the adjoint flux into each mesh. Among these, the contributon current distribution, which shows how particles from the source travel inside medium and reach the detector, provides the physical insight for analysis. More details about this mesh-based weight window can be found elsewhere[2][3].

The current version of the patch file can support neutron-only, gamma-ray-only and coupled neutron gamma-ray transport. Also time- dependency can be included in the importance map function. Extensive tests have been conducted for neutron-only and gamma-ray-only cases. Further tests for neutron gamma-ray coupled and time-dependent cases are underway. Users are encouraged to report any bugs or to send suggestions to the author by e-mail for further improvement.

2. Input File Preparation

The input file for the modified MCNP 4B is the same as for standard MCNP, except for the presence of the two user-specified input cards (idum and rdum), which are provided by MCNP for non- standard input control as explained in MCNP manual[1]. In addition to these two cards, three standard MCNP input cards (wwg, wwge:n and wwge:p) are also used for the fine-mesh-based weight window.

The first two cards are used to specify the spatial mesh (X-Y-Z) structure and time bins for weight window map. "idum" is to be followed by integer entries and "rdum" by float entries. The maximum entries for each card are 50.

The meanings of these two cards are explained below.

idum(1) -- weight window control (kww). 0 -- the new WW is not invoked. Normal MCNP option. 1 -- turn on the fine-mesh weight window game without WW map generation. (Input WW map file "imp.in' is required) 2 -- play WW game when WW map is available and turn on WW generator. No WW map resetting. (use with care) 3 -- play WW game when WW map is available and turn on WW generator. WW map is reset after each MCNP print-out. (recommended over 2) idum(2) -- number of uniform regions for X (nxr=abs(idum(2)). idum(3),..., idum(2+nxr) -- mesh number in each of the nxr regions. idum(3+nxr) -- number of uniform regions for Y (nyr=abs(idum(3+nxr)). idum(4+nxr),..., idum(3+nxr+nyr) -- mesh number in each of the nyr regions. idum(4+nxr+nyr) -- number of uniform regions for Z (nzr=abs(idum(4+nxr+nyr)). idum(5+nxr+nyr),..., idum(4+nxr+nyr+nzr) -- mesh number in each of the nzr regions. idum(5+nxr+nyr+nzr) -- number of uniform regions for time (ntr=abs(idum(5+nxr+nyr+nzr)). idum(6+nxr+nyr+nzr),..., idum(5+nxr+nyr+nzr+ntr) -- mesh number in each of the ntr regions. rdum(1),..., rdum(1+nxr) -- region boundaries for each of the uniform regions for X. rdum(2+nxr),..., rdum(2+nxr+nyr) -- region boundaries for each of the uniform regions for Y. rdum(3+nxr+nyr),..., rdum(3+nxr+nyr+nzr) -- region boundaries for each of the uniform regions for Z. rdum(4+nxr+nyr+nzr),..., rdum(4+nxr+nyr+nzr+ntr) -- region boundaries for each of the uniform regions for time.

Standard MCNP wwg card is used to choose the tally number for which the importance map is optimized. The wwg card has seven entries,

wwg It Ic Wg Xr Yr Zr J

In this application, only the first entry has effect.

    The weight window generation energies card wwge works the same way 
as for standard MCNP.



wwge:n E1 E2... Ei... Ej; j<16



where Ei is the upper energy bound for weight window group and n is the 
particle type indicator. Only neutron (n=n) and photon (n=p) types are 
allowed here.



    When the entered number of uniform regions for a dimension is negative, 
it indicates a symmetry with respect to the leftmost plane of the 
boundaries and only the meshes on the right side of the symmetry plane 
are specified by input. The meshes on the left is generated 
automatically according to the symmetry and the actual mesh number of 
the dimension is doubled. Using geometric symmetry may accelerate the 
convergence of WW map.



CAUTION



    Special care needs to be taken to assure of no overlapping between the 
non-physical dividing planes and the planes used for physical geometric
description.



    In the current patch file, the maximum mesh number allowed is 6 
(maxne=6) energy groups, 30 (maxnx=30) X meshes, 30 (maxny=30) Y meshes,
30 (maxnz=30) Z meshes and 4 (maxnt=30) time meshes. If a mesh structure
more than this maximum is desired, the following statement appearing in 
the patch file needs to be modified and MCNP pre-processing has to be 
redone.



      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)



The mesh structure is printed in standard MCNP output file (outp by 
default) and its verification is recommended. 



    Future development of this fine mesh weight window will replace the 
idum and rdum cards with wwgx, wwgy, wwgz and wwgt cards, to be 
consistent with standard MCNP card wwge.



 
3. WW Map File Format


    There are two binary WW map files involved when running this special 
version of MCNP, i.e. imp.in and imp.out. The first one is read in when 
required, and the second is written out. These files can contain very 
useful global distribution information for physical analysis. 



    The data in the files can be read out by the following FORTRAN 
statements.




      logical x_symm,y_symm,z_symm
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      dimension rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $          cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $          emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $          zmesh(maxnz+1),tmesh(maxnt+1),nemesh(2)
c
      open(1,file='imp.in',status='old',form='unformatted')
      read(1)x_symm,y_symm,z_symm
      read(1)nxmesh,nymesh,nzmesh,nemesh(1),nemesh(2),ntmesh
      read(1)(xmesh(i1),  i1=1,nxmesh+1)
      read(1)(ymesh(i2),  i2=1,nymesh+1)
      read(1)(zmesh(i3),  i3=1,nzmesh+1)
      read(1)(emesh(i4,1),i4=1,nemesh(1)+1)
      read(1)(emesh(i4,2),i4=1,nemesh(2)+1) 
      read(1)(tmesh(i5),  i5=1,ntmesh+1)
      read(1)frmsrc,rspns
      read(1)(((((rflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $         i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)

      read(1)(((((cflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $         i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)
      close(1)



where the array rflux corresponds to the regular flux and cflux to the 
contributon flux. The mesh structure is defined by arrays, xmesh
(spatial--X), ymesh(spatial--Y), zmesh (spatial--Z), emesh(*,1)

(neutron energy), emesh(*,2) (gamma-ray energy), and tmesh (time).






4. Installation of the Special MCNP 4B




    The installation process is the same as for standard MCNP, except the 
patch file used when pre-processing with PRPR program. For standard MCNP 
pre-processing, the following commands are executed on UNIX workstation.



> cp patchf patch
> cp mcnpf.id codef
> prpr



Here the patch file (patchf) is simply a machine specifier. When preparing 
the special version of MCNP, the following commands need to be executed 
instead. 



> cat patchf patch.cear > patch
> cp mcnpf.id codef
> prpr



The difference is the first command line, by which the patch file 
patch.cear is appended to the standard patch file.





5. Example 1 -- Neutron Porosity Benchmark




    The first example is a neutron porosity benchmark problem. The detailed 
specification can be found in Reference [4]. In both of these examples 
presented below, history number control is used to terminate the MCNP run, 
so as to make this simulation results repeatable. In real calculations, CPU 
time control may be used instead. 



    The input file for the second example is shown below 
(filename = test1.i).



NCSU neutron porosity logging tool -- limestone of 20 PU 
1     1   -0.000502   -1   7   -8                      $ near detector
2     1   -0.000502   -2   9   -10                     $ far  detector
3     2   -7.8        -3   6   -11  #1  #2             $ tool
4     3   -1.0        -4   3    6   -11                $ borehole
5     4   -2.36       -5   4    6   -11                $ formation
6     0                5:-6:11                         $ outside world



1     cz   1.27
2     cz   2.54
3     cz   3.81
4     c/z  -6.35  0.0   10.16
5     cz   1000.0
6     pz   -500.0
7     pz   15.24
8     pz   22.86
9     pz   38.10
10    pz   63.50
11    pz   800.0


mode   n
imp:n  1  4r  0
sdef   erg=d1
si1    h  0.0   0.15  0.50  1.67  2.24  3.03  3.22  4.09
          5.02  5.91  7.20  7.83  8.72  9.20  10.22
sp1    d  0.0   0.1052183  0.0811733  0.0335948  0.0581514
          0.1139154  0.1427353  0.1021487  0.1193724
          0.0765689  0.0492838  0.0632674  0.0187585
          0.0153479  0.0204639   
m1     2003.50c    1.0
m2     26000.55c   1.0
m3     1001.50c   0.6667   8016.50c   0.3333 
m4     1001.50c   -0.0094     6000.50c   -0.1098  
       8016.50c   -0.5147     20000.51c  -0.3661  
mt3    lwtr.01t
mt4    lwtr.01t
idum    3      1     25
              -1     13
               1     29
               1     1
rdum          -50         60.3
               0.001      51
              -40         90
               0          1e20
wwg    14   3   0.0   0.0   0.0  0.0
wwge:n  1e-7   1e-5   1e-3  2e-1  12
f4:n    1
fm4    -38.6111  1   103  
f14:n   2
fm14   -514.815  1   103 
prdmp  -400       -2000
nps     400000




    The input cards related to the mesh-based weight window are underlined 
in the above input file. The first entry of idum (=3) indicates the option
of WW game and WW map generation with map resetting at each output print. 
The X variable is divided into 25 meshes uniformly from -50.0 to 60.3. 
In Y axis, there is symmetry with respect to Y=0.001 (to avoid overlapping 
with physical plane) plane. To the right of the plane there are 13 meshes 
in the range from 0.001 to 51. There are also 13 meshes on the left of the 
Y=0.001 plane symmetrically. In Z dimension, there are 29 meshes

uniformly distributed  from -40 to 95.1. One time bin from 0 to 1e20 
ignores time-dependency. The importance map is optimized for tally 14, 
i.e. the H-3(n,p) reaction rate. Five energy groups are used for 
importance map generation.



    The test example is calculated by the modified MCNP 4B in an iterative 
fashion. The script file of running the first example is listed below 
(filename = job1). Note that since "imp.in" file does not exist in the 
first iteration, WW game is not played. 




rm -f imp.in
mcnpcear  inp=test1.i outp=test1.o1 runtpe=test1.r
rm -f test1.r
mv imp.out imp.in
mcnpcear  inp=test1.i outp=test1.o2  runtpe=test1.r
rm -f test1.r 
mv imp.out imp.in
mcnpcear  inp=test1.i outp=test1.o3 runtpe=test1.r
rm -f test1.r
mv imp.out imp.in
mcnpcear  inp=test1.i outp=test1.o4 runtpe=test1.r
rm -f test1.r 
mv imp.out imp.in
mcnpcear  inp=test1.i outp=test1.o5 runtpe=test1.r
rm -f test1.r
mv imp.out imp.in
mcnpcear  inp=test1.i outp=test1.o6 runtpe=test1.r
rm -f test1.r
mv imp.out test1.imp  



    The test results on Digital workstations and SUN Ultra 1 workstations 
are presented in Table 1.

Table 1: Test results for the neutron benchmark problem (nps=400,000)
Iteration #	Far detector response	Relative Error	Dec Alpha		SUN Ultra 1		Digital WS 433
Iteration #	Far detector response	Relative Error	CPU time	FOM	CPU time	FOM	CPU time	FOM
1	1.94e-4	0.1015	306.22	0.32	128.21	0.76	29.57	3.28
2	1.95e-4	0.0536	270.32	1.29	113.51	3.07	26.37	13.20
3	2.00e-4	0.0141	232.34	21.0	92.25	54.50	22.02	228.43
4	1.98e-4	0.0109	239.51	35.10	94.50	89.10	22.04	381.89
5	1.99e-4	0.0107	243.89	35.80	96.54	90.5	22.83	382.58
6	1.98e-4	0.0104	208.83	44.30	109.76	84.2	22.07	418.92

It should be pointed out that even though the simulation results (responses, relative errors, etc.) are repeatable the CPU time spent on each iteration may not be the same due to computer performance fluctuation, so the FOM's may not be precisely repeatable.

6. Example 2 -- Gamma-Ray Density Benchmark

The second example is a gamma-ray density benchmark, whose specification can also be found in Reference [4]. The input file for the second example is shown below (filename = test2.i).

Benchmark problem calculation by MCNP-- limestone 20 PU 1 1 -1.0 -2 1 -3 $source chamber 2 2 -3.67 -4 5 -6 $near NaI 3 1 -1.0 1 -7 4 -3 $near collimator 4 2 -3.67 -8 9 -10 $far NaI 5 1 -1.0 1 -11 8 -3 $far collimator 6 3 -19.3 14 -15 -3 #1 #2 #3 #4 #5 $tool 7 4 -1.0 3 -12 14 -15 $borehole 8 5 -2.36 14 -15 12 -13 $formation 9 0 -14:15:13 $outside world

1 px 0.0 2 gq 0.5 1.0 0.5 0.0 0.0 -1.0 0.0 0.0 0.0 -1.0 3 cz 5.0 4 cz 1.0 5 pz 19 6 pz 21 7 c/x 0.0 20.0 1.0 8 cz 2.0 9 pz 38 10 pz 42 11 c/x 0.0 40.0 2.0 12 c/z -5.0 0.0 10.0 13 c/z -5.0 0.0 200. 14 pz -100. 15 pz 200.

mode p imp:p 1 7r 0 sdef erg=0.662 pos=0.0001 0 0 dir=d1 vec=0.707 0 0.707 sb1 -31 3 idum 3 1 30 -1 15 1 30 1 1 rdum -31 60.1 0.001 45.0 -22 70.1 0 1e20 wwg 11 1 0 0 0 40 wwge:p 0.134 0.268 0.402 0.536 0.670 m1 8000.01p 1. m2 11000.01p 0.5 53000.01p 0.5 m3 26000.01p 1. m4 1000.01p 0.66667 8000.01p 0.33333 m5 1000.01p -0.01 6000.01p -0.11 8000.01p -0.514 20000.01p -0.366 f14:p 4 f24:p 4 e24 25e-3 50e-3 75e-3 100e-3 125e-3 150e-3 175e-3 200e-3 225e-3 250e-3 275e-3 300e-3 325e-3 350e-3 375e-3 400e-3 425e-3 450e-3 475e-3 500e-3 525e-3 550e-3 575e-3 600e-3 625e-3 650e-3 675e-3 f1:p 8 c1 0 1 f11:p 8 cut:p j 0.010 prdmp -900 -4000 nps 4000000

In this test problem, X variable is uniformly divided into 30 sections from -30 to 60.1. In Y dimension, there is also a symmetry with respect to Y=0.001 plane, and on either side of the symmetry plane there are 15 3-cm meshes. Z variable is discretized into 30 sections from -22 to 70.1. Again time-dependency is not considered. The WW map is optimized for tally 11, particle current into the cylindrical surface of the far detector. Five energy groups is used to describe the energy dependency of WW map.

The script file of running the second example is listed below (filename =job2).

rm -f runtpe imp.in mcnpcear inp=test2.i outp=test2.o1 runtpe=test2.r rm -f test2.r mv imp.out imp.in mcnpcear inp=test2.i outp=test2.o2 runtpe=test2.r rm -f test2.r mv imp.out imp.in mcnpcear inp=test2.i outp=test2.o3 runtpe=test2.r rm -f test2.r mv imp.out imp.in mcnpcear inp=test2.i outp=test2.o4 runtpe=test2.r rm -f test2.r mv imp.out imp.in mcnpcear inp=test2.i outp=test2.o5 runtpe=test2.r rm -f test2.r mv imp.out imp.in mcnpcear inp=test2.i outp=test2.o6 runtpe=test2.r rm -f test2.r mv imp.out test2.imp

The test results for the second example are given in Table 2.

Table 2: Test results for the gamma-ray benchmark problem (nps=4,000,000)
Iteration #	Far detector response	Relative Error	Dec Alpha		SUN Ultra 1		Digital WS 433
Iteration #	Far detector response	Relative Error	CPU time	FOM	CPU time	FOM	CPU time	FOM
1	8.49e-8	0.7071	144.88	0.014	59.71	0.033	16.55	0.121
2	8.04e-8	0.4669	350.58	0.013	144.88	0.032	37.82	0.121
3	3.26e-8	0.1079	270.57	0.32	106.80	0.80	27.79	3.09
4	3.35e-8	0.0237	443.12	4.0	170.88	10.4	43.39	41.03
5	3.25e-8	0.0094	432.48	26.2	168.64	67.1	42.91	263.75
6	3.26e-8	0.0081	423.90	36.0	166.43	91.6	41.93	363.50

The convergence of the importance map can be sped up by severely biasing the source particle in the direction of the source collimator (0.707 0 0.707), as demonstrated before[2].

7. File List in the Release

1. manual.ps -- PS file of the manual 2. test1.i -- input file for the neutron problem 3. job1 -- script file to run the first test 4. test2.i -- input file for the gamma-ray problem 5. job2 -- script file to run the second test 6. patch.cear -- patch file to create the special MCNP 4B 7. mcnpcear_sun -- SUN-executable of the modified MCNP (not available for the release) 8. mcnpcear_dec -- DEC Alpha-executable of the modified MCNP (not available for the release)

References

[1] Briesmeister, J. F., "MCNP - A General Monte Carlo N-Particle Transport Code, Version 4B," LA-12625-M, Los Alamos National Laboratory (1997).

[2] Liu, L. and R. P. Gardner, "A Geometry-Independent Fine Mesh-Based Monte Carlo Importance Generator," Nucl. Sci. End. 125, 188 (1997).

[3] Liu, L., "Self-Optimizing Monte Carlo Method for Nuclear Well Logging Simulation," Ph.D. Dissertation, North Carolina State University (1997).

[4] Gardner, R. P. and K. Verghese, "Monte Carlo Nuclear Well Logging Benchmark Problems with Preliminary Intercomparison Results," Nucl. Geophys. 4, 4 (1991).

Department of Nuclear Engineering
NC State University
Raleigh, North Carolina