Commit 65ab4f05 authored by Salko Jr, Robert's avatar Salko Jr, Robert
Browse files

Update Hdf5Tools after updating CTF HDF5 file format

Reformat the HDF5 file so it is more efficient.
parent 69ff526c
Loading
Loading
Loading
Loading
+106 −97
Original line number Diff line number Diff line
import h5py
import numpy as np
import sys
import os
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../../')))
from SubKit.utils.Timer import eTime

class Hdf5Interface(object):
    """ An interface to the SubKit formatted HDF5 file
@@ -26,7 +30,7 @@ class Hdf5Interface(object):
        me._crossSectionalArea = {} # Takes section number (1-based) and returns cross-sectional area
        me._pinNumRadial = {} # Takes 1-based pin index and returns number of radial nodes

        me.numSection = max(me.h5['/CORE/chan_sections'][()])
        me.numSection = max(me.h5['/CORE/chan_sections'][:])

        me.numState = 1
        # Get list of keys instead of constantly checking if dataset in h5 file because h5py contains is really slow
@@ -49,85 +53,80 @@ class Hdf5Interface(object):
        for state in range(1, me.numState+1):
           me._stateDatasets[state] = me.h5['STATE_{:04d}'.format(state)].keys()

        meshDatasets = sorted(me.h5['CORE/mesh'].keys())
        if 'PIN_' in meshDatasets[-1]:
           me.numPin = int(meshDatasets[-1][4:])
        if 'CORE/pin_numaxial' in me.h5:
            me.pinNumAxial = me.h5['CORE/pin_numaxial'][:]
            me.numPin = len(me.pinNumAxial)
            me.pinNumAxial = me.h5['CORE/pin_numaxial'][:]
            me.pinNumSector = me.h5['CORE/pin_num_azimuthal'][:]
            me.pinNumRadial = me.h5['CORE/pin_numradial'][:]
            me.pinAxialNodes = me.h5['CORE/pin_axialnodes'][:, :]
            me.pinAxialHeights = me.h5['CORE/pin_axialheights'][:, :]
            me.pinRadialNodes = me.h5['CORE/pin_radial_nodes'][:, :]
            me.pinChanConnections = me.h5['CORE/pin_chan_connections'][:, :, :]
            me.pinAzimuthalFraction = me.h5['CORE/pin_azimuthal_fraction'][:, :, :]
            pinMaxAxial = np.amax(me.pinNumAxial)
            pinMaxSect = np.amax(me.pinNumSector)
            pinMaxLevel = np.amax(me.pinNumAxial)
            pinMaxRadial = np.amax(me.pinNumRadial)
            me.surfMask = np.ones((pinMaxAxial, pinMaxSect, me.numPin), dtype=bool)
            me.surfMask[:, :, :] = True
            for pin in range(me.numPin):
               me.surfMask[0:me.pinNumAxial[pin], 0:me.pinNumSector[pin], pin] = False
            # Create a mask for pin internal data
            me.pinInternalMask = np.ones((pinMaxSect, pinMaxLevel, pinMaxRadial, me.numPin), dtype=bool)
            me.pinInternalMask[:, :, :, :] = True
            for pin in range(me.numPin):
               nlev = me.pinNumAxial[pin]
               nsec = me.pinNumSector[pin]
               nrad = me.pinNumRadial[pin]
               me.pinInternalMask[0:nsec, 0:nlev, 0:nrad, pin] = False
            # Create a mask for surface values from the pin_temp dataset
            me.pinTempSurfMask = np.ones((pinMaxSect, pinMaxLevel, pinMaxRadial, me.numPin), dtype=bool)
            me.pinTempSurfMask[:, :, :, :] = True
            for pin in range(me.numPin):
                nlev = me.pinNumAxial[pin]
                nsec = me.pinNumSector[pin]
                # Only the surface value is unmasked
                me.pinTempSurfMask[0:nsec, 0:nlev, -1, pin] = False
            # Create a mask for pin axial data
            me.pinAxialMask = np.ones((pinMaxLevel, me.numPin), dtype=bool)
            me.pinAxialMask[:, :] = True
            for pin in range(me.numPin):
                nlev = me.pinNumAxial[pin]
                me.pinAxialMask[0:nlev, pin] = False
        else:
            me.numPin = 0

        # Takes pin index (1-based) and returns the number of levels in that pin
        me.pinNumAxial = {}
        for pin in range(me.numPin):
            me.pinNumAxial[pin+1] = len(me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(pin+1)])

        # Takes pin index (1-based) and returns number of sectors
        me.pinNumSector = {}
        for pin in range(me.numPin):
           me.pinNumSector[pin+1] = me.h5['CORE/mesh/PIN_{:05d}/azimuthal_fraction'.format(pin+1)].shape[1]
        me.chanSections = me.h5['CORE/chan_sections'][:]

        me.numCh = len(me.chanSections)

        # Create a mask for pin surface datasets, which can be used to extract relevant data
        # True means the data is irrelevant
        # For now we assum that dnbr will always be in the file (if there are pins in the model)
        # If we ever make it optional, will need to do a search for a surf dataset
        if 'STATE_0001/pin_surf_dnbr' in me.h5:
           me.surfMask = np.ones(me.h5['STATE_0001/pin_surf_dnbr'].shape, dtype=bool)
           for pin in range(1, me.numPin+1):
              nlev = me.getPinNumLev(pin)
              nsec = me.getPinNumSector(pin)
              me.surfMask[0:nlev, 0:nsec, pin-1] = False
        me.sectionHeights = me.h5['CORE/section_heights'][:, :]
        me.sectionBounds = me.h5['CORE/section_bounds'][:, :]
        me.sectionLevelCenters = me.h5['CORE/section_level_centers'][:, :]
        me.sectionNumLevel = me.h5['CORE/section_num_level'][:]

        lastCh = meshDatasets[-me.numPin-1]
        assert 'CHAN_' in lastCh
        me.numCh = int(lastCh[5:])
        me.sectionBottom = np.array([0.0]*me.numSection)
        me.sectionTop = np.array([0.0]*me.numSection)
        for secID in range(me.numSection):
            me.sectionBottom[secID] = me.sectionBounds[0, secID]
            me.sectionTop[secID] = me.sectionBounds[me.sectionNumLevel[secID]-2, secID]

        # Create a mask for channel datasets
        # Need one for scalar data and one for momentum data
        me.chanMomMask = np.ones(me.h5['STATE_0001/chan_temp'].shape, dtype=bool)
        me.chanScalarMask = np.ones(me.h5['STATE_0001/chan_pressure'].shape, dtype=bool)
        for chan in range(1, me.numCh+1):
              nlev = len(me.h5['CORE/mesh/CHAN_{:05d}/scalar_bounds'.format(chan)])
              me.chanMomMask[0:nlev, chan-1] = False
              me.chanScalarMask[0:nlev+1, chan-1] = False

        # Create a mask for pin internal data
        if me.numPin>0:
           me.pinInternalMask = np.ones(me.h5['STATE_0001/pin_temp'].shape, dtype=bool)
           me.pinInternalMask[:, :, :, :] = True
           for pin in range(1, me.numPin+1):
              nlev = me.getPinNumLev(pin)
              nsec = me.getPinNumSector(pin)
              nrad = me.getPinNumRadial(pin)
              me.pinInternalMask[0:nsec, 0:nlev, 0:nrad, pin-1] = False

        # Create a mask for surface values from the pin_temp dataset
        if 'STATE_0001/pin_temp' in me.h5:
           me.pinTempSurfMask = np.ones(me.h5['STATE_0001/pin_temp'].shape, dtype=bool)
           me.pinTempSurfMask[:, :, :, :] = True
           for pin in range(1, me.numPin+1):
              nlev = me.getPinNumLev(pin)
              nsec = me.getPinNumSector(pin)
              # Only the surface value is unmasked
              me.pinTempSurfMask[0:nsec, 0:nlev, -1, pin-1] = False

        # Create a mask for pin axial data
        if 'STATE_0001/pin_axial_centerline_temp' in me.h5:
           me.pinAxialMask = np.ones(me.h5['STATE_0001/pin_axial_centerline_temp'].shape, dtype=bool)
           for pin in range(1, me.numPin+1):
              nlev = me.getPinNumLev(pin)
              me.pinAxialMask[0:nlev, pin-1] = False

        # Location the axial locations of the section boundaries
        me.sectionBounds = [0.0]
        for section in list(range(1, me.getNumSections()+1)):
            chans = me.getChansInSection(section)
            me.sectionBounds.append(me.getChanHeight(chans[0])+me.sectionBounds[section-1])
        me.chanMomMask = np.ones((np.amax(me.sectionNumLevel)-1, me.numCh), dtype=bool)
        me.chanScalarMask = np.ones((np.amax(me.sectionNumLevel), me.numCh), dtype=bool)
        for chan in range(me.numCh):
            secID = me.chanSections[chan]
            nlev = me.sectionNumLevel[secID-1]
            me.chanMomMask[0:nlev-1, chan-1] = False
            me.chanScalarMask[0:nlev, chan-1] = False

        # Get the mesh type for each channel dataset in the model
        me.chanDatasetMeshTypes = {}
        for state in range(1, me.numState+1):
            for datasetName in me.h5['STATE_{:04d}'.format(state)]:
                dataset = me.h5['STATE_{:04d}/{:s}'.format(state, datasetName)]
        for state in range(me.numState):
            for datasetName in me.h5['STATE_{:04d}'.format(state+1)]:
                dataset = me.h5['STATE_{:04d}/{:s}'.format(state+1, datasetName)]
                if 'mesh_type' in dataset.attrs:
                    me.chanDatasetMeshTypes[datasetName] = dataset.attrs['mesh_type'][0]

@@ -248,7 +247,7 @@ class Hdf5Interface(object):

        """
        me._pinIndexValid(pin)
        return me.pinNumAxial[pin]
        return me.pinNumAxial[pin-1]

    def getPinDz(me, pin):
        """ Returns a vector of the heights of each level in the pin mesh
@@ -258,7 +257,7 @@ class Hdf5Interface(object):

        """
        me._pinIndexValid(pin)
        return me.h5['/CORE/mesh/PIN_{:05d}/axial_heights'.format(pin)][()]
        return me.pinAxialHeights[0:me.pinNumAxial[pin-1], pin-1]

    def getPinZ(me, pin):
       """ Returns a vector of the axial locations of each level in the pin mesh
@@ -290,7 +289,7 @@ class Hdf5Interface(object):
        """
        me._pinLevelValid(pin, level)
        section = me.getPinSection(pin, level)
        return me.h5['CORE/mesh/PIN_{:05d}/azimuthal_fraction'.format(pin)][section-1, :][()]
        return me.pinAzimuthalFraction[:, section-1, pin-1]

    def getDatasetUnits(me, dataset):
        """ Returns the units for a passed dataset name in the file.
@@ -320,8 +319,15 @@ class Hdf5Interface(object):
        """
        me._chanIndexValid(chan)
        me._datasetValid(dataset)
        secID = me.chanSections[chan-1]
        numLev = me.sectionNumLevel[secID-1]
        meshType = me.chanDatasetMeshTypes[dataset]
        return len(me.h5['CORE/mesh/CHAN_{:05d}/{:s}'.format(chan, meshType)])
        if meshType=='scalar_centers':
            return numLev
        elif meshType=='scalar_bounds':
            return numLev-1
        else:
            assert False

    def getChanAxialMesh(me, chan, dataset):
        """ Gets the axial mesh (a vector of axial locations) that is applicable for passed chan/dataset
@@ -332,7 +338,11 @@ class Hdf5Interface(object):

        """
        meshType = me.chanDatasetMeshTypes[dataset]
        return me.h5['CORE/mesh/CHAN_{:05d}/{:s}'.format(chan, meshType)][()]
        secID = me.chanSections[chan-1]
        if meshType=='scalar_centers':
            return me.sectionLevelCenters[0:me.sectionNumLevel[secID-1], secID-1]
        else:
            return me.sectionBounds[0:me.sectionNumLevel[secID-1]-1, secID-1]

    def getPinNumRadial(me, pin):
        """ Returns the number of radial nodes in the passed pin index
@@ -342,9 +352,7 @@ class Hdf5Interface(object):

        """
        me._pinIndexValid(pin)
        if pin not in me._pinNumRadial:
           me._pinNumRadial[pin] = len(me.h5['CORE/mesh/PIN_{:05d}/radial_nodes'.format(pin)])
        return me._pinNumRadial[pin]
        return me.pinNumRadial[pin-1]

    def getPinRadialNodes(me, pin):
        """ Returns the radial node locations in the pin
@@ -354,7 +362,7 @@ class Hdf5Interface(object):

        """
        me._pinIndexValid(pin)
        return me.h5['CORE/mesh/PIN_{:05d}/radial_nodes'.format(pin)][()]
        return me.pinRadialNodes[0:me.pinNumRadial[pin-1], pin-1]

    def getPinNumSector(me, pin):
        """ Returns the number of azimuthal sectors in passed pin index
@@ -364,7 +372,7 @@ class Hdf5Interface(object):

        """
        me._pinIndexValid(pin)
        return me.pinNumSector[pin]
        return me.pinNumSector[pin-1]

    def _allChanSameLevels(me):
        """ Returns True if all channels in the model have the same number of levels."""
