hawcstab2.py

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 14 14:12:58 2014

@author: dave
"""
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import range
from io import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import object

import os
import re

import numpy as np
import pandas as pd

from wetb.prepost import (mplutils, misc)


class dummy(object):
    def __init__(self, name='dummy'):
        self.__name__ = name

regex_units = re.compile('(\\[.*?\\])')


def ReadFileHAWCStab2Header(fname):
    """
    Read a file with a weird HAWCStab2 header that starts with a #, and
    includes the column number and units between square brackets.
    """

    def _read(fname, header=0, widths=[20]*15, skipfooter=0):
        df = pd.read_fwf(fname, header=header, widths=widths,
                         skipfooter=skipfooter)
        units = regex_units.findall(''.join(df.columns))
        return df, units

    with open(fname) as f:
        line = f.readline()

    # when gradients are included in the output
    if len(line) > 800:
        df, units = _read(fname, header=1, widths=[30]*27)
        # column name has the name, unit and column number in it...
        df.columns = [k[:-2].replace('#', '').strip() for k in df.columns]
        return df, units
    elif len(line) > 200:
        df, units = _read(fname, header=0, widths=[20]*15)
        # column name has the name, unit and column number in it...
        df.columns = [k[:-2].replace('#', '').strip() for k in df.columns]
        return df, units
    # older versions of HS2 seem to have two columns less
    else:
        df, units = _read(fname, header=0, widths=[14]*13)
        df.columns = [k.replace('#', '').strip() for k in df.columns]
        return df, units


class InductionResults(object):
    """Column width can vary between versions and with/withouth gradient in
    output. Use get_col_width() for automatic detection.
    """
    def __init__(self, colwidth=14):
        """with gradients currently ind has columns width of 28 instead of 14!
        """
        self.cw = colwidth

    def get_col_width(self, fname):
        # figure out column width
        with open(fname) as fid:
            fid.readline()
            line2 = fid.readline()
        cols = misc.remove_items(line2.split(' '), '')
        if len(cols[0]) > 15:
            self.cw = 28

    def read(self, fname):
        self.data = np.loadtxt(fname)
        self.wsp = int(fname.split('_u')[-1][:-4]) / 1000.0
        try:
            self.df_data = pd.read_fwf(fname, header=0, widths=[self.cw]*38)
        except:
            self.df_data = pd.read_fwf(fname, header=0, widths=[self.cw]*34)
        # sanitize the headers
        cols = self.df_data.columns
        self.df_data.columns = [k[:-2].replace('#', '').strip() for k in cols]


class results(object):
    """
    Loading HAWCStab2 result files
    """

    def __init__(self):
        pass

    def load_pwr(self, fname):
        pwr = np.loadtxt(fname)

        res = dummy()

        res.wind = pwr[:,0]
        res.power = pwr[:,1]
        res.thrust = pwr[:,2]
        res.cp = pwr[:,3]
        res.ct = pwr[:,4]
        res.pitch_deg = pwr[:,8]
        res.rpm = pwr[:,9]

        return res

    def load_pwr_df(self, fname):
        return ReadFileHAWCStab2Header(fname)

    def load_cmb(self, fname):
        # aero-(servo)-elastic results for HS2>=2.14 have real_eig as 3th set
        with open(fname) as f:
            header = f.readline().replace('\n', '')
            cmb = np.loadtxt(f)
        # first column header has the units of the type of opearing points
        cols = misc.remove_items(header.split(']')[1].split(' '), '')
        nrmodes = np.array([int(k) for k in cols]).max()
        # when there is only data for one operating condition we only have one
        # row and consequently only a 1D array
        if len(cmb.shape) == 1:
            cmb = cmb.reshape( (1, cmb.shape[0]) )
        wind = cmb[:,0]
        freq = cmb[:,1:nrmodes+1]
        damp = cmb[:,nrmodes+1:nrmodes*2+1]
        real_eig = None
        if cmb.shape[1] > nrmodes*2+1:
            real_eig = cmb[:,nrmodes*2+1:]

        return wind, freq, damp, real_eig

    def load_cmb_df(self, fname):
        # index name can be rotor speed or wind speed
        with open(fname) as f:
            header = f.readline()
        oper_name = header.split('1')[0].strip().replace('#', '').lower()
        oper_name = oper_name.replace(' ', '').replace('[', '_')[:-1]
        oper_name = oper_name.replace('/', '')

        speed, freq, damp, real_eig = self.load_cmb(fname)
        mods = freq.shape[1]
        ops = freq.shape[0]

        df = pd.DataFrame(columns=[oper_name, 'Fd_hz', 'damp_ratio', 'mode'])
        df['Fd_hz'] = freq.flatten()
        df['damp_ratio'] = damp.flatten()
        if real_eig is not None:
            df['real_eig'] = real_eig.flatten()
        # now each mode number is a row so that means that each operating
        # point is now repeated as many times as there are modes
        df[oper_name] = speed.repeat(mods)
        modes = np.arange(1, mods+1, 1)
        df['mode'] = modes.reshape((1,mods)).repeat(ops, axis=0).flatten()

        return df

    def load_frf(self, fname, nr_inputs=3):
        frf = np.loadtxt(fname)

        self.nr_outputs = ((frf.shape[1] - 1) / 2) / nr_inputs
        self.nr_inputs = nr_inputs

        return frf

    def load_ind(self, fname, colwidth=None):
        """for results withouth gradients, colwidth=14, otherwise 28. Set to
        None to derive automatically.
        """
        self.ind = InductionResults(colwidth=colwidth)
        if colwidth is None:
            self.ind.get_col_width(fname)
        self.ind.read(fname)

    def load_amp(self, fname):

        with open(fname) as f:
            line = f.readline()

        width = 14
        nrcols = int((len(line)-1)/width)
        # first columns has one extra character
        # col nr1: rotor speed, col nr2: radius
        widths = [width+1] + [width]*(nrcols-1)
        # last line is empty
        df = pd.read_fwf(fname, header=2, widths=widths, skipfooter=1)
        units = regex_units.findall(''.join(df.columns))
        # no column number in the column name
        # since U_x, u_y, phase and theta will be repeated as many times as
        # there are modes, add the mode number in the column name
        columns = [k.replace('#', '').strip() for k in df.columns]
        nrmodes = int((len(columns) - 2 )/6)
        for k in range(nrmodes):
            for i in range(6):
                j = 2+k*6+i
                columns[j] = columns[j].split('.')[0] + ' nr%i' % (k+1)
        df.columns = columns

        return df, units

    def load_operation(self, fname):

        operation = np.loadtxt(fname, skiprows=1)
        # when the array is empty, set operation to an empty DataFrame
        if len(operation) == 0:
            cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm']
            self.operation = pd.DataFrame(columns=cols)
            return
        # when there is only one data point, the array is 1D, we allways need
        # a 2D array otherwise the columns become rows in the DataFrame
        elif len(operation.shape) == 1:
            operation = operation.reshape((1, operation.shape[0]))
        try:
            cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm']
            self.operation = pd.DataFrame(operation, columns=cols)
        except ValueError:
            cols = ['windspeed', 'pitch_deg', 'rotorspeed_rpm', 'P_aero',
                    'T_aero']
            self.operation = pd.DataFrame(operation, columns=cols)

    def load_matrices(self, fpath, basename, operating_point=1,
                      control_mat=False, local_wind_mat=False):
        """Load HAWCStab2 State Space system matrices

        The general file name format is:
        BASENAMETYPE_ase_ops_OPERATING_POINT_NUMBER.dat

        Where TYPE can be of the following:
            * amat, bmat, bvmat, cmat, dmat, dvmat, emat, fmat, fvmat

        Additionally, when including the control matrices:
            * BASENAMETYPE_ops_OPERATING_POINT_NUMBER.dat
            * TYPE: acmat, bcmat, ccmat, dcmat

        Or when including local wind speed
            * BASENAMETYPE_ase_ops_OPERATING_POINT_NUMBER.dat
            * TYPE: bvmat_loc_v, dvmat_loc_v, fvmat_loc_v


        Parameters
        ----------

        fpath : str

        basename : str

        operating_point : int, default=1


        Returns
        -------

        matrices : dict

        """
        mnames = ['amat', 'bmat', 'bvmat', 'cmat', 'dmat', 'dvmat', 'emat',
                  'fmat', 'fvmat']
        mnames_c = ['acmat', 'bcmat', 'ccmat', 'dcmat']
        mnames_v = ['bvmat_loc_v', 'dvmat_loc_v', 'fvmat_loc_v']

        if control_mat:
            mnames += mnames_c
        if local_wind_mat:
            mnames += mnames_v

        matrices = {}

        ase_template = '{:s}{:s}_ase_ops_{:d}.dat'
        ops_template = '{:s}{:s}_ops_{:d}.dat'

        for mname in mnames:
            rpl = (basename, mname, operating_point)
            template = ase_template
            if mname in mnames_c:
                template = ops_template
            fname = os.path.join(fpath, template.format(*rpl))
            matrices[mname] = np.loadtxt(fname)

        return matrices

    def write_ae_sections_h2(self):
        """
        Get the aerosection positions from the HS2 ind result file and
        write them as outputs for HAWC2
        """
        self.ind

    def plot_pwr(self, figname, fnames, labels=[], figsize=(11,7.15), dpi=120):

        results = []
        if isinstance(fnames, list):
            if len(fnames) > 4:
                raise ValueError('compare up to maximum 4 HawcStab2 cases')
            for fname in fnames:
                results.append(self.load_pwr(fname))
                # if the labels are not defined, take the file name
                if len(labels) < len(fnames):
                    labels.append(os.path.basename(fname))
        else:
            results.append(self.load_pwr(fname))

        colors = list('krbg')
        symbols = list('o<+x')
        alphas = [1.0, 0.9, 0.8, 0.75]

        fig, axes = mplutils.subplots(nrows=2, ncols=2, figsize=figsize,
                                       dpi=dpi, num=0)
        for i, res in enumerate(results):
            ax = axes[0,0]
            ax.plot(res.wind, res.power, color=colors[i],
                    label='Power %s ' % labels[i],
                    marker=symbols[i], ls='-', alpha=alphas[i])
            ax.set_title('Aerodynamic Power [kW]')#, RPM')

            ax = axes[0,1]
            ax.plot(res.wind, res.pitch_deg, color=colors[i],
                    label='Pitch %s' % labels[i],
                    marker=symbols[i], ls='-', alpha=alphas[i])
            ax.plot(res.wind, res.rpm, color=colors[i],
                    label='RPM %s ' % labels[i],
                    marker=symbols[i], ls='--', alpha=alphas[i])
            ax.set_title('Pitch [deg], RPM')

            ax = axes[1,0]
            ax.plot(res.wind, res.thrust, color=colors[i], label=labels[i],
                    marker=symbols[i], ls='-', alpha=alphas[i])
            ax.set_title('Thrust [kN]')

            ax = axes[1,1]
            ax.plot(res.wind, res.cp, label='$C_p$ %s ' % labels[i], ls='-',
                    color=colors[i], marker=symbols[i], alpha=alphas[i])
            ax.plot(res.wind, res.ct, label='$C_t$ %s ' % labels[i], ls='--',
                    color=colors[i], marker=symbols[i], alpha=alphas[i])
            ax.set_title('Power and Thrust coefficients [-]')

        for ax in axes.ravel():
            ax.legend(loc='best')
            ax.grid(True)
            ax.set_xlim([res.wind[0], res.wind[-1]])
        fig.tight_layout()

        print('saving figure: %s ... ' % figname, end='')
        figpath = os.path.dirname(figname)
        if not os.path.exists(figpath):
            os.makedirs(figpath)
        fig.savefig(figname)
        fig.clear()
        print('done!')


class ReadControlTuning(object):

    def __init__(self):
        """
        """
        pass

    def parse_line(self, line, controller):

        split1 = line.split('=')
        var1 = split1[0].strip()
        try:
            val1 = float(split1[1].split('[')[0])
            attr = getattr(self, controller)
            setattr(attr, var1, val1)

            if len(split1) > 2:
                var2 = split1[1].split(',')[1].strip()
                val2 = float(split1[2].split('[')[0])
                setattr(attr, var2, val2)
        except IndexError:
            pass

    def read_parameters(self, fpath):
        """
        Read the controller tuning file
        ===============================

        """

        with open(fpath, "r") as f:
            for i, line in enumerate(f):
                if i == 0:
                    controller = 'pi_gen_reg1'
                    setattr(self, controller, dummy())
                elif i == 2:
                    controller = 'pi_gen_reg2'
                    setattr(self, controller, dummy())
                elif i == 6:
                    controller = 'pi_pitch_reg3'
                    setattr(self, controller, dummy())
                elif i == 10:
                    controller = 'aero_damp'
                    setattr(self, controller, dummy())
                else:
                    self.parse_line(line, controller)

        # set some parameters to zero for the linear case, or when aerodynamic
        # gain scheduling is not used
        if not hasattr(self.pi_gen_reg2, 'Kd'):
            setattr(self.pi_gen_reg2, 'Kd', 0.0)
        if not hasattr(self.pi_pitch_reg3, 'Kd'):
            setattr(self.pi_pitch_reg3, 'Kd', 0.0)
        if not hasattr(self.pi_pitch_reg3, 'K2'):
            setattr(self.pi_pitch_reg3, 'K2', 0.0)
        if not hasattr(self.aero_damp, 'Kp2'):
            setattr(self.aero_damp, 'Kp2', 0.0)
        if not hasattr(self.aero_damp, 'Ko1'):
            setattr(self.aero_damp, 'Ko1', 0.0)
        if not hasattr(self.aero_damp, 'Ko2'):
            setattr(self.aero_damp, 'Ko2', 0.0)

    def parameters2tags(self):
        """Convert the tuning parameters into a dictionary whos keys are
        compatible with tag names in a HAWC2 master file.
        """

        tune_tags = {}

        tune_tags['[pi_gen_reg1.K]'] = self.pi_gen_reg1.K

        tune_tags['[pi_gen_reg2.I]'] = self.pi_gen_reg2.I
        tune_tags['[pi_gen_reg2.Kp]'] = self.pi_gen_reg2.Kp
        tune_tags['[pi_gen_reg2.Ki]'] = self.pi_gen_reg2.Ki
        tune_tags['[pi_gen_reg2.Kd]'] = self.pi_gen_reg2.Kd

        tune_tags['[pi_pitch_reg3.Kp]'] = self.pi_pitch_reg3.Kp
        tune_tags['[pi_pitch_reg3.Ki]'] = self.pi_pitch_reg3.Ki
        tune_tags['[pi_pitch_reg3.Kd]'] = self.pi_pitch_reg3.Kd
        tune_tags['[pi_pitch_reg3.K1]'] = self.pi_pitch_reg3.K1
        tune_tags['[pi_pitch_reg3.K2]'] = self.pi_pitch_reg3.K2

        tune_tags['[aero_damp.Kp2]'] = self.aero_damp.Kp2
        tune_tags['[aero_damp.Ko1]'] = self.aero_damp.Ko1
        tune_tags['[aero_damp.Ko2]'] = self.aero_damp.Ko2

        return tune_tags


if __name__ == '__main__':

    pass