Particle lost without warning/error or being killed, stalling the run

Hi All,

I was running a bioshield calculation when I noticed that my run was stalling on a single particle.
I’ve reduced the problem to a minimal working example where the issue is still present (python script attached).

The MWE is a large void filled box with reflective boundary conditions at x=0 and y=0 (I had a symmetrical geometry), a 2.45 MeV monoenergetic isotropic neutron source and a single mesh tally element. The MWE stalls on particle 23228 – when printing the trace of that particle it enters the cell but then gets stuck without performing any interactions, crossing any surfaces, or dying.

I can’t print the particle track while the simulation is running to see where it is disappearing to (I can only view particle tracks when the simulation is finished, which doesn’t help when it has stalled) so I have no idea where or why the stall is occurring.

Does anyone have any ideas/suggestions on what’s causing this bug/how to fix it?

Kind Regards
Matt

My OpenMC version specs are:
OpenMC version 0.13.3
Git SHA1: 50e39a4e20dc9e0f3d7ccf07333f6a5e6c797c8c
Copyright (c) 2011-2023 MIT, UChicago Argonne LLC, and contributors
MIT/X license at https://docs.openmc.org/en/latest/license.html
Build type: Release
Compiler ID: GNU 11.3.0
MPI enabled: no
Parallel HDF5 enabled: no
PNG support: yes
DAGMC support: yes
libMesh support: no
MCPL support: no
NCrystal support: no
Coverage testing: no
Profiling flags: no
Bioshield.py (1.9 KB)

1 Like

Hi Matt,

I think you’ve found a problem with our regular mesh tally routine for tracklength tallies. I say this because if I either a) change the mesh dimensions slightly or b) change the tally estimator from “tracklength” (default) to “collision”, then the simulation runs successfully.

This is a great minimal working example of the bug and I’ll take a look at a fix very soon, but in the meantime maybe this is enough information to proceed with your work. I’ll make sure to link the fix here once it’s been incorporated into OpenMC.

-Patrick

Hi Patrick

Thanks for having a look at this. I think it might also have something to do with an odd interaction with the reflective boundary conditions. I’ve implemented a (less than ideal) workaround of simulating my entire geometry without reflective boundaries and this seems to work reliably. It’d be great to have a patch though :slight_smile:

Cheers,
Matt

1 Like

Hi Patrick,
I’m facing the same problem described by @MIB101.
However changing the estimator or the dimensions of the mesh didn’t help.
The simultation still “freezes”.
Is there any more insights about this bug?

Thanks

Hello

I am trying to track this problem down with openmc-0.15.0

In file mesh.cpp:669 there is the follwing loop:

       // For all directions outside the mesh, find the distance that we need to
       // travel to reach the next surface. Use the largest distance, as only
       // this will cross all outer surfaces.
       int k_max {0};
       for (int k = 0; k < n; ++k) {
         if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
             (distances[k].distance > traveled_distance)) {
           traveled_distance = distances[k].distance;
           k_max = k;
         }
       }

For particle 23228, traveled_distance is always equal to 0, so the particle will never advance in the subsequent code, causing the problem to hang because of an upstream infinite loop.

For particle 23228 the involved variables in this loop have the following values:

  • ijk = {-7, 1, 1}
  • shape_ = {1, 1, 1}
  • distances =
    {next_index = -6, max_surface = true, distance = 0},
    {next_index = 2, max_surface = true, distance = 1614.1689884883433},
    {next_index = 0, max_surface = false, distance = 4916.5287060087403}

So, in the case of particle 23228 the traveled_distance will be only updated when k==1, but the distance in this case is 0, so it won’t go into the if clause.

I am not sure if the problem is in ijk or in the distances variables.

distance is calculated in mesh.cpp:1050 as:

d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];

where:

  • positive_grid_boundary(ijk, i) == 100
  • r0[i] == 100

giving as a result a distance of 0.