openmc-dev / openmc

OpenMC Monte Carlo Code
https://docs.openmc.org
Other
757 stars 486 forks source link

Mesh Tally Crashes the Simulation #3057

Open itay-space opened 3 months ago

itay-space commented 3 months ago

Bug Description

There is a serious bug in the calculations involving a mesh tally, specifically with the track-length estimator. The simulation freezes without producing any additional statepoints. This issue is evident because when it occurs, the number of CPUs in use drops from utilizing all available OpenMP threads to just one. Similar behavior has been reported multiple times in the OpenMC discussion forum: openmc-freezes-when-simulating-batches simulation-freezing particle-lost-without-warning-error-or-being-killed-stalling-the-run

Steps to Reproduce

Attached is a straightforward demonstration of this problem. The input file can be generated using the following Python code. This code will create a wall made of the H1 isotope and a mesh filter surrounding it. Additionally, there is an attached figure illustrating the simple geometry, with an arrow indicating the beam direction. In this example (with seed=17) the simulation will crash around the 6th batch.

image


H = openmc.Material(name='Hydrogen')
H.add_nuclide('H1', 1)
materials = openmc.Materials([H])
materials.export_to_xml()

z1 = openmc.ZPlane( z0 = 0)
z2 = openmc.ZPlane( z0 = 200)
cyl = openmc.ZCylinder(r = 200)
wall_reg = -cyl & +z1 & -z2
wall_cell = openmc.Cell(name = "wall" ,  region= wall_reg , fill = H)
world_sphere = openmc.Sphere(r=1000,name="world_sphere",boundary_type='vacuum')
world = openmc.Cell(region=-world_sphere &(~wall_reg) , name = 'world')
univ = openmc.Universe(cells=[wall_cell , world])
geometry = openmc.Geometry(univ)
geometry.export_to_xml()

settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.particles = 1000000
bat = 200
settings.batches = bat
settings.statepoint['batches'] = list(range(bat+1))[::1]
angle = openmc.stats.Monodirectional((0,0,1))
energy = openmc.stats.Discrete(10e6,1)
point = openmc.stats.Point((0, 0, -0.5))
source = openmc.Source(space=point,angle=angle,energy=energy,particle="neutron")
settings.source = source
settings.seed = 17
settings.export_to_xml()

mesh = openmc.RegularMesh()
mesh.dimension = [50, 50,50]
mesh.lower_left = [-200,-200,-10]
mesh.upper_right = [200,200,300]
mesh_filter = openmc.MeshFilter(mesh)
tallies = openmc.Tallies()
mesh_tally = openmc.Tally(name="mesh_tally")
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ["flux"]
mesh_tally.estimator = 'tracklength'
tallies.append(mesh_tally)
tallies.export_to_xml()

Environment

The tests were done using the latest docker image:

  "Id": "sha256:1b0e57ac2bdddda83097c9e3370c173f5dd9e92b4f4c818fc31d441a4dc25bbb",
  "Digest": null,
  "RepoDigests": [
    "openmc/openmc@sha256:56a36944da2cf5f2c2cb4054c8389344ec3068f36b8134af77c63e5f9abc6f57"
  ],
  "Labels": null
}

And some more details:

OpenMC version 0.14.0
Git SHA1: fa2330103de61a864c958d1a7250f11e5dd91468
Copyright (c) 2011-2023 MIT, UChicago Argonne LLC, and contributors
MIT/X license at <https://docs.openmc.org/en/latest/license.html>
Build type:            RelWithDebInfo
Compiler ID:           GNU 10.2.1
MPI enabled:           yes
Parallel HDF5 enabled: yes
PNG support:           yes
DAGMC support:         no
libMesh support:       no
MCPL support:          no
NCrystal support:      no
Coverage testing:      no
Profiling flags:       no

Data library: fendl-3.2-hdf5

pshriwise commented 3 months ago

Hi @itay-space! Thank you very much for reporting this and for the simple working example. I'll dig into it as soon as I can.

itay-space commented 2 months ago

Hi @pshriwise

I wanted to follow up on the issue I reported last week. Were you able to reproduce the problem?

Your assistance is greatly appreciated. Thank you for your help!

pshriwise commented 2 months ago

Hi @itay-space!

Thanks for your patience. I found the cause of the problem for the example you've provided above. The underlying issue was actually resolved in https://github.com/openmc-dev/openmc/pull/2811, which should be included in the latest OpenMC release. In this case, it was causing nan values in a particle position that led to an infinite loop in the mesh filter routine.

This doesn't seem to be the cause of all the hanging simulations you linked above, however, so I'm going to be looking into those as well separately, but I wanted to update you on this one sooner rather than later.