Report

PHITS Multi-Purpose Particle and Heavy Ion Transport code System Advanced Lecture (II): variance reduction techniques to improve efficiency of calculation Jun. 2015 revised title 1 Contents of Lecture １．Introduction ２．Neutron deep penetration calculation [importance] Geometry splitting and Russian Roulette [weight window] Weight windows ３．Calculation of particle production in thin target [forced collision] 2 Neutron deep penetration calculation Calculate neutron transport in thick shield and dose rate distribution to depth Execute “imp.inp” and measure computational time (you can find CPU time in phits.out “total CPU time”) 14 MeV neutron with 1 cm radius maxcas = maxbch = 1500 1 Concrete 50cm radius x 180 cm thick cylinder Air 3 Neutron deep penetration calculation Normal calculation using a single CPU (1.5GHz AMD) imp-dose-xz.eps Number of history= 1,500 total cpu time = 19.53 sec imp-dose-xz.eps Number of history= 270,000 total cpu time = 660.05 sec Need to improve the efficiency of Monte Carlo simulation! Use variance reduction techniques 4 Concept of weight in Monte Carlo calculation Example：track length tally Weight： Importance of the particle in Monte Carlo simulation always to be 1 for normal calculation* Merit of controlling weight Artificially increase the probability of rare event occurrences Kill events that are not so important Demerit of controlling weight Inadequate weight control induces wrong simulation results Frequency distribution per a history cannot be calculated, e.g. [t-deposit] with output = deposit, NO MORE event generator! *Not the case for low-energy neutron transport simulation 5 Contents of Lecture １．Introduction ２．Neutron deep penetration calculation [importance] Geometry splitting and Russian Roulette [weight window] Weight windows ３．Calculation of particle production in thin target [forced collision] 6 Cell Importance Method Set important I to each cell. When a particle passes through the boarder of cell 1 and cell 2, multiple its weight by I1/I2 For I1 < I2, split the particle into I2/I1, and multiple its weight by I1/I2 e.g.1 I2/I1 = 3 (integer) I1=1 W=1 I2=3 W=1/3 W=1/3 W=1/3 • Always split into 3 • Weights of all split particles are 1/3 e.g. 2 I2/I1 = 2.75 (not an integer) • Split into 3 by 75% • Split into 2 by 25% • Weights of all split particles are 1/2.75 For I1 > I2 , play Russian Roulette, and multiple its weight by I1/I2 I1=1 I2=1/3 W=1 W=1 W=1 W=3 e.g.3 I2/I1 =0.33 • 33% of particles survives, rests are killed • Weights of all survived particle are 3 Check the Trajectory of Single Particle Execute “imp-hist.inp”, and check the neutron trajectory 空気 5MeV I =1.0 Neutron 10 コンクリート （R= 20cm， Depth = 8cm×2） I1=1.0 I2=1.0 Importance for all cells is set to 1.0 in the default setting imp-trajectory.eps 8 Let’s Increase Importance imp-hist.inp I10=1 [importance] part = neutron reg imp 10 1 1 2 2 4 I1=2 Reaction occurs here I2=4 divided into 2 neutrons here imp-trajectory.eps 9 Let’s Decrease Importance imp-hist.inp I10=1 [importance] part = neutron reg imp 10 1 1 1/2 2 1/4 I1=1/2 I2=1/4 1/2 neutron survives 1/2 neutron is killed imp-trajectory.eps 10 Let’s Increase Importance EXTREMELY! imp-hist.inp I10=1 Particles are divided too much You waste your machine time without improving statistics! [importance] part = neutron reg imp 10 1 1 2 2 100 I1=2 I2=100 imp-trajectory.eps It is better to set 2~3 for max importance ratio between neighboring cells. “A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380 11 Example of calculation using [importance] Activate the [importance] section in “imp.inp”, and execute PHITS 14 MeV neutron with 1 cm radius Concrete 50cm radius x 180 cm thick cylinder Ii+1/Ii = 2.5 [importance] part = neutron reg imp 1 2.5**0 2 2.5**1 3 2.5**2 4 2.5**3 5 2.5**4 6 2.5**5 7 2.5**6 8 2.5**7 9 2.5**8 10 2.5**9 11 2.5**10 12 2.5**11 Thickness of 1 cell is 15 cm 12 Example of calculation using [importance] 1.50GHz, single [importance] off total history = 1,500 total cpu time = 19.53 sec 1 0.1 Dose [importance] total cpu time = 50.49 sec Relative error 0.01 1 0.1 0.01 More neutrons penetrate to deeper locations ・Check statistical uncertainty13 Red:1.0, Yellow: ~0.1, Green: ~0.01 Example of calculation using [importance] Original data：imp-dose-reg.out Original data：imp-energy-reg.out Dose rate in each cell (15 cm thickness) Neutron energy fluence in cell 3 Agreement in neutron fluence in each cell indicate the adequacy of importance setting 14 Important Notice of Setting [importance] source Good example 1 2 4 8 16 32 8 8 8 32 Bad example 1 1 It is better to set 2~3 for max importance ratio between neighboring cells. Reference: 15 “A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380 Bad Example for using [importance] Bad example: Very large importance gap between cell 5 & 6 Let’s try a bad example! Dose Relative error [[importance] part = neutron reg imp 1 2.5**0 2 2.5**0 3 2.5**0 4 2.5**0 5 2.5**0 6 2.5**5 7 2.5**6 8 2.5**7 9 2.5**8 10 2.5**9 11 2.5**10 12 2.5**11 Relative errors are too large in comparison to previous setting! Contents of Lecture １．Introduction ２．Neutron deep penetration calculation [importance] Geometry splitting and Russian Roulette [weight window] Weight windows ３．Calculation of particle production in thin target [forced collision] 17 Difference between cell importance and weight window methods • • “Cell importance method” assign a single value of weight for each cell. “Weight window method” assign allowed weight range (window) for each cell and each energy group. Efficient simulation with focusing on the important energy region, such as high-energy neutron. W2=1/3 region1 region2 Number of particle 1 → I2/I1=3 Cell importance method W1=1 Weight W1=1 WU=2.5 WU=1.25 W=1 Weight Weight I1=1 WU=2.5 I2=3 WL=0.5 WL=0.5 WL=0.25 region1 energy group1 W2=0.5 region2 energy group2 Number of particle 1 → 1 WU=0.75 WL=0.15 1→ W1/W2=2 Weight window method 18 Parameters Related to Weight Window wupn: Ratio of the upper and lower limits of allowed weight (D=5.0) mxspln: Maximum number of split per event (D=5.0) mwhere: Timing for considering weight window (D=0) -1: reaction, 0: both, 1: crossing surface See manual “4.2.4 Cut off time, cut off weight, and weight window” in more detail 19 Check the Trajectory of Single Particle Let’s execute “weight-hist.inp”, and check the neutron trajectory 5MeV neutron Air Water （r= 10cm， Depth = 5cm×2） weight-trajectory.eps 20 Set [Weight Window] Section Activate the 1st [weight window] section in weight-hist.inp, and execute PHITS! Lower limit of allowed weights, WL, for each cell is defined in [weight window] Upper limit of allowed weights, WU, is automatically set to WL x wupn (D=5) weight-hist.inp [weight window] part = neutron reg ww1 1 0.25 2 0.25 Same weight window for all energies WU=5 x WL =1.25 Weight W1=1 WL=0.25 W =0.25 L Cell 1 Cell 2 weight-trajectory.eps 21 Weight of neutron is within the weight window range -> No split and No Russian Roulette Let’s Decrease Weight Window in Cell 2 Decrease the lower limit of allowed weights in cell 2, and execute PHITS! weight-hist.inp [weight window] part = neutron reg ww1 1 0.25 2 0.133 WU=1.25 Weight W1=1 WU=0.67 WL=0.25 WL=0.133 Cell 1 Cell 2 Particle Split weight-trajectory.eps Neutron is divided into 2 because its weight (=1) is higher than the upper limit (=0.67)22 Let’s Decrease Weight Window in Cell 2 EXTREMELY! Decrease the lower limit of allowed weights in cell 2 EXTREMELY, and execute PHITS! weight-hist.inp [weight window] part = neutron reg ww1 1 0.25 2 0.005 WU=1.25 Weight W1=1 WL=0.25 WU=0.025 WL=0.005 Cell 1 Cell 2 weight-trajectory.eps Divided into 40 neutrons (but maximum number of split per event = 5) It is better to set 2~3 for max weight window ratio between neighboring cells 23 Set different Weight Window for Each Energy Group weight-hist.inp Let’s “off” the 1st [weight window], and activate the 2nd one & execute PHITS 2 energy groups are defined But the weight window are the same for each group WU=1.25 WU=1.25 W1=1 W1=1 WL=0.25 WL=0.133 Cell 1 Cell 2 Particle Split E < 0.01 MeV WU=0.67 Weight WU=0.67 Weight [weight window] Maximum energy part = neutron for each group eng = 2 0.01 1.0e5 reg ww1 ww2 1 0.25 0.25 2 0.133 0.133 WL=0.25 WL=0.133 Cell 1 Cell 2 Particle Split 0.01 MeV < E < 105MeV weight-trajectory.eps Neutron is divided into 2 because its weight (=1) is higher than the upper limit (=0.67) (Same as 2 pages before) 24 Increase Weight Window Only for Low-Energy Neutrons weight-hist.inp Increase weight window for low-energy neutrons by 10 times, and execute PHITS WU=12.5 WU=6.65 WL=2.5 W1=1 WU=1.25 WL=1.33 W1=1 Weight Weight WU=0.67 Cell 1 E < 0.01 MeV Low-energy neutron is killed WL=0.25 WL=0.133 Cell 2 Russian Roulette [weight window] part = neutron eng = 2 0.01 1.0e5 reg ww1 ww2 1 0.25*10 0.25 2 0.133*10 0.133 Cell 1 Cell 2 Particle Split 0.01 MeV < E < 105MeV weight-trajectory.eps Reduce the production of low-energy neutrons, and concentrate on high-energy neutron transport simulation 25 Get poor and better statistics at shallower and deeper locations, respectively Example of calculation using [weight window] weight.inp Execute PHITS using weight.inp, and check your CPU time [weight window] set: c10[1.0] part = neutron eng = 2 1.0e-3 reg ww1 1 c10*2.5**(-1) 2 c10*2.5**(-2) 3 c10*2.5**(-3) ・・・ Geometry, source, history are the same as imp.inp Cylindrical Concrete (r = 50cm, depth = 1.8 m) 1e5 ww2 2.5**(-1) 2.5**(-2) 2.5**(-3) 14MeV neutron Depth of each cell = 15 cm total cpu time = 25.89 sec Decrease weight window of neighboring cells by 1/2.5 for particle split Set same weight window for all energy groups (because c10=1.0) weight-dose-xz.eps weight-dose-xz_err.eps 26 Example of calculation using [weight window] weight.inp Increase the weight window for low-energy neutrons by 10 times, and execute PHITS [weight window] set: c10[10.0] part = neutron eng = 2 1.0e-3 reg ww1 1 c10*2.5**(-1) 2 c10*2.5**(-2) 3 c10*2.5**(-3) ・・・ Most low-energy neutrons are immediately killed when they produced, because of Russian Roulette 1e5 ww2 2.5**(-1) 2.5**(-2) 2.5**(-3) weight-dose-xz.eps Low-energy neutrons have short range, so statistical uncertainties of deeper locations do not change that much total cpu time = 25.89 -> 14.53 sec weight-dose-xz_err.eps 27 Comparison of CPU time Normal calculation (No [importance], No [weight window]) total cpu time = 19.53 sec Use [importance] total cpu time = 50.49 sec Use [weight window]; same value for all energies total cpu time = 25.89 sec Use [weight window]; different values for low- and high-energy neutrons total cpu time = 14.53 sec 28 Dose Reduction Rate in Concrete 90 cm 180 cm Good agreement except for deeper location doses obtained from default setting You can improve the efficiency of your calculation by appropriately setting weight window for each neutron energy group 29 Contents of Lecture １．Introduction ２．Neutron deep penetration calculation [importance] Geometry splitting and Russian Roulette [weight window] Weight windows ３．Calculation of particle production in thin target [forced collision] 30 Forced collision The forced collision is useful for analyzing secondary particles generated from a thin target Incoming weight Wi Split into two particles d: distance across cell Uncollided particle Collided particle Forced collision cell Weight of uncollided particle： Wi×exp(-σd) Weight of collided particle： Wi×{1-exp(-σd)} σ: macroscopic cross section Collide position is decided by cross section and random number. Forced collision factor |fcl| 1 fcl = 0: no forced collision in cell, |fcl| = 1: 100% forced collision •fcl < 0 = applies only to particles entering the cell (weight cut-off is not applied) •fcl > 0 = applies to all particles surviving weight cutoff (weight cut-off is applied) 31 Example of forced collision Energy distribution of neutron and alpha produced by reaction of 100MeV proton incidence on 1mm thick Si. force.inp maxcas = 50000 maxbch = 2 Secondary particle flux calculated without forced collision： No collision occurred in such thin target [forced collisions] part = proton reg fcl 1 1.0 Secondary particle flux calculated with forced collision： You can get good statistic data!! 32 Summary Cell importance and weight window methods are effective in deep penetration calculations. Ratios of importance and weight window between neighboring cells are better to be be less than 3. Weight window method with energy dependence is more efficient than other methods for deep penetration calculations. Forced collision method is effective to calculate particle production in a thin target. 33