search

nrc-cnrc / EGSnrc

Toolkit for Monte Carlo simulation of ionizing radiation — Trousse d'outils logiciels pour la simulation Monte Carlo du rayonnement ionisant
http://nrc-cnrc.github.io/EGSnrc
GNU Affero General Public License v3.0
216 stars 146 forks source link

dosxyznrc n_split segfault #566

Closed jw3126 closed 1 week ago

jw3126 commented 4 years ago

DOSXYZnrc seems to segfault if n_split > 1 and the phase space plane is outside of the phantom.

some_title                                                                       #!GUI1.0
1
H2O521ICRU
0.521, 0.01, 0, 0
-1, -1, -1, 0
-1
1, 2
-1
1, 2
-1
1, 2
0, 0, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 0, 0
0, 0, 0, 0, 0, 0, 0, 0
0, 2, 0, 0, 0, 180, 0, 1.1, 180, 0, 0, 0, 0, 0
2, 0, 0, 1000, 0, , ,
$EGS_HOME/dosxyznrc/tiny.egsphsp
100000, 0, 999, 0, 1, 100, 0, 0, 0, 0, , 0, 1, 0, 10, 1, 0
 #########################
:start MC Transport Parameter:
    Global ECUT = 0.521
    Global PCUT = 0.01
    Global SMAX = 5
    ESTEPE = 0.25
    XIMAX = 0.5
    Boundary crossing algorithm = PRESTA-I
    Skin depth for BCA = 0
    Electron-step algorithm = PRESTA-II
    Spin effects = On
    Brems angular sampling = Simple
    Brems cross sections = BH
    Bound Compton scattering = Off
    Compton cross sections = default
    Pair angular sampling = Simple
    Pair cross sections = BH
    Photoelectron angular sampling = Off
    Rayleigh scattering = Off
    Atomic relaxations = Off
    Electron impact ionization = Off
    Photon cross sections = xcom
    Photon cross-sections output = Off
:stop MC Transport Parameter:
 #########################
$ dosxyznrc -i segfault -p 521icru
================================================================================
EGSnrc version 4 for x86_64-unknown-linux-gnu          Thu Jan  2 07:08:36 2020
================================================================================
configuration...............................................ansible
user code...................................................dosxyznrc
pegs file...................................................521icru on HEN_HOUSE
using host..................................................ptwlinuxvm
input file..................................................segfault
output file(s)..............................................segfault
================================================================================

 Begin execution with large arrays being zeroed
 Thisis only needed for Linux g77 compiler - comment
 this code near the top of dosxyznrc.mortran if you are
 not using a linux g77 compiler

 *******************************************************************************
 NRCC/UW EGSnrc user-code DOSXYZnrc
 ON linux (ansible)                                        07:08:36 Jan  2 2020
 *******************************************************************************
 **                                                                           **
