periodic_surface_id is lost when reading a statepoint to perform a restart simulation

[ecasglez]

Bug Description

I have run some simulations that make use of 4 surfaces with periodic boundary conditions, so I have used periodic_surface to define the periodic surfaces. The geometry.xml file is properly generated including periodic_surface_id with the id of the periodic surfaces, and openmc runs well.

I am trying to do a restart starting from the results of those simulations but the periodic_surface_id is not read from the statepoint file for some reason. A new geometry.xml is generated but it does not contain the id of the periodic surfaces and the code stops.

I do not know if this parameter is not read from the statepoint when trying to restart, or if it is missing from the file. An example that illustrates this behavior is given below.

Steps to Reproduce

This input reproduces the bug. First, it configures a simple case with periodic boundary conditions, generates geometry.xml and runs openmc generating a statepoint file. Then it loads this statepoint file and exports the read geometry to geometry_after_reading_statepoint.xml file. The differences are appreciated while comparing those files as shown below:

import openmc

#Materials
mat_uo2 = openmc.Material()
mat_uo2.add_element('U', 1.0, 'ao', enrichment=4.35, enrichment_type='wo')
mat_uo2.add_element('O', 2.0, 'ao')
mat_uo2.set_density('g/cm3', 10.5156)
mat_uo2.temperature = 900
mats = openmc.Materials([mat_uo2])
mats.export_to_xml()

#Geometry
min_x = openmc.XPlane(x0=-11.5, boundary_type = 'periodic')
max_x = openmc.XPlane(x0= 11.5, boundary_type = 'periodic')
min_x.periodic_surface = max_x
min_y = openmc.YPlane(y0=-11.5, boundary_type = 'periodic')
max_y = openmc.YPlane(y0= 11.5, boundary_type = 'periodic')
min_y.periodic_surface = max_y

pincell = openmc.Cell(name='pincell')
pincell.region = +min_x & -max_x & +min_y & -max_y
pincell.fill = mat_uo2

pinuniv = openmc.Universe(name='pincell', cells=[pincell])

geom = openmc.Geometry()
geom.root_universe = pinuniv
geom.export_to_xml()

#Settings
point = openmc.stats.Point((0, 0, 0))
src = openmc.IndependentSource(space=point)
settings = openmc.Settings()
settings.source = src
settings.batches = 10     # Number of generations
settings.inactive = 2
settings.particles = 500 # Number of particles per generation

settings.export_to_xml()

openmc.run()

#Restart loading
statepoint = 'statepoint.10.h5'
with openmc.StatePoint(statepoint) as sp:
    geometry = sp.summary.geometry
geometry.export_to_xml('geometry_after_reading_statepoint.xml')

diff geometry.xml geometry_after_reading_statepoint.xml
3,7c3,7
<   <cell id="1" material="1" name="pincell" region="1 -2 3 -4" universe="1"/>
<   <surface boundary="periodic" coeffs="-11.5" id="1" periodic_surface_id="2" type="x-plane"/>
<   <surface boundary="periodic" coeffs="11.5" id="2" periodic_surface_id="1" type="x-plane"/>
<   <surface boundary="periodic" coeffs="-11.5" id="3" periodic_surface_id="4" type="y-plane"/>
<   <surface boundary="periodic" coeffs="11.5" id="4" periodic_surface_id="3" type="y-plane"/>
---
>   <cell id="1" material="1" name="pincell" region="1 -2 3 -4" temperature="900.0" universe="1"/>
>   <surface boundary="periodic" coeffs="-11.5" id="1" type="x-plane"/>
>   <surface boundary="periodic" coeffs="11.5" id="2" type="x-plane"/>
>   <surface boundary="periodic" coeffs="-11.5" id="3" type="y-plane"/>
>   <surface boundary="periodic" coeffs="11.5" id="4" type="y-plane"/>

Environment

OpenMC 0.14 compiled from sources in a Linux computer.