Simulations.py

    def cp(self, power, wind, A, rho=1.225):
        return power / (0.5 * rho * A * wind * wind * wind)

    def shaft_power(self):
        """
        Return the mechanical shaft power based on the shaft torsional loading
        """
        try:
            i = self.res.ch_dict['bearing-shaft_rot-angle_speed-rpm']['chi']
            rads = self.res.sig[:,i]*np.pi/30.0
        except KeyError:
            try:
                i = self.res.ch_dict['bearing-shaft_rot-angle_speed-rads']['chi']
                rads = self.res.sig[:,i]
            except KeyError:
                i = self.res.ch_dict['Omega']['chi']
                rads = self.res.sig[:,i]
        try:
            nn_shaft = self.config['nn_shaft']
        except:
            nn_shaft = 4
        itorque = self.res.ch_dict['shaft-shaft-node-%3.3i-momentvec-z'%nn_shaft]['chi']
        torque = self.res.sig[:,itorque]
        # negative means power is being extracted, which is exactly what a wind
        # turbine is about, we call that positive
        return -1.0*torque*rads

    def calc_torque_const(self, save=False, name='ojf'):
        """
        If we have constant RPM over the simulation, calculate the torque
        constant. The current loaded HAWC2 case is considered. Consequently,
        first load a result file with load_result_file

        Parameters
        ----------

        save : boolean, default=False

        name : str, default='ojf'
            File name of the torque constant result. Default to using the
            ojf case name. If set to hawc2, it will the case_id. In both
            cases the file name will be extended with '.kgen'

        Returns
        -------

        [windspeed, rpm, K] : list

        """
        # make sure the results have been loaded previously
        try:
            # get the relevant index to the wanted channels
            # tag: coord-bodyname-pos-sensortype-component
            tag = 'bearing-shaft_nacelle-angle_speed-rpm'
            irpm = self.res.ch_dict[tag]['chi']
            chi_rads = self.res.ch_dict['Omega']['chi']
            tag = 'shaft-shaft-node-001-momentvec-z'
            chi_q = self.res.ch_dict[tag]['chi']
        except AttributeError:
            msg = 'load results first with Cases.load_result_file()'
            raise ValueError(msg)

#        if not self.case['[fix_rpm]']:
#            print
#            return

        windspeed = self.case['[windspeed]']
        rpm = self.res.sig[:,irpm].mean()
        # and get the average rotor torque applied to maintain
        # constant rotor speed
        K = -np.mean(self.res.sig[:,chi_q]*1000./self.res.sig[:,chi_rads])

        result = np.array([windspeed, rpm, K])

        # optionally, save the values and give the case name as file name
        if save:
            fpath = self.case['[post_dir]'] + 'torque_constant/'
            if name == 'hawc2':
                fname = self.case['[case_id]'] + '.kgen'
            elif name == 'ojf':
                fname = self.case['[ojf_case]'] + '.kgen'
            else:
                raise ValueError('name should be either ojf or hawc2')
            # create the torque_constant dir if it doesn't exists
            try:
                os.makedirs(fpath)
            except OSError:
                pass