**                                   DOSXYZnrc                                **
 **                              Z pronounced zed                             **
 **                                                                           **
 **      Code developed at the National Research Council of Canada and        **
 **           University ofWisconsin as part of the OMEGA project             **
 **                                                                           **
 **                                                                           **
 **                                                                           **
 *******************************************************************************

     The following parameters may be adjusted in dosxyz_user_macros.mortran
 $MXMED:    Max number of media:  7
 $MXSTACK:  Max stack size:     10000
 $IMAX,etc: Max dose scoring regions in x,y,z directions:  360  360  360
 $MAXDOSE:  Max dose scoring regions consistent with above:*******
 $DOSEZERO(=1) 1=> all doses with uncert > 50% are zeroed in .3ddose file

 The following parameters may be adjusted in srcxyz.macros
 $INVDIM:   number of elements in inverse CPD for input energy spectra = 1000
 $NENSRC:   number of bins in input energy spectrum =  200

 ===============================================================================
 -------------------------------------------------------------------------------

 Title:  some_title                                                                      
 -------------------------------------------------------------------------------
 ===============================================================================

 Number of media (min = 1, max =   7, 0 => CT data):                1
 Medium  1:                     H2O521ICRU              

 ECUTIN,PCUTIN,(ESTEPE,SMAX--DUMMY INPUTS): 
              0.521     0.010     0.000     0.000

 # regions in x (max= 360),y (max= 360),z (max= 360) directions
 (if<0,implies # groups of reg), IPHANT (1 to output a .egsphant
 file for dosxyz_show, 0[default] to not output this file)
     :      -1    -1    -1     0

 Input boundaries in the x-direction
 -----------------------------------
 Initial boundary:        -1.000
 Width in this group, number of regions in group:         1.000    2
 Boundaries
      -1.000       0.000       1.000

 Input boundaries in the y-direction
 -----------------------------------
 Initial boundary:        -1.000
 Width in this group, number of regions in group:         1.000    2
 Boundaries
      -1.000       0.000       1.000

 Input boundaries in the z-direction
 -----------------------------------
 Initial boundary:        -1.000
 Width in this group, number of regions in group:         1.000    2
 Boundaries
      -1.000       0.000       1.000

 Total # regions including exterior =         9

 Input groups of regions for which density and medium are not defaults
 Lower,upper i, j, k,  MEDIUM, DENSITY
    Found blank line => end of this input

 Input groups of regions for which ECUT and PCUT are not defaults
 NB This option is disabled, just input 8 zeros.
 Dummy values of lower,upper i, j, k,  ECUT, PCUT
    Found blank line => end of this input

 Enter 8 numbers on one line
   3 pairs defininglower,upper x,y,z indicies of dose regions
                 for which results are to be output
   IZSCAN:       non-zero for z-scan/page
   MAX20:        if any one = 1, output summary of max 20 doses.
   end signaled by first pair both zero
   forno dose printed, MAX20 is still read from first line

    Found blank line => end of this input
 No doses will be output to egslst file

 Source configuration

      (0) Parallel, rectangular beam incident from the front
          Requires 9 inputs:
          charge (-1,0,1),
          0 (mandatory, to identify source type),
   lower x-coordinate of the beam (cm),
          upper x-coordinate of the beam (cm),
          lower y-coordinate of the beam (cm),
          upper y-coordinate of the beam (cm),
         angle of beam with respect to the positive x-axis (degrees),
          angle of beam with respect to the positive y-axis (degrees),
          angle of beam with respect to the negative z-axis (degrees) 
            (angles default to 90,90,0--incident on front of phantom)

   or (1) Parallel, rectangular beam incident from any direction
          Requires 10 inputs:
          charge (-1,0,1),
          1 (mandatory, to identify source type),
       x-coordinate of the isocenter (cm),
          y-coordinate of the isocenter (cm),
          z-coordinate of the isocenter (cm),
          angle between +z direction and the line joining the   
             center of the beam (collimator) to the isocenter
             --called the polar angle(degrees),
          angle between +x direction and the projection of the
             line joining the center of the beam (collimator)
             to the isocenter on the xy plane--called the azimuthal
             angle (degrees),
          total x-widthof the beam in the plane perpendicular
             to the beam direction (cm),
          total y-width of the beam in the plane perpendicular
             to the beam direction (cm),
          angle by which the collimator is rotated in the
           collimator plane perpendicular to the beam
     direction (degrees),                     
             (+ve rotation is counterclockwise looking along
              the beam direction

  or  (2) Full phase-space of each particle
Requires 9 inputs plus data stored in units 43 and 44:
  charge (-1 electron,0 photon,1 positron, 2 all),
          2 (mandatory, to identify source type),
          x-coordinateof the isocenter (cm),
          y-coordinate of the isocenter (cm),
          z-coordinate of the isocenter (cm),
     angle between +z direction and the line joining the   
          origin in the phase space plane to the isocenter
          --called the polar angle(degrees),
          angle between +x direction and the projection of the
             line joining the origin in the phase space plane
             tothe isocenter on the xy plane--called the azimuthal
  angle (degrees),
          absolute distance from the isocenter to the origin
             in the phase space plane               
          angle by which the source is rotated in the
             phase space plane perpendicular to the beam
             direction (degrees),                     
             (+ve rotation is counterclockwise looking down
            from the origin in the phase space plane),
   i_dbs--set to 1 if DBS was used in BEAM simulation used
          to generate the phsp source and you want to reject fat
             photons, 0 otherwise,
          DBS splitting radius (cm),
          SSD at which splitting radius defined (cm),
          Z at which phsp source collected (cm),
  No. of times to split charged particles.

   or (3) Point, rectangular beam incident from the front
          Requires 7 inputs:
          charge (-1,0,1),
          3 (mandatory, to identify source type),
lower x-coordinate of the beam (cm),
          upper x-coordinate of the beam (cm),
          lower y-coordinate of the beam (cm),
          upper y-coordinate of the beam (cm),
      distance to the plane (cm),

  or  (6) Uniform isotropically radiating parallelepiped within
          the phantom
          Requires 8 inputs:
          charge (-1,0,1),
          6 (mandatory, to identify source type),
          lower x-coordinate of active volume (cm)
          upper x-coordinate of active volume (cm),
       lower y-coordinate of active volume (cm) 
          upper y-coordinate of active volume (cm),
          lower z-coordinate of active volume (cm)  
          upper z-coordinate of active volume (cm) 

  or  (7) Parallel beam incident from multiple, user-selected angles
          Requires 9 inputs on this line:
   charge (-1,0,1),
          7 (mandatory, to identify source type),
          x-coordinate of the isocenter (cm),
      y-coordinate of the isocenter (cm),
          z-coordinate of the isocenter (cm),
          number of incident theta-phi pairs or -ve number of
             groups of incident theta-phi pairs where, within a group
             only theta or phi can vary, the varying angles are
             evenly distributed and have equal probability,
          total x-width of the beam in the plane perpendicular
             to the beamdirection (cm),
          total y-width of the beam in the plane perpendicular
             to the beam direction (cm),
          angle by which the collimator is rotated in the
         collimator plane perpendicular to the beam
   direction (degrees),                     
             (+ve rotation is counterclockwise looking along
              thebeam direction

  or  (8) Full phase-space incident from multiple angles
          Requires 8 inputs on this line and data stored in units 43,44:
          charge (-1 electron,0 photon,1 positron, 2all),
          2 (mandatory, to identify source type),
        x-coordinate of the isocenter (cm),
          y-coordinate of the isocenter (cm),
          z-coordinate of the isocenter (cm),
          number of incident theta-phi pairs or -ve number of
             groups of incident theta-phi pairs where, within a group
             only theta or phi can vary, the varying angles are
             evenly distributed and have equal probability,
          absolute distance from the isocenter to the origin
             in the phase space plane                  
          angle by which the source is rotated in the
             phase space plane perpendicular to the beam
             direction (degrees),                     
             (+ve rotation is counterclockwise looking down
              from the origin in the phase space plane),
      i_dbs--set to 1 if DBS was used in BEAM simulation used
             to generate the phsp source and you want to reject fat
             photons, 0 otherwise,
          DBS splitting radius (cm),
          SSD at which splitting radius defined(cm),
          Z at which phsp source collected (cm)
     No. of times to split charged particles.

  or  (9) BEAM simulation of treatment head
Requires 11 inputs plus name of accelerator simulation,
   input file used in accelerator simulation, and pegs4
   data used in accelerator simulation:
          charge (-1 electron,0 photon,1 positron, 2 all),
          9 (mandatory, to identify source type),
          x-coordinate of the isocenter (cm),
          y-coordinate of the isocenter (cm),
      z-coordinate of the isocenter (cm),
          angle between beam central axis and +z axis in DOSXYZ
             geometry--called the polar angle(degrees),
          angle between +x direction in DOSXYZ geometry and
             beam central axis projected on the DOSXYZ xy plane
             --called the azimuthal angle (degrees),
          absolute distance from the isocenter to centre of
             scoring plane in BEAM simulation,
          angle to rotate BEAM simulation about its central
             axis (degrees) (+ve rotation is counterclockwise
             looking down the axis),
          i_dbs--set to 1 if DBS is being used in BEAM simulation
    and you want to reject fat photons, 0 otherwise,
No. of times to split charged particles.

  or  (10) BEAM simulation source incident from multiple angles
          Requires 10 inputs plus name of accelerator simulation,
          input file used in accelerator simulation,and pegs4
          data used in accelerator simulation:
       charge (-1 electron,0 photon,1 positron, 2 all),
    9 (mandatory, to identify source type),
          x-coordinate of the isocenter (cm),
          y-coordinate of the isocenter (cm),
          z-coordinate of the isocenter (cm),
          number of incident theta-phi pairs or -ve number of
             groups of incident theta-phi pairs where, within a group
             only theta or phi can vary, the varying angles are
             evenly distributed and have equal probability,
          absolute distance from the isocenter to centre of
             scoring plane in BEAM simulation,
   angle to rotate BEAM simulation about its central
    axis (degrees) (+ve rotation is counterclockwise
   looking down the axis),
          i_dbs--set to 1 if DBS is being used in BEAM simulation
             and you want to reject fat photons, 0 otherwise,
          No. of times to splitcharged particles.

  or  (20) Phase Space Incident from multiple settings 
  optionally through an MLC or through a BEAM accel.
 Requires 6 inputs plus name of the input file used
       forBEAM/vcu SIM. and the BEAM/VCU code if used
          data used in simulation:
          charge (-1 electron,0 photon,1 positron, 2 all),
          20 (mandatory, to identify source type),
          number of control points, 
          i_dbs: set to 1 if DBS is being used in simulation
             and you want to reject fat photons, 0 otherwise,
          r_dbs: radius of DBS splitting field in original
             BEAM simulation,
          ssd_dbs: SSD of DBS splitting field, 
   z_dbs: Z position where phase space was scored 
 in original BEAM simulation, 
          No. of times to splitcharged particles,
          i_muidx_out: Set to 1 to include fractional MU index
             in output phase space (i_phsp_out=1 or 2)
          calflag: Set to 1 to skip the calibrationrun performed
             to refine the estimate of NRCYCL.

  or  (21) BEAM simulation of treatment head will multiple settings
  optionally through a MLC
          Requires  name of accelerator simulation,
          input file used in accelerator simulation, and pegs4
          data used in accelerator simulation, name of input file
       for vcu SIM. and the VCU code (particleDmlc) if using MLC
          charge (-1 electron,0 photon,1 positron, 2 all),
          21 (mandatory, to identify source type),
          number of control points 
         i_dbs: set to 1 if DBS is being used in BEAM simulation
             and you want to reject fat photons, 0 otherwise,
          No. of times to split charged particles.
i_muidx_out: Set to 1 to include fractional MU index
   in output phase space (i_phsp_out=1 or 2)

 All inputs on one line:
 Full phase space from an arbitrary direction
 x-coordinate of the isocenter:                           0.0000
 y-coordinate of the isocenter:                           0.0000
 z-coordinate of the isocenter:                           0.0000
 Polar angle of origin in source plane:                 180.0000
 Azimuthal angle of origin in source plane:               0.0000
 Distance from isocenter to origin in source plane:       1.1000
 Source rotation angle:                                 180.0000
 No. of times to split e+/e-:                                  0

 Charge of incident particles to use:  0

 Enflag(0=mono-E,1=spectr,2=phsp or full BEAM sim.,3=dosecomp,4=beam model),
 Mode(0,2), medsur(0 = vacuum), dsurround(1), dflag(0=1 dsurround,
 1=4 dsurrounds), dsurround(2), dsurround(3), dsurround(4)
 :           2           0           01000.000           0

 The material in the region outside the phantom is vacuum.
 The thickness of this region (in x, y & z direction) is:1000.000 cm

 Input name of phase space file (A256)
  : 
 Particles will be read from file: 
$EGS_HOME/dosxyznrc/tiny.egsphsp                                                                                                                                                                                                                                

 Retrieving file: /home/ptwlinuxvm/EGSnrc/egs_home//dosxyznrc/tiny.egsphsp

 ******WARNING*****
 PHASE SPACE SOURCE WAS GENERATED USING AN OLDER
 VERSION OF BEAM IN WHICH # OFPARTICLES INCIDENT FROM
 ORIGINAL SOURCE WAS NOT STORED.  THISNUMBER WILL
 BE SET EQUAL TO THE NUMBER OF PARTICLES IN THE PHASE
 SPACE SOURCE, AND DOSE WILL BE NORMALIZED
 WITH RESPECTTO THE # OF HISTORIES RUN...NOT
 THE # OF INCIDENT PARTICLES FROM THE ORIGINAL SOURCE

 Total number of particles in file      :            1
 Total number of photons                :            1
The rest are electrons/positrons.

 Maximum kinetic energyof the particles:             1.000 MeV
 Minimum kinetic energy ofthe electrons:          Infinity MeV
 # of particles incident fromoriginal source:          1.0

 NCASE,IWATCH,TIMMAX,INSEED1,INSEED2,BEAM_SIZE,ISMOOTH,IRESTART,IDAT,
 IREJECT,ESAVE_GLOBAL,NRCYCL,IPARALLEL,PARNUM,n_split,ihowfarless,i_phsp_out
 :  First RN seed outside allowed range and default value set

      100000   0 999.00      1802         1 100.00      0    0   0   0   0.00****   1   0  10   1   0

***************** Warning: 
File /home/ptwlinuxvm/EGSnrc/egs_home/dosxyznrc/egsrun_15748_segfault_ptwlinuxvm/segfault.errors
 is already opened and connected to unit  15
 Will not try to re-open this file, assuming it has been opened
 by specifying it in the .io file.

 Bound Compton start region
  Setting all to            0
 Rayleigh start region
  Setting all to            0
 Relaxations start region
  Setting all to            0
 PE sampling start region
  Setting all to            0

 Call hatch
 ----------

 ===> Photonuclear flag:   0
 Rayleigh data available for medium  1 in PEGS4 data set.

(Re)-initializing photon cross sections with files from the series: xcom

 Compton cross sections: default
 Using Compton cross sections from /home/ptwlinuxvm/EGSnrc/HEN_HOUSE/data/compton_sigma.data
 Working on medium   1 ... OK
  old PRESTA calculates default min. step-size for BCA: 
      minimum ECUT found:   0.52100000000000002     
      default BLCMIN is:    3.2034336048218699     
      this corresponds to    24.616909991620847       elastic MFPs 

Reading screened Rutherford MS data ...............  done 

Reading spin data base from /home/ptwlinuxvm/EGSnrc/HEN_HOUSE/data/spinms.data
EGSnrc spin data, version 2.0   
Data generated on a machine with 1234 endianess
The endianess of this CPU is 1234
Ranges:      1.00   100.00  0.30054  1.00000

  medium    1 .....................  done

  Medium            1  sige =    6.3784170227787529        6.1807995067678911       monotone =  F F

  Initializing tmxs for estepe =   0.25000000000000000       and ximax =   0.50000000000000000     

Output from subroutine EDGSET:
==============================
 Atomic relaxations not requested! 

 Bound Compton scattering not requested! 

 EGSnrc SUCCESSFULLY 'HATCHED' FOR ONE MEDIUM.

================================================================================

                   Electron/Photon transport parameter

================================================================================

 Photon cross sections                                      xcom
 Compton cross sections                                     default
 Photon transport cutoff(MeV)                                    0.1000E-01
 Pair angular sampling                                       SIM
 Pair cross sections                                         BH 
 Triplet production                                          Off
 Bound Compton scattering                                    OFF           
 Radiative Compton corrections                               Off           
 Rayleigh scattering                                         OFF           
 Atomic relaxations                                          OFF           
 Photoelectron angular sampling                              OFF           

 Electron transport cutoff(MeV)                               0.5210
 Bremsstrahlung cross sections                              BH  
 Bremsstrahlung angular sampling                             SIM
 Spin effects                                                On
 Electron Impact Ionization                                  Off
 Maxium electron step in cm (SMAX)                                5.000    
 Maximum fractional energy loss/step (ESTEPE)                0.2500
 Maximum 1st elastic moment/step (XIMAX)                     0.5000
 Boundary crossing algorithm                                 PRESTA-I  
 Skin-depth for boundary crossing (MFP)                      24.62    
 Electron-step algorithm                                     PRESTA-II 

================================================================================

 ****WARNING****
 ECUTIN > ECUT input in EGSnrc parameters (     0.5210 MeV).
 ECUT defaults to ECUTIN.

 Will recyle each phase space particle **** times before going on 
 to next particle.
 Restarted phase space files and recycled particles are not
 redistributed

 Starting a new calculation

 Store intermediate files for each batch

 *******************************************************************************

         Summary of source parameters (srcxyznrc)
 *******************************************************************************

               Full phase space input for each incident particle

          x-coordinate of the isocenter,                         0.0000 cm
          y-coordinate of the isocenter,                         0.0000 cm
          z-coordinate of the isocenter,                         0.0000 cm
          Polar angle of source plane:                         180.0000 degrees
          Azimuthal angle of origin in source plane:             0.0000 degrees
          Distance from isocenter to origin in source plane:    -1.1000 cm
          Source plane rotation angle,                         180.0000 degrees
          Total number of particles in phase space file:                1

 Particles to be simulated: photon only

 Medium                AE        AP
 H2O521ICRU          0.521     0.010

 No range rejection.

 Photons will be split   10 times

 ***************************************************************

  Histories to be simulated for this run       100000

  Histories to be analyzed after this run      100000

 ***************************************************************
   Elapsed wall clock time to this point=       0.400 s

   CPU time so far for this run =       0.400 s

 BATCH #  TIME-ELAPSED  TOTAL CPUTIME  RATIO  TIME OF DAY  RNG pointers

     1          0.0            0.0      0.00    07:08:36   ixx jxx =   97  33 

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:
#0  0x7f323b1b231a
#1  0x7f323b1b1503
#2  0x7f323a82ef1f
#3  0x55a5e6d1222d
#4  0x55a5e6d12f3c
#5  0x55a5e6d79d86
#6  0x55a5e6cd980e
#7  0x7f323a811b96
#8  0x55a5e6cd9839
#9  0xffffffffffffffff
Segmentation fault (core dumped)
ojalaj commented 4 years ago

You seem to have a warning of using an old phsp file format and I'm not sure is this a problem or not: "Minimum kinetic energy ofthe electrons: Infinity MeV".

Could you try with any "normal" phsp file?

jw3126 commented 4 years ago

@ojalaj good point. I tested and it also segafaults with a "normal" phsp file.

ojalaj commented 4 years ago

@jw3126 OK...does simulation run otherwise normally with nsplit=1? I see that you have increased IMAX,JMAX,KMAX from default values to 360. Your dsurround is quite large as well - 1000cm. Could these have some effect. Could you run your simulation starting with some of the example files and one by one changing the values towards your input file above.

jw3126 commented 4 years ago

Yeah it also segfaults if I don't touch IMAX etc. Starting with an example file is a good idea.

blakewalters commented 1 week ago

@rtownson, @jw3126: I'm unable to reproduce this bug using the current version of DOSXYZnrc, even in the diabolical case where my phase space file consists of 1 photon. @jw3126, can you still reproduce this issue? If yes, please send me your .egsinp and .egsphsp1 files. Otherwise, let's close this issue.