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

Dose deposited in vacuum in DOSRZnrc #1061

Open dworogers opened 1 year ago

dworogers commented 1 year ago

Using DOSRZnrc with a vacuum region on the surface, I occasionally get dose deposited in the vacuum region which is impossible. The attached file takes some time to run but shorter runs don't see the dose deposited. Blake looked into this in Oct 2022 and in a series of emails isolated the problem but the solution wasn't obvious. I can forward emails if not found at NRC.

vac_bug.txt

blakewalters commented 8 months ago

Looking back through our emails, @dworogers, I see that I'd traced this to electrons reflected back to the vacuum region and then immediately being discarded. Moreover, I was able to replicate this issue using dosxyznrc. After looking into it a bit more, I became fixated on this block of coding from subroutine ELECTR:

        IF( irnew = irl & eie <= ecut(irl)) [
           go to :ECUT-DISCARD:;
        ]

        medold = medium;
        IF(medium ~= 0)
        [
            ekeold = eke; eke = eie - rm; "update kinetic energy
            elke   = log(eke);
            $SET INTERVAL elke,eke; "Get updated interval
        ]

        IF(irnew ~= irold) [ $electron_region_change; ]

        "After transport call to user scoring routine

        $AUSCALL($TRANAUSA);

        IF(eie <= ecut(irl)) [
           go to :ECUT-DISCARD:;
        ]

The first condition:

            IF( irnew = irl & eie <= ecut(irl)) [
               go to :ECUT-DISCARD:;
            ]

only allows ECUT discard after the step if there has been no region change (i.e. irnew=irl). But then it seems the second condition:

            IF(eie <= ecut(irl)) [
               go to :ECUT-DISCARD:;
            ]

retroactively allows ECUT discard anyway, and energy deposition the next time through the loop, even if the region has been changed (in your example, to the vacuum region).

This gives rise to a couple of questions:

0) Am I missing something in this logic? 1) If E <= ECUT at the end of a charged particle step to region boundary, does it logically and physically make sense to dump the particle's energy in the new region? It seems like this is how we've been doing things all along. 2) If so, how should we deal with the case where this new region is vacuum? Should we give the region unit mass, as dosxynrc does, or should we just zero the dose? Should this be handled in the user code or in egsnrc.mortran?

ftessier commented 8 months ago

my take on question 2: Zero the does, and this ought to be low level in egsnrc.mortran. If the medium is vacuum, then any dose deposition ought to be identically 0.0, acrosss the board.

blakewalters commented 3 months ago

Okay, I'll go ahead and do this, @ftessier. As a result of this fix, though, some user codes--I'm looking at you, dosxyznrc--may need to be modified to handle (zero energy/zero mass) better, because (finite energy/zero mass) = Inf, which seems to be handled okay, while (zero energy/zero mass) = NAN, which causes a crash at the end of the run.