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
233 stars 146 forks source link

egs++ cd geometry inscribed in envelope geometry occasionally fails #367

Open ojalaj opened 7 years ago

ojalaj commented 7 years ago

Referencing to https://plus.google.com/108037604399060498077/posts/ZY3qZ9imN2B from my PhD student, the example below works with parallel source, but not with IAEA phsp sources. Do we have a bug here or a user error from our side?

###############################################################################
# Testing CD Geometry
# No VRT in use
###############################################################################

###################################
# Define the simulation geometries.
################################### 

:start geometry definition:

    # The base geometry, this will be the Chopping Device (CD)
    # The base geometry can be any geometry, even a composite one
    :start geometry:
        name        = my_cd_planes
        library     = egs_planes
        type        = EGS_Zplanes
        positions   = 3 9 11
        # No media required
    :stop geometry:

    :start geometry:
        name        = my_cd_cylinder
        library     = egs_cylinders
        type        = EGS_ZCylinders
        radii       = 1.6 2
        :start media input:
            media = AIR521ICRU H2O521ICRU
            set medium = 1 1
        :stop media input:
    :stop geometry:

    :start geometry:
        name        = my_cd_sphere
        library     = egs_spheres
        midpoint = 0 0 9
        radii = 1.6 2
        :start media input:
            media = AIR521ICRU H2O521ICRU
            set medium = 1 1
        :stop media input:
    :stop geometry:

    # The composite geometry
    :start geometry:
        name            = chamber1
        library         = egs_cdgeometry
        base geometry   = my_cd_planes
        # set geometry = 1 geom means:
        # "in region 1 of the basegeometry, use geometry "geom"
        set geometry   = 0 my_cd_cylinder
        set geometry   = 1 my_cd_sphere
        # The final region numbers are attributed by the cd geometry object;
        # Use the viewer to determine region numbers
    :stop geometry:

#----------------------------------------------------------------

##### Surrounded by a huge air region #####
    :start geometry:
       name = world_of_air
       library  = egs_box
       box size = 200
       :start media input:
              media = AIR521ICRU
       :stop media input:
    :stop geometry:

##### chamber geometry inscribed into water phantom #####

    :start geometry:
    library = egs_genvelope
    name = dose_to_chamber1
    base geometry = world_of_air
    inscribed geometries = chamber1
     :stop geometry:

simulation geometry = dose_to_chamber1 

:stop geometry definition:

###########################################
#  The source
###########################################

:start source definition:

#    :start source:
#        library = egs_parallel_beam
#        name    = the_source
#        charge  = 0
#        :start shape:
#            library = egs_rectangle
#            rectangle = -5 -5  5  5
#            :start transformation:
#                translation = 0 0 -20
#            :stop transformation:
#        :stop shape:
#        :start spectrum:
#            type = monoenergetic
#            energy = 2
#        :stop spectrum:
#    :stop source:

# simulation source = the_source

    :start source:
        library = iaea_phsp_source
        name    = BEAM_iX_6MV
        iaea phase space file = /../VarianClinaciX_6MV_10x10_w1
        particle type = all
    :stop source:

   simulation source = BEAM_iX_6MV

:stop source definition:

############################################ 
# Run control
############################################
:start run control:

    ncase = 1000
:stop run control:

#############################################
# Scoring
############################################
:start scoring options:

    calculation type = dose

    # smaller cavity
    :start calculation geometry:
    geometry name = dose_to_chamber1
    cavity regions = 1
    cavity mass = 0.002987025891
    :stop calculation geometry: 

:stop scoring options:

############################################
# Variance reduction
############################################

    :start variance reduction:
    cs enhancement = 0
    TmpPhsp = 0
    :stop variance reduction:

############################################
# Transport Parameters
############################################
:start MC transport parameter:

 Global ECUT= 0.521
 Global PCUT= 0.001
 Global SMAX= 1e10
 ESTEPE= 0.25
 XImax= 0.5
 Skin depth for BCA= 3
 Boundary crossing algorithm= EXACT
 Electron-step algorithm= PRESTA-II
 Spin effects= on
 Brems angular sampling= KM
 Brems cross sections= NIST
 Photon cross sections= xcom
 Electron Impact Ionization= On
 Triplet production= On
 Radiative Compton corrections= On
 Bound Compton scattering= On
 Pair angular sampling= KM
 Pair cross sections= NRC
 Photoelectron angular sampling= On
 Rayleigh scattering= On
 Atomic relaxations= On
 Photonuclear attenuation= On

:stop MC transport parameter:
 ########################
mchamberland commented 7 years ago

I can reproduce this in tutor7pp, even by moving the z=11 plane to 12 (to avoid overlap with the tip of the sphere).

mchamberland commented 7 years ago

This can also be reproduced when electron transport is turned off. As far as I can tell, this will require actual debugging and following the particle steps through the geometry. : /

ojalaj commented 7 years ago

@mchamberland, thanks for debugging! Do you have non-IAEA phsp files to test if this happens also with them? I apologize @blakewalters that once again we are producing IAEA phsp source related issues :)

mchamberland commented 7 years ago

Unfortunately, I don't since I only ever use IAEA phsp files.

mchamberland commented 7 years ago

I think the IAEA phase space source is a red herring. I've followed a photon through the CD geometry and this seems like just another edge case that the CD geometry howfar logic fails to properly account for. Thankfully, the geometry is very simple to follow so we should be able to find a solution for this...

(But I'm actually on vacation, so I won't have much time to spend on this in the next 2 weeks...)

mchamberland commented 7 years ago

As far as I can tell, the howfar of CD geometry is okay and working as intended. It's the floating point errors that end up being relatively large compared to other situations I've encountered before: ~1e-5. So, using a larger boundary tolerance than usual (say, 1e-4) seems to fix the issue (at least, up to 1e8 histories).

blakewalters commented 7 years ago

I tested this geometry (actually, a simplified version of it) with both the IAEA phase space source you're using and an EGS format version of it, and both produce the CD geometry error (albeit fewer errors in the case of the EGS format phase space). And while this may be an ideal situation in which to (further) firm up the roundoff error checking in EGS_CDGeometry, in the meantime, I would follow the suggestion of @mchamberland and increase the boundary tolerance (using the "boundary tolerance =" input for your CD geometry).

ojalaj commented 7 years ago

Thank you @mchamberland and @blakewalters for further debugging. As you may guess, this is just an oversimplified example of the geometry in which we originally encountered the error, but I'll let my PhD student to know this workaround and for now, let's hope it works for our case.

ftessier commented 3 years ago

Does anyone know if this issue has been resolved through "recent" commit dealing with boundary tolerance and CD geometry errors since 2017? For example: #383, #388, fd008c25c2a33ac86fdca2946dabbca8edc86b8a.