Escaping Neutron with USRMED.F Reflection

Adobe_Achromat · 21 April 2026 23:50

Versions

Please provide the used software versions.

FLUKA: 4-5.1
Flair: 3.4-5.2
Operating system: Ubuntu on Windows Subsystem for Linux
Compiler: GNU Fortran (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0

Description

Hello, I am trying to incorporate specular reflecting boundary conditions to a micron scale geometry. I have been adapting Dávid Horváth’s usrmed.f script for neutron reflections to my purposes.

The geometry consists of a square region, CELLVOID, filled with 200 smaller cylinders (on the order of micron diameter). Collectively, we’ll refer to the CELLVOID and tungsten cylinders (named F01Y0101 to F01X1010) as the simulation “cell”. This cell is illuminated by a 440 GeV/c proton beam in the z direction, which has a uniform square shape with width about half that of the cell’s transverse directions.

The cell is contained in a sphere (the blackhole region), with an upstream and downstream part delineated by the faces ±z faces of the cubic cell. Between the ±z faces of the cell, the region outside the cell, extending in the X,Y directions is called the “boundary”. The idea is that any neutrons produced by interactions of the proton beam with the cylinders in the cell would be reflected back into the cell as soon as they tried to cross into the boundary region. In this way, the neutrons would only be able to exit downstream or upstream.

…Long story short, I’ve been having a lot of trouble with getting this to work perfectly. It seems that, somehow, a neutron is able to escape into the boundary region sometimes, and perhaps more strangely, will even have a step in the boundary region, whose material is vacuum. We’d expect the neutron to be transported directly to the bounding sphere of blackhole (correct me if I’m wrong on that part, but in any case, the neutron escaping into the boundary mysteriously is the bigger problem for me). I’ve made lots of edits and added bits to debug the usrmed.f script, as you can see, but I can’t seem to identify where the single escaping neutron is coming from, and I figured it would be good to ask for help at this point.

STEPSIZE is already set to be much lower than the problem dimensions, and the upstream and downstream regions (except for a small “vestibule” upstream for the primary proton beam) have their material set to blackhole so it can’t be making it way in by somehow scattering off something else. The script is set to reflect the neutron if it is coming from any of the regions directly in contact with the boundary region. I’ve also debugged by having USRMED raise errors in various cases to try to tell me where this neutron is coming from, but I can’t seem to catch it before it escapes and raises the MREG=NEWREG=BOUNDARY issue.

To replicate the problem: use the errorReport.flair file attached (the other .flair is just geometry). Compile a custom exe with usrmed_reflect.f (attached), and spawn 4 cores and run for at least 20,000 histories (or just use what is entered in the .flair project currently as the default run). The second spawn should crash and print the following in the .log file:

ERROR, Current Region is: 4
X,Y,Z = 3.3279637819895939E-003 -2.8039045997420905E-003 2.6820878649234598E-003

TX,TY,TZ = 0.48656139870009008 -0.86797836181848853 -9.9355768377609191E-002

STOP MREG and NEWREG equal in vacuum region

…With region 4 being the boundary region that has the reflection imposed by MATPROP and usrmed_reflect.f (and again, the material of that reguin is vacuum). The attached image, “escapedNeutron.png”, shows a cylindrical usrbin of neutron fluence that shows the simulation cell full of neutrons, but with a single neutron escaping into the boundary region. I’m trying to prevent that last part! The neutrons should only be able to escape upstream or downstream, ideally.

Input files

Please upload all relevant files. (FLUKA input file, Flair project file, user routines, and data files)

errorReport.flair → flair project that has all supporting cards for setting up the simulation

geometry.flair → just read for the simulation cell geometry (it’s generated by a python script, too monotonous to build by hand)

usrmed_reflect.f → the user routine that tries to make all neutrons reflect specularly if they try to escape the simulation cell in the X,Y directions. It’s a terribly mangled version of Dávid Horváth’s script for a neutron reflector, in short!

Thanks very much for your time and help, and hopefully I’m just overlooking something!

–Will

usrmed_reflect.f (4.9 KB)

geometry.flair (38.5 KB)

errorReport.flair (2.6 KB)

Adobe_Achromat · 5 May 2026 18:47

Hi, I understand that this is a really lengthy one, but I’ll just bump it once here in case that stirs any interest. Thanks!

kltaylor · 8 May 2026 14:03

Dear Will,

Apologies for the late reply, I am looking into this and will be get to you shortly

Thank you kindly for understanding!

dcalzola · 23 May 2026 09:55

Hello Will,

Sorry for the long delay in replying. I was able to reproduce your problem, and indeed I observed the crash. I strongly believe this is a numerical issue (e.g. neutrons generated exactly on the boundary, or small machine-precision errors introduced during the reflection procedure itself, which may place them incorrectly).

For this problem, I would suggest a few possible approaches:

Kill all non-neutron particles leaving the geometry of interest (and which are not reflected). This way, you can ensure that additional neutrons are not generated by decays occurring outside the geometry of interest. In my case, this already fixed the issue.
In the case of very rare spurious crashes, it could be useful to tolerate a limited number of numerical issues. For instance, you could raise a warning every time this issue occurs, but abort only after the fifth occurrence.

On a side note, gently shifting the reflected particles to ensure they fall inside the correct region can sometimes help (although not in this specific case):

            XX = XX + TXX * 1D-10
            YY = YY + TYY * 1D-10
            ZZ = ZZ + TZZ * 1D-10

Adobe_Achromat · 26 May 2026 19:39

Hello Daniele,

Thank you very much, I really appreciate your time and help with the problem. I can definitely see how this is numerical after what you pointed out. The upshot is that I was going to switch to all particles being reflected anyways (focusing on neutrons was just supposed to be an easy way to test the reflection script. That way if I score protons and neutrons, I see the neutrons are reflected back and the protons exit out of the cell freely. I neglected the fact that a small number of events could result in neutrons from decays in the boundary region!). So it’s possible that once I get rid of the IJ .eq. 8 requirement, it will solve it just like you noted in bulletpoint 1. And for bulletpoint 2, that’s a very good suggestion for the production runs. If it only happens very rarely, I really shouldn’t be worried about the result. Oh and for the gentle shift inward—I was even thinking of trying a spatial translation for the periodicity instead of a reflection eventually, and that would maybe allow the little ``nudge’’ to help with the problem too.

Thanks again for taking a look, I really appreciate the help and guidance!

--Will