#            print('gen K saving at:', fpath+fname
            np.savetxt(fpath+fname, result, header='windspeed, rpm, K')

        return result

    def compute_envelope(self, sig, ch_list, int_env=False, Nx=300):
        """
        The function computes load envelopes for given signals and a single
        load case. Starting from Mx and My moments, the other cross-sectional
        forces are identified.

        Parameters
        ----------

        sig : list, time-series signal

        ch_list : list, list of channels for enevelope computation

        int_env : boolean, default=False
            If the logic parameter is True, the function will interpolate the
            envelope on a given number of points

        Nx : int, default=300
            Number of points for the envelope interpolation

        Returns
        -------

        envelope : dictionary,
            The dictionary has entries refered to the channels selected.
            Inside the dictonary under each entry there is a matrix with 6
            columns, each for the sectional forces and moments

        """

        envelope= {}
        for ch in ch_list:
            chi0 = self.res.ch_dict[ch[0]]['chi']
            chi1 = self.res.ch_dict[ch[1]]['chi']
            s0 = np.array(sig[:, chi0]).reshape(-1, 1)
            s1 = np.array(sig[:, chi1]).reshape(-1, 1)
            # Compute a Convex Hull, the vertices number varies according to
            # the shape of the poligon
            cloud =  np.append(s0, s1, axis=1)
            hull = scipy.spatial.ConvexHull(cloud)
            closed_contour = np.append(cloud[hull.vertices,:],
                                       cloud[hull.vertices[0],:].reshape(1,2),
                                       axis=0)

            # Interpolate envelope for a given number of points
            if int_env:
                _,_,_,closed_contour_int = int_envelope(closed_contour[:,0],
                                                        closed_contour[:,1],Nx)


            # Based on Mx and My envelope, the other cross-sectional moments
            # and forces components are identified and appended to the initial
            # envelope
            for ich in range(2, len(ch)):
                chix = self.res.ch_dict[ch[ich]]['chi']
                s0 = np.array(sig[hull.vertices, chix]).reshape(-1, 1)
                s1 = np.array(sig[hull.vertices[0], chix]).reshape(-1, 1)
                s0 = np.append(s0, s1, axis=0)
                closed_contour = np.append(closed_contour, s0, axis=1)
                if int_env:
                    _,_,_,extra_sensor = int_envelope(closed_contour[:,0],
                                                       closed_contour[:,ich],Nx)
                    es = np.atleast_2d(np.array(extra_sensor[:,1])).T
                    closed_contour_int = np.append(closed_contour_int,es,axis=1)

            if int_env:
                envelope[ch[0]] = closed_contour_int
            else:
                envelope[ch[0]] = closed_contour

        return envelope

    def envelope(self, silent=False, ch_list=[], append=''):
        """
        Calculate envelopes and save them in a table.

        Parameters
        ----------


        Returns
        -------


        """
        # get some basic parameters required to calculate statistics
        try:
            case = list(self.cases.keys())[0]
        except IndexError:
            print('no cases to select so no statistics, aborting ...')
            return None

        post_dir = self.cases[case]['[post_dir]']
        sim_id = self.cases[case]['[sim_id]']

        if not silent:
            nrcases = len(self.cases)
            print('='*79)
            print('statistics for %s, nr cases: %i' % (sim_id, nrcases))

        fname = os.path.join(post_dir, sim_id + '_envelope' + append + '.h5')
        h5f = tbl.openFile(fname, mode="w", title=str(sim_id),
                           filters=tbl.Filters(complevel=9))

        # Create a new group under "/" (root)
        for ii, (cname, case) in enumerate(self.cases.items()):

            groupname = str(cname[:-4])
            groupname = groupname.replace('-', '_')
            h5f.createGroup("/", groupname)
            ctab = getattr(h5f.root, groupname)

            if not silent:
                pc = '%6.2f' % (float(ii)*100.0/float(nrcases))
                pc += ' %'
                print('envelope progress: %4i/%i %s' % (ii, nrcases, pc))

            self.load_result_file(case)

            envelope = self.compute_envelope(self.sig, ch_list)

            for ch_id in ch_list:
                h5f.createTable(ctab, str(ch_id[0].replace('-', '_')),
                                EnvelopeClass.section,
                                title=str(ch_id[0].replace('-', '_')))
                csv_table = getattr(ctab, str(ch_id[0].replace('-', '_')))
                tablerow = csv_table.row
                for row in envelope[ch_id[0]]:
                    tablerow['Mx'] = float(row[0])
                    tablerow['My'] = float(row[1])
                    if len(row)>2:
                        tablerow['Mz'] = float(row[2])
                        if len(row)>3:
                            tablerow['Fx'] = float(row[3])
                            tablerow['Fy'] = float(row[4])
                            tablerow['Fz'] = float(row[5])
                        else:
                            tablerow['Fx'] = 0.0
                            tablerow['Fy'] = 0.0
                            tablerow['Fz'] = 0.0
                    else:
                        tablerow['Mz'] = 0.0
                        tablerow['Fx'] = 0.0
                        tablerow['Fy'] = 0.0
                        tablerow['Fz'] = 0.0
                    tablerow.append()
                csv_table.flush()
        h5f.close()

    def force_lower_case_id(self):
        """Keep for backwards compatibility with the dlctemplate.py
        """
        msg = "force_lower_case_id is depricated and is integrated in "
        msg += "Cases.createcase() instead."
        warnings.warn(msg, DeprecationWarning)

        tmp_cases = {}
        for cname, case in self.cases.items():
             tmp_cases[cname.lower()] = case.copy()
        self.cases = tmp_cases


