FLUKA Experts,
I am attempting to simulate planetary albedo radiation by impinging a spherical body with the GCR spectrum (with no atmosphere) and then examining the particles that return from the planetary body to space.
I am running into three issues with which I would greatly appreciate your assistance.
The albedo photons do not appear to produce well-defined peaks that I would expect to see based on the reactions occurring in the regolith. Is there some other way that I should be configuring my physics settings to calculate the expected gammas?
My input results in the “geofar” error when certain particles are transported across the planetary surface. You have indicated in previous responses that this error is generally traceable to geometry errors; however, my geometry is quite simple in this case. Is there an issue with the precision of the geometry due to the size of the planetary body? I have not been able to eliminate these errors while maintaining an acceptable size for the problem I am attempting to solve.
Here is an example of the error I am receiving below:
Geofar: Particle in region 1 (cell # 0) in position -5.246185808E+05 -4.504638359E+05 1.089533346E+05 is now causing trouble, requesting a step of 1.378900290E+00 cm to direction -2.992359831E-01 -5.294804002E-01 7.937936332E-01 end position -5.246189934E+05 -4.504645660E+05 1.089544291E+05 R2: 6.914783603E+05 R3: 7.000093942E+05 cm error count: 0
XU (2D): 3.954965289E+05 XU (3D): 4.819829922E+05 cm
XUOLD(2D): 1.317632191E+05 XUOLD(3D): 7.253226538E+04 cm
Kloop: 182379, Irsave: 1, Irsav2: 1, error code: -3 Nfrom: 5000
old direction -6.770690511E-01 4.960215844E-01 -5.436359883E-01, lagain, lstnew, lsense, lsnsct F F F T
Particle index 7 total energy 9.941665052E-05 GeV Nsurf 0
We succeeded in saving the particle: current region is n. 2 (cell # 0)
I have attached a simplified version of my input file for reference. Thank you in advance for your help and expertise. albedo-0465.inp (9.8 KB)
Sorry for the late reply, and thank you for providing additional material concerning your question.
Deuterons
In order to have more accurate results for deuterons, you need to activate the following physics card : SDUM = COALESCE, otherwise, you do not produce deuterons at all above 50MeV. Once you have activated them, you will see that your yield will be higher for the 50 - 500 MeV/n range. Please note that deuterons interactions below 150 MeV/n are not yet available in FLUKA (as of version 4.1.1). The IONSPLIT physics card can be used for a coarse attempt at reproducing them.
Geofar Issue
I ran into geofar issues like you while running your input file. As you saw, those are not fatal and I would posit that they are caused by numerical issues when manipulating large numbers. Taking into account the size of the region (10^6 cm) and the typical step length you encounter (10^-3 cm), and comparing with the 15 stable decimal places of a double, we can expect numerical trouble when squaring and subtracting
Gamma spectrum
The well-defined peaks seen in the gamma spectrum are coming from (low-energy-)neutron capture events. FLUKA follows a group-wise treatment for such interactions. This means that the secondary gammas production is performed group-wise as well (details about gamma groups are available in the manual). Sharp peaks at certain energies are therefore spread out uniformly across the width of the corresponding group → we observe fatter and lower peaks
To illustrate this, please look at the figure below showing, at production, the gamma yield coming off of neutron projectiles. In orange are the photons coming from low-energy neutron capture and in red, the ones produced as nuclear inelastic secondaries. We can verify that the groups with the highest intensity are indeed the ones where the most important lines are.
Now that is a picture at production level. At the location of your scoring (exit of the planet region) the gamma spectrum is composed of annihilation photons, a background coming most likely from Compton scattering at those energies, and the photons from the spectrum above or at least those that make it out, considering their mean path length is of the order of tens of cm at these energies.
In the case of a group-wise treatment, the peaks are not visible on the exit spectrum because they are flatter at production and therefore cannot be resolved at exit.
We checked that last point by running a simulation where we modify the energy of produced photons after neutron capture. The energy of photons in the 7.5 → 8 MeV group is set at 7.638 MeV, tentatively reproducing the double Iron line around this energy. In that case, the scoring at exit indeed shows the peak emerging from the background.
