From dadc4fe4184e66848a60a052c00d0b4cfeebb247 Mon Sep 17 00:00:00 2001
From: GuySten <62616591+GuySten@users.noreply.github.com>
Date: Thu, 5 Jun 2025 03:07:59 +0300
Subject: [PATCH] Fix no serialization of periodic_surface_id bug (#3421)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 include/openmc/boundary_condition.h |  4 +++
 openmc/summary.py                   | 12 +++++++++
 src/surface.cpp                     | 13 +++++++++
 tests/unit_tests/test_summary.py    | 41 +++++++++++++++++++++++++++++
 4 files changed, 70 insertions(+)
 create mode 100644 tests/unit_tests/test_summary.py

diff --git a/include/openmc/boundary_condition.h b/include/openmc/boundary_condition.h
index cd8988162..af40131f1 100644
--- a/include/openmc/boundary_condition.h
+++ b/include/openmc/boundary_condition.h
@@ -111,6 +111,10 @@ public:
 
   std::string type() const override { return "periodic"; }
 
+  int i_surf() const { return i_surf_; }
+
+  int j_surf() const { return j_surf_; }
+
 protected:
   int i_surf_;
   int j_surf_;
diff --git a/openmc/summary.py b/openmc/summary.py
index 6423f8ce1..ca5cfadd7 100644
--- a/openmc/summary.py
+++ b/openmc/summary.py
@@ -127,11 +127,23 @@ class Summary:
             self._fast_materials[material.id] = material
 
     def _read_surfaces(self):
+        periodic_surface_ids = set()
         for group in self._f['geometry/surfaces'].values():
             surface = openmc.Surface.from_hdf5(group)
             # surface may be None for DAGMC surfaces
             if surface:
                 self._fast_surfaces[surface.id] = surface
+                if surface.boundary_type == "periodic":
+                    periodic_surface_ids.add(surface.id)
+
+        # Assign periodic surfaces when information is in file
+        for surface_id in periodic_surface_ids:
+            group = self._f[f'geometry/surfaces/surface {surface_id}']
+            surface = self._fast_surfaces[surface_id]
+            if 'periodic_surface_id' in group:
+                periodic_surface_id = int(group['periodic_surface_id'][()])
+                surface.periodic_surface = self._fast_surfaces[periodic_surface_id]
+
 
     def _read_cells(self):
 
diff --git a/src/surface.cpp b/src/surface.cpp
index 0002275c3..ea19514a3 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -172,6 +172,19 @@ void Surface::to_hdf5(hid_t group_id) const
     if (bc_) {
       write_string(surf_group, "boundary_type", bc_->type(), false);
       bc_->to_hdf5(surf_group);
+
+      // write periodic surface ID
+      if (bc_->type() == "periodic") {
+        auto pbc = dynamic_cast<PeriodicBC*>(bc_.get());
+        Surface& surf1 {*model::surfaces[pbc->i_surf()]};
+        Surface& surf2 {*model::surfaces[pbc->j_surf()]};
+
+        if (id_ == surf1.id_) {
+          write_dataset(surf_group, "periodic_surface_id", surf2.id_);
+        } else {
+          write_dataset(surf_group, "periodic_surface_id", surf1.id_);
+        }
+      }
     } else {
       write_string(surf_group, "boundary_type", "transmission", false);
     }
diff --git a/tests/unit_tests/test_summary.py b/tests/unit_tests/test_summary.py
new file mode 100644
index 000000000..315ac0ba5
--- /dev/null
+++ b/tests/unit_tests/test_summary.py
@@ -0,0 +1,41 @@
+import openmc
+
+
+def test_periodic_surface_roundtrip(run_in_tmpdir):
+    # Create a simple model with periodic surfaces
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+    mat.set_density('g/cm3', 1.0)
+    cyl = openmc.ZCylinder(r=1.0)
+    x0 = openmc.XPlane(-5.0, boundary_type='periodic')
+    y0 = openmc.YPlane(-5.0, boundary_type='periodic')
+    z0 = openmc.ZPlane(-5.0, boundary_type='periodic')
+    x1 = openmc.XPlane(5.0, boundary_type='periodic')
+    y1 = openmc.YPlane(5.0, boundary_type='periodic')
+    z1 = openmc.ZPlane(5.0, boundary_type='periodic')
+    x0.periodic_surface = x1
+    y0.periodic_surface = y1
+    z0.periodic_surface = z1
+    cell1 = openmc.Cell(fill=mat, region=-cyl)
+    cell2 = openmc.Cell(fill=mat, region=+cyl & +x0 & -x1 & +y0 & -y1 & +z0 & -z1)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2])
+    model.settings.particles = 100
+    model.settings.batches = 1
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(1.0e4)
+    )
+
+    # Run model
+    model.run()
+
+    # Load summary data and check periodic surfaces
+    summary = openmc.Summary('summary.h5')
+    surfs = summary.geometry.get_all_surfaces()
+    for s in [x0, y0, z0, x1, y1, z1]:
+        assert surfs[s.id].boundary_type == 'periodic'
+    pairs = [(x0, x1), (y0, y1), (z0, z1)]
+    for s0, s1 in pairs:
+        assert surfs[s0.id].periodic_surface == surfs[s1.id]
+        assert surfs[s1.id].periodic_surface == surfs[s0.id]