EGSnrc
EGSnrc copied to clipboard
Dose deposited in vacuum in DOSRZnrc
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.
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:
- Am I missing something in this logic?
- 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.
- 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?
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.
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.