class EnvelopeClass(object):
    """
    Class with the definition of the table for the envelope results
    """
    class section(tbl.IsDescription):

        Mx = tbl.Float32Col()
        My = tbl.Float32Col()
        Mz = tbl.Float32Col()
        Fx = tbl.Float32Col()
        Fy = tbl.Float32Col()
        Fz = tbl.Float32Col()


# TODO: implement this
class Results(object):
    """
    Move all Hawc2io to here? NO: this should be the wrapper, to interface
    the htc_dict with the io functions

    There should be a bare metal module/class for those who only want basic
    python support for HAWC2 result files and/or launching simulations.

    How to properly design this module? Change each class into a module? Or
    leave like this?
    """

    # OK, for now use this to do operations on HAWC2 results files

    def __init___(self):
        """
        """
        pass

    def m_equiv(self, st_arr, load, pos):
        r"""Centrifugal corrected equivalent moment

        Convert beam loading into a single equivalent bending moment. Note that
        this is dependent on the location in the cross section. Due to the
        way we measure the strain on the blade and how we did the calibration
        of those sensors.

        .. math::

            \epsilon = \frac{M_{x_{equiv}}y}{EI_{xx}} = \frac{M_x y}{EI_{xx}}
            + \frac{M_y x}{EI_{yy}} + \frac{F_z}{EA}

            M_{x_{equiv}} = M_x + \frac{I_{xx}}{I_{yy}} M_y \frac{x}{y}
            + \frac{I_{xx}}{Ay} F_z

        Parameters
        ----------

        st_arr : np.ndarray(19)
            Only one line of the st_arr is allowed and it should correspond
            to the correct radial position of the strain gauge.

        load : list(6)
            list containing the load time series of following components
            .. math:: load = F_x, F_y, F_z, M_x, M_y, M_z
            and where each component is an ndarray(m)

        pos : np.ndarray(2)
            x,y position wrt neutral axis in the cross section for which the
            equivalent load should be calculated

        Returns
        -------

        m_eq : ndarray(m)
            Equivalent load, see main title

        """

        F_z = load[2]
        M_x = load[3]
        M_y = load[4]

        x, y = pos[0], pos[1]

        A = st_arr[ModelData.st_headers.A]
        I_xx = st_arr[ModelData.st_headers.Ixx]
        I_yy = st_arr[ModelData.st_headers.Iyy]

        M_x_equiv = M_x + ( (I_xx/I_yy)*M_y*(x/y) ) + ( F_z*I_xx/(A*y) )
        # or ignore edgewise moment
        #M_x_equiv = M_x + ( F_z*I_xx/(A*y) )

        return M_x_equiv


