Adding support for generating a WeightWindows class from a wwinp file.

openmc/weight_windows.py

@@ -5,7 +5,7 @@ from xml.etree import ElementTree as ET
 import numpy as np
 
 from openmc.filter import _PARTICLES
-from openmc.mesh import MeshBase, UnstructuredMesh
+from openmc.mesh import MeshBase, RectilinearMesh, UnstructuredMesh
 import openmc.checkvalue as cv
 
 from ._xml import get_text
@@ -410,3 +410,161 @@ class WeightWindows(IDManagerMixin):
             weight_cutoff=weight_cutoff,
             id=id
         )
+
+    @staticmethod
+    def wwinp(filename):
+        """
+        Returns the next value in the wwinp file.
+
+        filename : str or pathlib.Path
+            Location of the wwinp file
+        """
+        fh = open(filename, 'r')
+
+        # read the first line of the file and
+        # keep only the first four entries
+        while(True):
+            line = next(fh)
+            if line and not line.startswith('c'):
+                break
+
+        values = line.strip().split()[:4]
+        for value in values:
+            yield value
+
+        # the remainder of the file can be read as
+        # sequential values
+        while(True):
+            line = next(fh)
+            # skip empty or commented lines
+            if not line or line.startswith('c'):
+                continue
+            values = line.strip().split()
+            for value in values:
+                yield value
+
+    @classmethod
+    def from_wwinp(cls, filename):
+        """Reads a wwinp file into WeightWindowDomain's
+
+        Parameters
+        ----------
+        path : str
+            Path to the wwinp file.
+
+        Returns
+        -------
+        list of openmc.WeightWindows
+        """
+        # create generator for getting the next parameter from the file
+        wwinp = WeightWindows.wwinp(filename)
+
+        # first parameter, if, of wwinp file is unused
+        next(wwinp)
+
+        # check time parameter, iv
+        if int(float(next(wwinp))) > 1:
+            raise ValueError('Time-dependent weight windows are not yet supported.')
+
+        # number of particles, ni
+        ni = int(float(next(wwinp)))
+
+        # read the mesh type, nr
+        nr = int(float(next(wwinp)))
+
+        if nr != 10:
+            # TODO: read the first entry by default and display a warning
+            raise ValueError('Cylindrical meshes are not currently supported')
+
+        # read the number of energy groups for each particle, ne
+        nes = [int(next(wwinp)) for _ in range(ni)]
+
+        if len(nes) == 1:
+            particles = ['neutron']
+        elif len(nes) == 2:
+            particles = ['neutron', 'photon']
+        else:
+            msg = ('More than two particle types are present. '
+                   'Only neutron and photon weight windows will be read.')
+            raise Warning(msg)
+
+        # read number of fine mesh elements in each coarse
+        # element: nfx, nfy, nfz
+        nfx = int(float(next(wwinp)))
+        nfy = int(float(next(wwinp)))
+        nfz = int(float(next(wwinp)))
+
+        # read the mesh origin: x0, y0, z0
+        llc = tuple(float(next(wwinp)) for _ in range(3))
+
+        # read the number of coarse mesh elements, ncx, ncy, ncz
+        ncx = int(float(next(wwinp)))
+        ncy = int(float(next(wwinp)))
+        ncz = int(float(next(wwinp)))
+
+        # skip the value defining the geometry type, nwg, we already know this
+        next(wwinp)
+
+        def _read_mesh_coords(wwinp, n_coarse_bins):
+            coords = [float(next(wwinp))]
+
+            for _ in range(n_coarse_bins):
+                # TODO: These are setup to read according to the MCNP5 format
+                sx = int(float(next(wwinp))) # number of fine mesh elements in between, sx
+                px = float(next(wwinp))  # value of next coordinate, px
+                qx = next(wwinp)  # this value is unused, qx
+                print(qx)
+
+                # append the fine mesh coordinates for this coarse element
+                coords += list(np.linspace(coords[-1], px, sx + 1))[1:]
+
+            return np.asarray(coords)
+
+        # read the coordinates for each dimension into a rectilinear mesh
+        mesh = RectilinearMesh()
+        mesh.x_grid = _read_mesh_coords(wwinp, ncx)
+        mesh.y_grid = _read_mesh_coords(wwinp, ncy)
+        mesh.z_grid = _read_mesh_coords(wwinp, ncz)
+
+        dims = ('x', 'y', 'z')
+        # check consistency of mesh coordinates
+        mesh_llc = mesh_val = (mesh.x_grid[0], mesh.y_grid[0], mesh.z_grid[0])
+        for dim, header_val, mesh_val in zip(dims, llc, mesh_llc):
+            if header_val != mesh_val:
+                msg = ('The {} corner of the mesh ({}) does not match '
+                       'the value read in block 1 of the wwinp file ({})')
+                raise ValueError(msg.format(dim, mesh_val, header_val))
+
+        mesh_dims = mesh.dimension
+        for dim, header_val, mesh_val in zip(dims, (nfx, nfy, nfz), mesh_dims):
+            if header_val != mesh_val:
+                msg = ('Total number of mesh elements read in the {} '
+                       'direction ({}) is inconsistent with the '
+                       'number read in block 1 of the wwinp file ({})')
+                raise ValueError(msg.format(dim, mesh_val, header_val))
+
+        # total number of fine mesh elements, nft
+        nft = nfx * nfy * nfz
+        # read energy bins and weight window values for each particle
+        wws = []
+        for particle, ne in zip(particles, nes):
+            # read energy
+            e_groups = np.asarray([float(next(wwinp)) for _ in range(ne)])
+
+            # adjust energy from MeV to eV
+            e_groups *= 1E6
+
+            # create an array for weight window lower bounds
+            ww_lb = np.zeros((ne, nft))
+            for e in range(ne):
+                ww_lb[e, :] = [float(next(wwinp)) for _ in range(nft)]
+
+            settings = WeightWindows(id=None,
+                                     mesh=mesh,
+                                     lower_ww_bounds=ww_lb.flatten(),
+                                     upper_bound_ratio=5.0,
+                                     energy_bins=e_groups,
+                                     particle_type=particle)
+            wws.append(settings)
+
+        return wws