@@ -533,7 +541,7 @@ class Hdf5Interface(object):
        me._stateValid(state)
        me._datasetValid(dataset)
        me._pinIndexValid(pin)
        levels = list(range(me.pinNumAxial[pin]))
        levels = list(range(me.pinNumAxial[pin-1]))
        return me.h5['/STATE_{:04d}/{:s}'.format(state, dataset)][levels, pin-1][()]

    def getPinSurfaceDataset(me, dataset, state, pin):
@@ -550,7 +558,7 @@ class Hdf5Interface(object):
        me._stateValid(state)
        me._datasetValid(dataset)
        me._pinIndexValid(pin)
        n = me.pinNumAxial[pin]
        n = me.pinNumAxial[pin-1]
        s = me.getPinNumSector(pin)
        return me.h5['/STATE_{:04d}/{:s}'.format(state, dataset)][0:n, 0:s, pin-1][()]

@@ -569,7 +577,7 @@ class Hdf5Interface(object):

        me._stateValid(state)
        me._pinIndexValid(pin)
        n = me.pinNumAxial[pin]
        n = me.pinNumAxial[pin-1]
        s = me.getPinNumSector(pin)
        return me.h5['/STATE_{:04d}/pin_temp'][0:s, 0:n, -1, pin-1][()]

