From e5d1b88ec7d9099aeab19a7ff3f76e1af18aca20 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 2 Mar 2022 10:12:10 -0600 Subject: [PATCH] Adding support for generating a WeightWindows class from a wwinp file. --- openmc/weight_windows.py | 160 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 159 insertions(+), 1 deletion(-) diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index c10bb7bf41..74019a6ace 100644 --- a/openmc/weight_windows.py +++ b/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 \ No newline at end of file