class MannTurb64(prepost.PBSScript):
    """
    alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
    mann_turb_x64.exe fname 1.0 29.4 3.0 1209 256 32 32 2.0 5 5 true.

    Following tags have to be defined:
        * [tu_model]
        * [turb_base_name]
        * [MannAlfaEpsilon]
        * [MannL]
        * [MannGamma]
        * [tu_seed]
        * [turb_nr_u]
        * [turb_nr_v]
        * [turb_nr_w]
        * [turb_dx]
        * [turb_dy]
        * [turb_dz]
        * [high_freq_comp]
    """

    def __init__(self, silent=False):
        super(MannTurb64, self).__init__()
        self.exe = 'time wine mann_turb_x64.exe'
        self.winefix = 'winefix\n'
        # PBS configuration
        self.umask = '0003'
        self.walltime = '00:59:59'
        self.queue = 'workq'
        self.lnodes = '1'
        self.ppn = '1'
        self.silent = silent
        self.pbs_in_dir = 'pbs_in_turb/'

    def gen_pbs(self, cases):
        """
        Parameters
        ----------

        cases : dict of dicts
            each key holding a dictionary with tag/value pairs.
        """

        case0 = cases[list(cases.keys())[0]]
        # make sure the path's end with a trailing separator, why??
        self.pbsworkdir = os.path.join(case0['[run_dir]'], '')
        if not self.silent:
            print('\nStart creating PBS files for turbulence with Mann64...')
        for cname, case in cases.items():

            # only relevant for cases with turbulence
            if '[tu_model]' in case and int(case['[tu_model]']) == 0:
                continue
            if '[turb_base_name]' not in case:
                continue

            base_name = case['[turb_base_name]']
            # pbs_in/out dir can contain subdirs, only take the inner directory
            out_base = misc.path_split_dirs(case['[pbs_out_dir]'])[0]
            turb = case['[turb_dir]']

            self.path_pbs_e = os.path.join(out_base, turb, base_name + '.err')
            self.path_pbs_o = os.path.join(out_base, turb, base_name + '.out')
            self.path_pbs_i = os.path.join(self.pbs_in_dir, base_name + '.p')

            # apply winefix
            self.prelude = self.winefix
            # browse to scratch dir
            self.prelude += 'cd {}\n'.format(self.scratchdir)

            self.coda = '# COPY BACK FROM SCRATCH AND RENAME, remove _ at end\n'
            # copy back to turb dir at the end
            if case['[turb_db_dir]'] is not None:
                dst = os.path.join('$PBS_O_WORKDIR', case['[turb_db_dir]'],
                                   base_name)
            else:
                dst = os.path.join('$PBS_O_WORKDIR', case['[turb_dir]'],
                                   base_name)
            # FIXME: Mann64 turb exe creator adds an underscore to output
            for comp in list('uvw'):
                src = '{}_{}.bin'.format(base_name, comp)
                dst2 = '{}{}.bin'.format(dst, comp)
                self.coda += 'cp {} {}\n'.format(src, dst2)

            # alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
            rpl = (float(case['[MannAlfaEpsilon]']),
                   float(case['[MannL]']),
                   float(case['[MannGamma]']),
                   int(case['[tu_seed]']),
                   int(case['[turb_nr_u]']),
                   int(case['[turb_nr_v]']),
                   int(case['[turb_nr_w]']),
                   float(case['[turb_dx]']),
                   float(case['[turb_dy]']),
                   float(case['[turb_dz]']),
                   int(case['[high_freq_comp]']))
            params = '%1.6f %1.6f %1.6f %i %i %i %i %1.4f %1.4f %1.4f %i' % rpl
            self.execution = '%s %s %s' % (self.exe, base_name, params)
            self.create(check_dirs=True)


def eigenbody(cases, debug=False):
    """
    Read HAWC2 body eigenalysis result file
    =======================================

    This is basically a cases convience wrapper around Hawc2io.ReadEigenBody

    Parameters
    ----------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case.

    Returns
    -------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case. For each
        case, it is updated with the results, results2 of the eigenvalue
        analysis performed for each body using the following respective
        tags: [eigen_body_results] and [eigen_body_results2].

    """

    #Body data for body number : 3 with the name :nacelle
    #Results:         fd [Hz]       fn [Hz]       log.decr [%]
    #Mode nr:  1:   1.45388E-21    1.74896E-03    6.28319E+02

    for case in cases:
        # tags for the current case
        tags = cases[case]
        file_path = os.path.join(tags['[run_dir]'], tags['[eigenfreq_dir]'])
        # FIXME: do not assuem anything about the file name here, should be
        # fully defined in the tags/dataframe
        file_name = tags['[case_id]'] + '_body_eigen'
        # and load the eigenfrequency body results
        results, results2 = windIO.ReadEigenBody(file_path, file_name,
                                                  nrmodes=10)
        # add them to the htc_dict
        cases[case]['[eigen_body_results]'] = results
        cases[case]['[eigen_body_results2]'] = results2

    return cases

def eigenstructure(cases, debug=False):
    """
    Read HAWC2 structure eigenalysis result file
    ============================================

    This is basically a cases convience wrapper around
    windIO.ReadEigenStructure

    Parameters
    ----------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case.

    Returns
    -------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case. For each
        case, it is updated with the modes_arr of the eigenvalue
        analysis performed for the structure.
        The modes array (ndarray(3,n)) holds fd, fn and damping.
    """

    for case in cases:
        # tags for the current case
        tags = cases[case]
        file_path = os.path.join(tags['[run_dir]'], tags['[eigenfreq_dir]'])
        # FIXME: do not assuem anything about the file name here, should be
        # fully defined in the tags/dataframe
        file_name = tags['[case_id]'] + '_strc_eigen'
        # and load the eigenfrequency structure results
        modes = windIO.ReadEigenStructure(file_path, file_name, max_modes=500)
        # add them to the htc_dict
        cases[case]['[eigen_structure]'] = modes

    return cases


if __name__ == '__main__':
    pass