@@ -832,7 +840,7 @@ class Hdf5Interface(object):
        """
        me._pinIndexValid(pin)
        me._pinLevelValid(pin, level)
        return me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(pin)][()][level]
        return me.pinAxialNodes[level, pin-1]

    def getRadialTemp(me, state, pin, level, sector=None, includeCenterline=False):
        """ Returns the radial temperature distribution across the pin radius
@@ -951,9 +959,9 @@ class Hdf5Interface(object):
        me._pinLevelValid(pin, level)
        if me.getNumSections()==1: return 1
        z = me.getPinAxial(pin, level)
        for i in range(len(me.sectionBounds)-1):
            if z>=me.sectionBounds[i] and z<me.sectionBounds[i+1]:
                section = i+1
        for secID in range(me.numSection):
            if z>=me.sectionBottom[secID] and z<me.sectionTop[secID]:
                section = secID+1
        return section

    def getChanSection(me, chan):
@@ -989,9 +997,8 @@ class Hdf5Interface(object):
        me._pinIndexValid(pin)
        me._pinLevelValid(pin, level)
        me._pinSectorValid(pin, sector)
        path = "CORE/mesh/PIN_{:05d}/chan_index".format(pin)
        section = me.getPinSection(pin, level)
        return me.h5[path][section-1, sector]
        return me.pinChanConnections[sector, section-1, pin-1]

    def getPinWithMinDNBR(me, state=None):
        """ Returns state, pin, sector index, level index experiencing minimum DNBR in the model
@@ -1045,7 +1052,7 @@ class Hdf5Interface(object):
            _sector = sector

        data = me.getPinSurfaceDataset('pin_surf_dnbr', _state, _pin)
        axial = me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(_pin)][()]
        axial = me.pinAxialNodes[0:me.pinNumAxial[_pin-1], _pin-1]

        return axial, data[:, _sector]

@@ -1064,7 +1071,7 @@ class Hdf5Interface(object):
        me._stateValid(state)
        me._pinIndexValid(pin)
        me._pinSectorValid(pin, sector)
        axial = me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(pin)][()]
        axial = me.pinAxialNodes[0:me.pinNumAxial[pin-1], pin-1]
        data = me.getPinSurfaceDataset('pin_surf_q_tot', state, pin)
        return axial, data[:, sector]

@@ -1127,7 +1134,7 @@ class Hdf5Interface(object):
        if sector is not None: _sector = sector

        data = me.h5['STATE_{:04d}/pin_temp'.format(_state)][()][_sector, 0:me.getPinNumLev(_pin), me.getPinNumRadial(_pin)-1, _pin-1]
        axial = me.h5['CORE/mesh/PIN_{:05d}/axial_nodes'.format(_pin)][()]
        axial = me.pinAxialNodes[0:me.pinNumAxial[_pin-1], _pin-1]

        return list(axial), list(data)

@@ -1139,7 +1146,8 @@ class Hdf5Interface(object):

        """
        me._chanIndexValid(chan)
        return me.h5['CORE/mesh/CHAN_{:05d}/scalar_centers'.format(chan)][()]
        secID = me.chanSections[chan-1]
        return me.sectionLevelCenters[0:me.sectionNumLevel[secID-1], secID-1]


    def getChanLiqEnth(me, chan, state=1):
@@ -1304,7 +1312,8 @@ class Hdf5Interface(object):
           chan (int) : Channel index (1-based)

        """
        return np.sum(me.h5['CORE/mesh/CHAN_{:05d}/scalar_heights'.format(chan)][()])
        secID = me.chanSections[chan-1]
        return np.sum(me.sectionHeights[0:me.sectionNumLevel[secID-1]-2, secID-1])

    def _getPinStart(me, pin):
        """ Gets the global axial location at the bottom of the pin
+7.17 KiB (1.24 MiB)

File changed.

No diff preview for this file type.

−360 B (96.1 KiB)

File changed.

No diff preview for this file type.

+3.49 KiB (133 KiB)

File changed.

No diff preview for this file type.