h2_vs_hs2.py

# -*- coding: utf-8 -*-
"""
Created on Mon Nov  2 15:23:15 2015

@author: dave
"""

from __future__ import division
from __future__ import print_function

import os

import numpy as np
#import scipy.interpolate as interpolate
import pandas as pd
from matplotlib import pyplot as plt

import Simulations as sim
import dlcdefs
import hawcstab2 as hs2
import mplutils


class Configurations:
    # HAWC2
    eigenan = {'[eigen_analysis]':True, '[time stop]':0,
               '[t0]'            :0,     '[output]'  :False}
    control = {'[Free shaft rot]':True,  '[dll]'        :True,
               '[fixed_op]'      :False, '[fixed_shaft]':False,
               '[init_wr]'       :0.5,   '[pitch_bearing]':True}
    opt_h2 = {'[output]'         :True,  '[hs2]'        :False,
              '[hawc2_exe]'      :'hawc2-latest',
              '[Case folder]'    :'HAWC2', '[hawc2]'    :True}
    fix_op = {'[Free shaft rot]' :False, '[dll]'        :False,
              '[fixed_op]'       :True,  '[fixed_shaft]':False,
              '[init_wr]'        :0.5,   '[fixed_omega]':0.5,
              '[pitch_bearing]'  :False}
    # HAWCStab2
    opt_hs2 = {'[output]'         :False, '[hs2]'        :True,
               '[Free shaft rot]' :True,  '[dll]'        :False,
               '[fixed_op]'       :False, '[fixed_shaft]':False,
               '[init_wr]'        :0.5,   '[fixed_omega]':0.5,
               '[pitch_angle]'    :0.0,   '[hawc2_exe]'  :'hs2cmd-latest',
               '[Case folder]'    :'HAWCStab2', '[hawc2]':False,
               '[pitch_bearing]'  :True}

    # AERODYNAMIC MODELLING OPTIONS
    aero_simple = {'[aerocalc]':1, '[Induction]':0, '[tip_loss]':0,
                   '[Dyn stall]':0, '[t0]':100, '[time stop]':150}
    # when induction is on, especially the low wind speeds will need more time
    # to reach steady state in HAWC2 compared to when there is no induction.
    aero_full = {'[aerocalc]':1, '[Induction]':1, '[tip_loss]':1,
                 '[Dyn stall]':1, '[t0]':500, '[time stop]':550}

    blade_stiff_pitchC4 = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                           '[blade_damp_z]':0.01, '[blade_set]':4,
                           '[blade_subset]':1, '[blade_posx]':-0.75}

    blade_stiff_pitchC2 = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                           '[blade_damp_z]':0.01, '[blade_set]':4,
                           '[blade_subset]':1, '[blade_posx]':0.0}

    blade_stiff_pitch3C4 = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                           '[blade_damp_z]':0.01, '[blade_set]':4,
                           '[blade_subset]':1, '[blade_posx]':0.75}

    blade_flex50_tstiff_C14 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                               '[blade_damp_z]':0.03, '[blade_set]':1,
                               '[blade_subset]':15, '[blade_posx]':-0.75,
                               '[blade_nbodies]': 17,
                               '[c12]':False, '[c14]':True}
    blade_flex50_tstiff_C12 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                               '[blade_damp_z]':0.03, '[blade_set]':1,
                               '[blade_subset]':13, '[blade_posx]':0.0,
                               '[blade_nbodies]': 17,
                               '[c12]':True, '[c14]':False}
    blade_flex50_C14 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                        '[blade_damp_z]':0.03, '[blade_set]':1,
                        '[blade_subset]':17, '[blade_posx]':-0.75,
                        '[blade_nbodies]': 17,
                        '[c12]':False, '[c14]':True}
    blade_flex50_C12 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                        '[blade_damp_z]':0.03, '[blade_set]':1,
                        '[blade_subset]':16, '[blade_posx]':0.0,
                        '[blade_nbodies]': 17,
                        '[c12]':True, '[c14]':False}
    blade_flex89_tstiff_C14 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':23, '[blade_posx]':-0.75,
                            '[blade_nbodies]': 17,
                            '[c12]':False, '[c14]':True}
    blade_flex89_tstiff_C12 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':22, '[blade_posx]':0.0,
                            '[blade_nbodies]': 17,
                            '[c12]':True, '[c14]':False}
    blade_flex89_t50_C14 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':19, '[blade_posx]':-0.75,
                            '[blade_nbodies]': 17,
                            '[c12]':False, '[c14]':True}
    blade_flex89_t50_C12 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':18, '[blade_posx]':0.0,
                            '[blade_nbodies]': 17,
                            '[c12]':True, '[c14]':False}
    blade_flex89_t50_C12_allstC14 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':21, '[blade_posx]':-0.75,
                            '[blade_nbodies]': 17,
                            '[c12]':False, '[c14]':True}
    blade_flex89_t50_C12_allstC12 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                            '[blade_damp_z]':0.03, '[blade_set]':1,
                            '[blade_subset]':19, '[blade_posx]':0.0,
                            '[blade_nbodies]': 17,
                            '[c12]':True, '[c14]':False}

    blade_flex89_t50_C12_cgshC14_eaC12 = {'[blade_damp_x]':0.03,
                                          '[blade_damp_y]':0.03,
                                          '[blade_damp_z]':0.03,
                                          '[blade_set]'   :1,
                                          '[blade_subset]':24,
                                          '[blade_posx]'  :0.0,
                                          '[blade_nbodies]': 17,
                                          '[c12]':True, '[c14]':False}

    blade_flex = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                  '[blade_damp_z]':0.01, '[blade_set]':1,
                  '[blade_subset]':1, '[blade_posx]':-0.75}

    blade_flex_allac = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                        '[blade_damp_z]':0.01, '[blade_set]':1,
                        '[blade_subset]':2, '[blade_posx]':-0.75}

    blade_flex_allac_11 = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                           '[blade_damp_z]':0.01, '[blade_set]':1,
                           '[blade_subset]':7, '[blade_posx]':-0.75}

    blade_flex_allac_33 = {'[blade_damp_x]':0.02, '[blade_damp_y]':0.02,
                           '[blade_damp_z]':0.02, '[blade_set]':1,
                           '[blade_subset]':8, '[blade_posx]':-0.75}

    blade_flex_allac_50 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                           '[blade_damp_z]':0.03, '[blade_set]':1,
                           '[blade_subset]':9, '[blade_posx]':-0.75}

    blade_flex_allac_50_pitchC2 = {'[blade_damp_x]':0.03, '[blade_damp_y]':0.03,
                                   '[blade_damp_z]':0.03, '[blade_set]':1,
                                   '[blade_subset]':9, '[blade_posx]':0.0}

    # configurations for the B-series (which has quite a few changes)
    # B0001
    stiff_pc14_cgsheac14 = {'[blade_damp_x]':0.01, '[blade_damp_y]':0.01,
                            '[blade_damp_z]':0.01, '[blade_set]':1,
                            '[blade_subset]':3,    '[blade_posx]':-0.75,
                            '[blade_nbodies]': 17,
                            '[st_file]'     :'blade_flex_rect.st',
                            '[c12]':False, '[c14]':True,
                            '[ae_tolrel]': 1e-7,
                            '[ae_itmax]' : 2000,
                            '[ae_1relax]': 0.2}
    # B0002
    flex_tstiff_pc14_cgsheac14 = {'[blade_damp_x]':0.13, '[blade_damp_y]':0.13,
                                  '[blade_damp_z]':0.15, '[blade_set]':1,
                                  '[blade_subset]':5,    '[blade_posx]':-0.75,
                                  '[blade_nbodies]': 17,
                                  '[st_file]'     :'blade_flex_rect.st',
                                  '[c12]':False, '[c14]':True,
                                  '[ae_tolrel]': 1e-7,
                                  '[ae_itmax]' : 2000,
                                  '[ae_1relax]': 0.7}
    # B0003
    flex_pc14_cgsheac14 = {'[blade_damp_x]':0.13, '[blade_damp_y]':0.13,
                           '[blade_damp_z]':0.15, '[blade_set]':1,
                           '[blade_subset]':6,    '[blade_posx]':-0.75,
                           '[blade_nbodies]': 17,
                           '[st_file]'     :'blade_flex_rect.st',
                           '[c12]':False, '[c14]':True,
                           '[ae_tolrel]': 1e-7,
                           '[ae_itmax]' : 2000,
                           '[ae_1relax]': 0.7}
    # B0004
    flex_pc12_cgsheac12 = {'[blade_damp_x]':0.15, '[blade_damp_y]':0.15,
                           '[blade_damp_z]':0.17, '[blade_set]':1,
                           '[blade_subset]':7,    '[blade_posx]':0.00,
                           '[blade_nbodies]': 17,
                           '[st_file]'     :'blade_flex_rect.st',
                           '[c12]':True, '[c14]':False,
                           '[ae_tolrel]': 1e-7,
                           '[ae_itmax]' : 2000,
                           '[ae_1relax]': 0.7}
    # B0005, B0006
    flex_pc12_cgshc14_eac12 = {'[blade_damp_x]':0.15, '[blade_damp_y]':0.15,
                               '[blade_damp_z]':0.17, '[blade_set]':1,
                               '[blade_subset]':8,    '[blade_posx]':0.00,
                               '[blade_nbodies]': 17,
                               '[st_file]'     :'blade_flex_rect.st',
                               '[c12]':True, '[c14]':False,
                               '[ae_tolrel]': 1e-7,
                               '[ae_itmax]' : 2000,
                               '[ae_1relax]': 0.98}


    def __init__(self):
        pass

    def opt_tags_h2_eigenanalysis(self, basename):
        """Return opt_tags suitable for a standstill HAWC2 eigen analysis.
        """
        opt_tags = [self.opt_h2.copy()]
        opt_tags[0].update(self.eigenan.copy())
        opt_tags[0]['[Case id.]'] = '%s_hawc2_eigenanalysis' % basename
        opt_tags[0]['[blade_damp_x]'] = 0.0
        opt_tags[0]['[blade_damp_y]'] = 0.0
        opt_tags[0]['[blade_damp_z]'] = 0.0
        opt_tags[0]['[blade_nbodies]'] = 1
        opt_tags[0]['[Windspeed]'] = 0.0
        opt_tags[0]['[init_wr]'] = 0.0
        opt_tags[0]['[operational_data]'] = 'case-turbine2-empty.opt'

        return opt_tags

    def opt_tags_hs_structure_body_eigen(self, basename):
        """Return opt_tags suitable for a standstill HAWCStab2 body eigen
        analysis, at 0 RPM.
        """
        opt_tags = [self.opt_hs2.copy()]
        opt_tags[0]['[Case id.]'] = '%s_hawc2_eigenanalysis' % basename
        opt_tags[0]['[blade_damp_x]'] = 0.0
        opt_tags[0]['[blade_damp_y]'] = 0.0
        opt_tags[0]['[blade_damp_z]'] = 0.0
        opt_tags[0]['[blade_nbodies]'] = 1
        opt_tags[0]['[Windspeed]'] = 0.0
        opt_tags[0]['[init_wr]'] = 0.0
        opt_tags[0]['[fixed_omega]'] = 0.0
        opt_tags[0]['[operational_data]'] = 'case-turbine2-empty.opt'

        return opt_tags

    def opt_tags_hs2(self, basename):

        opt_tags = [self.opt_hs2.copy()]
        opt_tags[0]['[Case id.]'] = '%s_hawcstab2' % basename
        return opt_tags

    def set_hs2opdata(self, master, basename):
        """Load the HS2 operational data file and create opt_tags for HAWC2
        cases.

        Returns
        -------
        opt_tags : list of dicts
        """
        fpath = os.path.join(master.tags['[data_dir]'],
                             master.tags['[operational_data]'])
        hs2_res = hs2.results()
        hs2_res.load_operation(fpath)
        omegas = hs2_res.operation.rotorspeed_rpm.values*np.pi/30.0
        winds = hs2_res.operation.windspeed.values
        pitchs = hs2_res.operation.pitch_deg.values

        return self.set_opdata(winds, pitchs, omegas, basename=basename)

    def set_opdata(self, winds, pitchs, omegas, basename=None):
        """Return opt_tags for HAWC2 based on an HAWCStab2 operational data
        file.

        Parameters
        ----------

        winds : ndarray(n)
            wind speed for given operating point [m/s]

        pitchs : ndarray(n)
            pitch angle at given operating point [deg]

        omegas : ndarray(n)
            rotor speed at given operating point [rad/s]

        basename : str, default=None
            If not None, the [Case id.] tag is composed out of the basename,
            wind speed, pitch angle and rotor speed. If set to None, the
            [Case id.] tag is not set.

        Returns
        -------
        opt_tags : list of dicts
        """

        # the HAWC2 cases
        opt_tags = []
        for wind, pitch, omega in zip(winds, pitchs, omegas):
            opt_dict = {}
            opt_dict.update(self.opt_h2.copy())
            opt_dict.update(self.fix_op.copy())
            rpl = (basename, wind, pitch, omega)
            if basename is not None:
                tmp = '%s_%02.0fms_%04.01fdeg_%04.02frads_hawc2' % rpl
                opt_dict['[Case id.]'] = tmp
            opt_dict['[Windspeed]'] = wind
            opt_dict['[pitch_angle]'] = pitch
            opt_dict['[fixed_omega]'] = omega
            opt_dict['[init_wr]'] = omega
#            opt_dict['[t0]'] = int(2000.0/opt_dict['[Windspeed]']) # or 2000?
#            opt_dict['[time stop]'] = opt_dict['[t0]']+100
#            opt_dict['[time_stop]'] = opt_dict['[t0]']+100
            opt_tags.append(opt_dict.copy())
        return opt_tags


class Sims(object):

    def __init__(self, sim_id, P_MASTERFILE, MASTERFILE, P_SOURCE, P_RUN,
                 PROJECT, POST_DIR):
        """
        Create HtcMaster() object
        =========================

        the HtcMaster contains all the settings to start creating htc files.
        It holds the master file, server paths and more.

        The master.tags dictionary holds those tags who do not vary for different
        cases. Variable tags, i.e. tags who are a function of other variables
        or other tags, are defined in the function variable_tag_func().

        It is considered as good practice to define the default values for all
        the variable tags in the master_tags

        Members
        -------

        Returns
        -------

        """

        self.sim_id = sim_id
        self.P_MASTERFILE = P_MASTERFILE
        self.MASTERFILE = MASTERFILE
        self.P_SOURCE = P_SOURCE
        self.P_RUN = P_RUN
        self.PROJECT = PROJECT
        self.POST_DIR = POST_DIR

        # TODO: write a lot of logical tests for the tags!!
        # TODO: tests to check if the dirs are setup properly (ending slahses)
        # FIXME: some tags are still variable! Only static tags here that do
        # not depent on any other variable that can change
        self.master = sim.HtcMaster()
        self.set_tag_defaults()

    def _var_tag_func(self, master, case_id_short=False):
        """
        Function which updates HtcMaster.tags and returns an HtcMaster object

        Only use lower case characters for case_id since a hawc2 result and
        logfile are always in lower case characters. Simulations.prepare_launch
        will force the value of the tags as defined in master.output_dirs
        to lower case.

        BE CAREFULL: if you change a master tag that is used to dynamically
        calculate an other tag, that change will be propageted over all cases,
        for example:
        master.tags['tag1'] *= master.tags[tag2]*master.tags[tag3']
        it will accumlate over each new case. After 20 cases
        master.tags['tag1'] = (master.tags[tag2]*master.tags[tag3'])^20
        which is not wanted, you should do
        master.tags['tag1'] = tag1_base*master.tags[tag2]*master.tags[tag3']

        """

        mt = master.tags

        dlc_case = mt['[Case folder]']
        mt['[data_dir]'] = 'data/'
        mt['[res_dir]'] = 'res/%s/' % dlc_case
        mt['[log_dir]'] = 'logfiles/%s/' % dlc_case
        mt['[htc_dir]'] = 'htc/%s/' % dlc_case
        mt['[case_id]'] = mt['[Case id.]']
        mt['[DLC]'] = dlc_case
        mt['[pbs_out_dir]'] = 'pbs_out/%s/' % dlc_case
        mt['[pbs_in_dir]'] = 'pbs_in/%s/' % dlc_case
        mt['[iter_dir]'] = 'iter/%s/' % dlc_case

        if mt['[eigen_analysis]']:
            rpl = os.path.join(dlc_case, mt['[Case id.]'])
            mt['[eigenfreq_dir]'] = 'res_eigen/%s/' % rpl

        # for HAWCStab2 certain things have to be done differently
        if mt['[hs2]']:
            mt['[htc_dir]'] = ''
            mt['[t0]'] = 0
            mt['[time stop]'] = 1
            mt['[hawc2]'] = False
            mt['[output]'] = False
            mt['[copyback_files]'] = ['./*.ind', './*.pwr', './*.log',
                                      './*.cmb', './*.bea']
            mt['[copyback_frename]'] = [mt['[res_dir]'], mt['[res_dir]'],
                                        mt['[log_dir]'], mt['[res_dir]'],
                                        mt['[res_dir]']]
            if mt['[hs2_bladedeform_switch]']:
                mt['[hs2_bladedeform]'] = 'bladedeform'
            else:
                mt['[hs2_bladedeform]'] = 'nobladedeform'

            if int(mt['[tip_loss]']) == 1:
                mt['[hs2_tipcorrect]'] = 'tipcorrect'
            else:
                mt['[hs2_tipcorrect]'] = 'notipcorrect'

            if int(mt['[Induction]']) == 1:
                mt['[hs2_induction]'] = 'induction'
            else:
                mt['[hs2_induction]'] = 'noinduction'

            if mt['[hs2_gradients_switch]']:
                mt['[hs2_gradients]'] = 'gradients'
            else:
                mt['[hs2_gradients]'] = 'nogradients'

        mt['[windspeed]'] = mt['[Windspeed]']
        mt['[time_stop]'] = mt['[time stop]']
        mt['[duration]'] = str(float(mt['[time_stop]']) - float(mt['[t0]']))

        return master

    def _set_path_auto_config(self, verbose=True):
        """
        auto configure directories: assume you are running in the root of the
        relevant HAWC2 model
        and assume we are in a simulation case of a certain turbine/project
        """
        (self.P_RUN, self.P_SOURCE, self.PROJECT,
             self.sim_id, self.P_MASTERFILE,
             self.MASTERFILE, self.POST_DIR) = dlcdefs.configure_dirs(verbose=verbose)

    def _set_path_config(self, runmethod='here'):
        """
        Set the path configuration into the tags
        """

        self.runmethod = runmethod

        if runmethod == 'here':
            self._set_path_auto_config()
        elif runmethod in ['local', 'local-script', 'none', 'local-ram']:
            self.p_root = '/home/dave/SimResults/h2_vs_hs2/'
        elif runmethod == 'windows-script':
            self.p_root = '/mnt/D16731/dave/Documents/_SimResults'
        elif runmethod == 'gorm':
            self.p_root = '/mnt/hawc2sim/h2_vs_hs2'
        elif runmethod == 'jess':
            self.p_root = '/mnt/hawc2sim/h2_vs_hs2'
        else:
            msg='unsupported runmethod, options: none, local, gorm or opt'
            raise ValueError, msg

        if not runmethod == 'here':
            self.P_RUN = os.path.join(self.p_root, self.PROJECT, self.sim_id)

        self.master.tags['[master_htc_file]'] = self.MASTERFILE
        self.master.tags['[master_htc_dir]'] = self.P_MASTERFILE
        # directory to data, htc, SOURCE DIR
        if self.P_SOURCE[-1] == os.sep:
            self.master.tags['[model_dir_local]']  = self.P_SOURCE
        else:
            self.master.tags['[model_dir_local]']  = self.P_SOURCE + os.sep
        if self.P_RUN[-1] == os.sep:
            self.master.tags['[run_dir]'] = self.P_RUN
        else:
            self.master.tags['[run_dir]'] = self.P_RUN + os.sep

        self.master.tags['[post_dir]'] = self.POST_DIR
        self.master.tags['[sim_id]'] = self.sim_id
        # set the model_zip tag to include the sim_id
        rpl = (self.PROJECT, self.master.tags['[sim_id]'])
        self.master.tags['[model_zip]'] = '%s_%s.zip' % rpl

    def set_tag_defaults(self):
        """
        Set the default values of the required master tags
        """
        mt = self.master.tags

        # other required tags and their defaults
        mt['[dt_sim]'] = 0.01
        mt['[hawc2_exe]'] = 'hawc2-latest'
        # convergence_limits  0.001  0.005  0.005 ;
        # critical one, risidual on the forces: 0.0001 = 1e-4
        mt['[epsresq]'] = '1.0' # default=10.0
        # increment residual
        mt['[epsresd]'] = '0.5' # default= 1.0
        # constraint equation residual
        mt['[epsresg]'] = '1e-8' # default= 1e-7
        # folder names for the saved results, htc, data, zip files
        # Following dirs are relative to the model_dir_server and they specify
        # the location of where the results, logfiles, animation files that where
        # run on the server should be copied to after the simulation has finished.
        # on the node, it will try to copy the turbulence files from these dirs
        mt['[animation_dir]'] = 'animation/'
        mt['[control_dir]']   = 'control/'
        mt['[data_dir]']      = 'data/'
        mt['[eigen_analysis]'] = False
        mt['[eigenfreq_dir]'] = False
        mt['[htc_dir]']       = 'htc/'
        mt['[log_dir]']       = 'logfiles/'
        mt['[meander_dir]']   = False
        mt['[opt_dir]']       = False
        mt['[pbs_out_dir]']   = 'pbs_out/'
        mt['[res_dir]']       = 'res/'
        mt['[iter_dir]']      = 'iter/'
        mt['[turb_dir]']      = 'turb/'
        mt['[turb_db_dir]']   = '../turb/'
        mt['[wake_dir]']      = False
        mt['[hydro_dir]']     = False
        mt['[mooring_dir]']   = False
        mt['[externalforce]'] = False
        mt['[Case folder]']   = 'NoCaseFolder'
        # zip_root_files only is used when copy to run_dir and zip creation, define
        # in the HtcMaster object
        mt['[zip_root_files]'] = []
        # only active on PBS level, so files have to be present in the run_dir
        mt['[copyback_files]'] = []   # copyback_resultfile
        mt['[copyback_frename]'] = [] # copyback_resultrename
        mt['[copyto_files]'] = []     # copyto_inputfile
        mt['[copyto_generic]'] = []   # copyto_input_required_defaultname

        # In master file tags within the HAWC2 vs HAWCStab2 context
        mt['[hawc2]'] = False
        mt['[output]'] = False
        mt['[eigen_analysis]'] = False
        mt['[system_eigen_analysis]'] = False
        mt['[operational_data]'] = 'case_name.opt'

        mt['[gravity]'] = 0.0
        mt['[shaft_tilt]'] = 0.0 # 5.0
        mt['[coning]'] = 0.0 # 2.5
        mt['[Windspeed]'] = 1.0
        mt['[wtilt]'] = 0.0
        mt['[wdir]'] = 0.0
        mt['[aerocalc]'] = 1
        mt['[Induction]'] = 0
        mt['[tip_loss]'] = 0
        mt['[Dyn stall]'] = 0
        mt['[tu_model]'] = 0
        mt['[shear_exp]'] = 0
        mt['[tower_shadow]'] = 0
        mt['[TI]'] = 1
        mt['[fixed_omega]'] = 1.0
        mt['[init_wr]'] = 0
        mt['[pc_file_name]'] = 'hawc_pc.mhh'
        mt['[ae_file_name]'] = 'hawc2_ae.mhh'
        mt['[nr_ae_sections]'] = 30
        mt['[use_nr_ae_sections]'] = True
        mt['[use_ae_distrb_file]'] = False
        mt['[ae_set_nr]'] = 1
        # tors_e output depends on the pitch axis configuration
        mt['[c12]'] = False
        mt['[c14]'] = False

        mt['[t0]'] = 500
        mt['[time stop]'] = 600

        mt['[hs2]'] = False
        mt['[nr_blade_modes_hs2]'] = 10
        mt['[stab_analysis]'] = False
        mt['[steady_states]'] = True
        mt['[hs2_bladedeform_switch]'] = True
        mt['[hs2_gradients_switch]'] = False
        # by default take the stiff set
        mt['[st_file]'] = 'hawc2_st.mhh'
        mt['[tower_set]'] = 4 # 1
        mt['[shaft_set]'] = 4 # 2
        mt['[blade_set]'] = 4 # 3
        mt['[tower_subset]'] = 1
        mt['[shaft_subset]'] = 1
        mt['[blade_subset]'] = 1
        mt['[blade_nbodies]'] = 1
        mt['[blade_posx]'] = -0.75
        mt['[blade_damp_x]'] = 0.01
        mt['[blade_damp_y]'] = 0.01
        mt['[blade_damp_z]'] = 0.01
        # HAWCStab2 convergence criteria
        mt['[bem_tol]'] = 1e-12
        mt['[bem_itmax]'] = 10000
        mt['[bem_1relax]'] = 0.02
        mt['[ae_tolrel]'] = 1e-7
        mt['[ae_itmax]'] = 2000
        mt['[ae_1relax]'] = 0.5
        mt['[tol_7]'] = 10
        mt['[tol_8]'] = 5
        mt['[tol_9]'] = 1e-8

        # =========================================================================
        # basic required tags by HtcMaster and PBS in order to function properly
        # =========================================================================
        # the express queue ('#PBS -q xpresq') has a maximum walltime of 1h
        mt['[pbs_queue_command]'] = '#PBS -q workq'
        # walltime should have following format: hh:mm:ss
        mt['[walltime]'] = '04:00:00'
        mt['[auto_walltime]'] = False

    def get_dlc_casedefs(self):
        """
        Create iter_dict and opt_tags based on spreadsheets
        """

        iter_dict = dict()
        iter_dict['[empty]'] = [False]

        # see if a htc/DLCs dir exists
        dlcs_dir = os.path.join(self.P_SOURCE, 'htc', 'DLCs')
        if os.path.exists(dlcs_dir):
            opt_tags = dlcdefs.excel_stabcon(dlcs_dir)
        else:
            opt_tags = dlcdefs.excel_stabcon(os.path.join(self.P_SOURCE, 'htc'))

        if len(opt_tags) < 1:
            raise ValueError('There are is not a single case defined. Make sure '
                             'the DLC spreadsheets are configured properly.')

        # add all the root files, except anything with *.zip
        f_ziproot = []
        for (dirpath, dirnames, fnames) in os.walk(self.P_SOURCE):
            # remove all zip files
            for i, fname in enumerate(fnames):
                if fname.endswith('.zip'):
                    fnames.pop(i)
            f_ziproot.extend(fnames)
            break
        # and add those files
        for opt in opt_tags:
            opt['[zip_root_files]'] = f_ziproot

        self.master.output_dirs.extend('[Case folder]')
        self.master.output_dirs.extend('[Case id.]')

        return iter_dict, opt_tags

    def create_inputs(self, iter_dict, opt_tags):

        sim.prepare_launch(iter_dict, opt_tags, self.master, self._var_tag_func,
                           write_htc=True, runmethod=self.runmethod, verbose=False,
                           copyback_turb=False, msg='', update_cases=False,
                           ignore_non_unique=False, run_only_new=False,
                           pbs_fname_appendix=False, short_job_names=False)

    def get_control_tuning(self, fpath):
        """
        Read a HAWCStab2 controller tuning file and return as tags
        """
        tuning = hs2.hs2_control_tuning()
        tuning.read_parameters(fpath)

        tune_tags = {}

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

        tune_tags['[pi_gen_reg2.I]'] = tuning.pi_gen_reg2.I
        tune_tags['[pi_gen_reg2.Kp]'] = tuning.pi_gen_reg2.Kp
        tune_tags['[pi_gen_reg2.Ki]'] = tuning.pi_gen_reg2.Ki

        tune_tags['[pi_pitch_reg3.Kp]'] = tuning.pi_pitch_reg3.Kp
        tune_tags['[pi_pitch_reg3.Ki]'] = tuning.pi_pitch_reg3.Ki
        tune_tags['[pi_pitch_reg3.K1]'] = tuning.pi_pitch_reg3.K1
        tune_tags['[pi_pitch_reg3.K2]'] = tuning.pi_pitch_reg3.K2

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

        return tune_tags

    def post_processing(self, statistics=True, resdir=None):
        """
        Parameters
        ----------

        resdir : str, default=None
            Defaults to reading the results from the [run_dir] tag.
            Force to any other directory using this variable. You can also use
            the presets as defined for runmethod in _set_path_config.
        """

        post_dir = self.POST_DIR

        # =========================================================================
        # check logfiles, results files, pbs output files
        # logfile analysis is written to a csv file in logfiles directory
        # =========================================================================
        # load the file saved in post_dir
        cc = sim.Cases(post_dir, self.sim_id, rem_failed=False)

        if resdir is None:
            # we keep the run_dir as defined during launch
            run_root = None
        elif resdir in ['local', 'local-script', 'none', 'local-ram']:
            run_root = '/home/dave/SimResults'
        elif resdir == 'windows-script':
            run_root = '/mnt/D16731/dave/Documents/_SimResults'
        elif resdir == 'gorm':
            run_root = '/mnt/hawc2sim/h2_vs_hs2'
        elif resdir == 'jess':
            run_root = '/mnt/hawc2sim/h2_vs_hs2'
        else:
            run_root = None
            cc.change_results_dir(resdir)

        if isinstance(run_root, str):
            forcedir = os.path.join(run_root, self.PROJECT, self.sim_id)
            cc.change_results_dir(forcedir)

        cc.post_launch()
        cc.remove_failed()

        if statistics:
            tags=['[windspeed]']
            stats_df = cc.statistics(calc_mech_power=False, ch_fatigue=[],
                                     tags=tags, update=False)
            ftarget = os.path.join(self.POST_DIR, '%s_statistics.xlsx')
            stats_df.to_excel(ftarget % self.sim_id)


class MappingsH2HS2(object):

    def __init__(self, chord_length=3.0):
        """
        """
        self.hs2_res = hs2.results()
        self.chord_length = chord_length

    def powercurve(self, h2_df_stats, fname_hs):

        self._powercurve_h2(h2_df_stats)
        self._powercurve_hs2(fname_hs)

    def _powercurve_h2(self, df_stats):

        mappings = {'Ae rot. power' : 'P_aero',
                    'Ae rot. thrust': 'T_aero',
                    'Vrel-1-39.03'  : 'vrel_39',
                    'Omega'         : 'rotorspeed',
                    'tower-tower-node-010-forcevec-y' : 'T_towertop',
                    'tower-shaft-node-003-forcevec-y' : 'T_shafttip'}

        df_stats.sort_values('[windspeed]', inplace=True)
        df_mean = pd.DataFrame()
        df_std = pd.DataFrame()

        for key, value in mappings.iteritems():
            tmp = df_stats[df_stats['channel']==key]
            df_mean[value] = tmp['mean'].values.copy()
            df_std[value] = tmp['std'].values.copy()

        # also add the wind speed
        df_mean['windspeed'] = tmp['[windspeed]'].values.copy()
        df_std['windspeed'] = tmp['[windspeed]'].values.copy()

        self.pwr_h2_mean = df_mean
        self.pwr_h2_std = df_std
        self.h2_df_stats = df_stats

    def _powercurve_hs2(self, fname):

        mappings = {u'P [kW]'  :'P_aero',
                    u'T [kN]'  :'T_aero',
                    u'V [m/s]' :'windspeed'}

        df_pwr, units = self.hs2_res.load_pwr_df(fname)

        self.pwr_hs = pd.DataFrame()
        for key, value in mappings.iteritems():
            self.pwr_hs[value] = df_pwr[key].values.copy()

    def blade_distribution(self, fname_h2, fname_hs2, h2_df_stats=None,
                           fname_h2_tors=None):

        self.hs2_res.load_ind(fname_hs2)
        self.h2_res = sim.windIO.ReadOutputAtTime(fname_h2)
        self._distribution_hs2()
        self._distribution_h2()
        if h2_df_stats is not None:
            self.h2_df_stats = h2_df_stats
            if fname_h2_tors is not None:
                self.distribution_torsion_h2(fname_h2_tors)

    def _distribution_hs2(self):
        """Read a HAWCStab2 *.ind file (blade distribution loading)
        """

        mapping_hs2 =  {u's [m]'       :'curved_s',
                        u'CL0 [-]'     :'Cl',
                        u'CD0 [-]'     :'Cd',
                        u'CT [-]'      :'Ct',
                        u'CP [-]'      :'Cp',
                        u'A [-]'       :'ax_ind',
                        u'AP [-]'      :'tan_ind',
                        u'U0 [m/s]'    :'vrel',
                        u'PHI0 [rad]'  :'inflow_angle',
                        u'ALPHA0 [rad]':'AoA',
                        u'X_AC0 [m]'   :'pos_x',
                        u'Y_AC0 [m]'   :'pos_y',
                        u'Z_AC0 [m]'   :'pos_z',
                        u'UX0 [m]'     :'def_x',
                        u'UY0 [m]'     :'def_y',
                        u'Tors. [rad]' :'torsion',
                        u'Twist[rad]'  :'twist',
                        u'V_a [m/s]'   :'ax_ind_vel',
                        u'V_t [m/s]'   :'tan_ind_vel',
                        u'FX0 [N/m]'   :'F_x',
                        u'FY0 [N/m]'   :'F_y',
                        u'M0 [Nm/m]'   :'M'}

        try:
            hs2_cols = [k for k in mapping_hs2]
            # select only the HS channels that will be used for the mapping
            std_cols = [mapping_hs2[k] for k in hs2_cols]
            self.hs_aero = self.hs2_res.ind.df_data[hs2_cols].copy()
        except KeyError:
            # some results have been created with older HAWCStab2 that did not
            # include CT and CP columns
            mapping_hs2.pop(u'CT [-]')
            mapping_hs2.pop(u'CP [-]')
            hs2_cols = [k for k in mapping_hs2]
            std_cols = [mapping_hs2[k] for k in hs2_cols]
            # select only the HS channels that will be used for the mapping
            self.hs_aero = self.hs2_res.ind.df_data[hs2_cols].copy()

        # change column names to the standard form that is shared with H2
        self.hs_aero.columns = std_cols
        self.hs_aero['AoA'] *= (180.0/np.pi)
        self.hs_aero['inflow_angle'] *= (180.0/np.pi)
        self.hs_aero['torsion'] *= (180.0/np.pi)
#        self.hs_aero['pos_x'] = (-1.0) # self.chord_length / 4.0

    def _distribution_h2(self):
        mapping_h2 =  { u'Radius_s'  :'curved_s',
                        u'Cl'        :'Cl',
                        u'Cd'        :'Cd',
                        u'Ct_local'  :'Ct',
                        u'Cq_local'  :'Cq',
                        u'Induc_RPy' :'ax_ind_vel',
                        u'Induc_RPx' :'tan_ind_vel',
                        u'Vrel'      :'vrel',
                        u'Inflow_ang':'inflow_angle',
                        u'alfa'      :'AoA',
                        u'pos_RP_x'  :'pos_x',
                        u'pos_RP_y'  :'pos_y',
                        u'pos_RP_z'  :'pos_z',
                        u'Secfrc_RPx':'F_x',
                        u'Secfrc_RPy':'F_y',
                        u'Secmom_RPz':'M'}
        h2_cols = [k for k in mapping_h2]
        std_cols = [mapping_h2[k] for k in h2_cols]

        # select only the h2 channels that will be used for the mapping
        h2_aero = self.h2_res[h2_cols].copy()
        # change column names to the standard form that is shared with HS
        h2_aero.columns = std_cols
        h2_aero['def_x'] = self.h2_res['Pos_B_x'] - self.h2_res['Inipos_x_x']
        h2_aero['def_y'] = self.h2_res['Pos_B_y'] - self.h2_res['Inipos_y_y']
        h2_aero['def_z'] = self.h2_res['Pos_B_z'] - self.h2_res['Inipos_z_z']
        h2_aero['ax_ind_vel'] *= (-1.0)
        h2_aero['pos_x'] += (self.chord_length / 2.0)
        h2_aero['F_x'] *= (1e3)
        h2_aero['F_y'] *= (1e3)
        h2_aero['M'] *= (1e3)
#        # HAWC2 includes root and tip nodes, while HAWC2 doesn't. Remove them
#        h2_aero = h2_aero[1:-1]
        self.h2_aero = h2_aero

    def distribution_torsion_h2(self, fname_h2):
        """Determine torsion distribution from the HAWC2 result statistics.
        tors_e is in degrees.
        """
        if not hasattr(self, 'h2_aero'):
            raise UserWarning('first run blade_distribution')

        # load the HAWC2 .sel file for the channels
        fpath = os.path.dirname(fname_h2)
        fname = os.path.basename(fname_h2)
        res = sim.windIO.LoadResults(fpath, fname, readdata=False)
        sel = res.ch_df[res.ch_df.sensortype == 'Tors_e'].copy()
        sel.sort_values(['radius'], inplace=True)
        self.h2_aero['Radius_s_tors'] = sel.radius.values.copy()
        self.h2_aero['tors_e'] = sel.radius.values.copy()
        tors_e_channels = sel.ch_name.tolist()

        # find the current case in the statistics DataFrame
        case = fname.replace('.htc', '')
        df_case = self.h2_df_stats[self.h2_df_stats['[case_id]']==case].copy()
        # and select all the torsion channels
        df_tors_e = df_case[df_case.channel.isin(tors_e_channels)].copy()
        # join the stats with the channel descriptions DataFrames, have the
        # same name on the joining column
        df_tors_e.set_index('channel', inplace=True)
        sel.set_index('ch_name', inplace=True)

        # joining happens on the index, and for which the same channel has been
        # used: the unique HAWC2 channel naming scheme
        df_tors_e = pd.concat([df_tors_e, sel], axis=1)
        df_tors_e.radius = df_tors_e.radius.astype(np.float64)
        # sorting on radius, combine with ch_df
        df_tors_e.sort_values(['radius'], inplace=True)

        # FIXME: what if number of torsion outputs is less than aero
        # calculation points?
#        df_tmp = pd.DataFrame()
        self.h2_aero['torsion'] = df_tors_e['mean'].values.copy()
        self.h2_aero['torsion_std'] = df_tors_e['std'].values.copy()
        self.h2_aero['torsion_radius_s'] = df_tors_e['radius'].values.copy()
#        df_tmp = pd.DataFrame()
#        df_tmp['torsion'] = df_tors_e['mean'].copy()
#        df_tmp['torsion_std'] = df_tors_e['std'].copy()
#        df_tmp['torsion_radius_s'] = df_tors_e['radius'].copy()
#        df_tmp.set_index('')

    def body_structure_modes(self, fname_h2, fname_hs):
        self._body_structure_modes_h2(fname_h2)
        self._body_structure_modes_hs(fname_hs)

    def _body_structure_modes_h2(self, fname):
        self.body_freq_h2 = sim.windIO.ReadEigenBody(fname)

        blade_h2 = self.body_freq_h2[self.body_freq_h2['body']=='blade1'].copy()
        # because HAWCStab2 is sorted by frequency
        blade_h2.sort_values('Fd_hz', inplace=True)
        # HAWC2 usually has a lot of duplicate entries
        blade_h2.drop_duplicates('Fd_hz', keep='first', inplace=True)
        # also drop the ones with very high damping, and 0 frequency
        query = '(log_decr_pct < 500 and log_decr_pct > -500) and Fd_hz > 0.0'
        self.blade_body_freq_h2 = blade_h2.query(query)

    def _body_structure_modes_hs(self, fname):
        self.body_freq_hs = hs2.results().load_cmb_df(fname)


class Plots(object):
    """
    Comparison plots between HACW2 and HAWCStab2. This is done based on
    the HAWC2 output output_at_time, and HAWCStab2 output *.ind
    """

    def __init__(self):

        self.h2c = 'b'
        self.h2ms = '+'
        self.h2ls = '-'
        self.hsc = 'r'
        self.hsms = 'x'
        self.hsls = '--'
        self.errls = '-'
        self.errc = 'k'
        self.errlab = 'diff [\\%]'
        self.interactive = False

        self.dist_size = (16, 11)
        self.dist_channels = ['pos_x', 'pos_y', 'AoA', 'inflow_angle',
                              'Cl', 'Cd', 'vrel', 'ax_ind_vel',
                              'F_x', 'F_y', 'M']

    def load_h2(self, fname_h2, h2_df_stats=None, fname_h2_tors=None):

        res = MappingsH2HS2()
        res.h2_res = sim.windIO.ReadOutputAtTime(fname_h2)
        res._distribution_h2()
        if h2_df_stats is not None:
            res.h2_df_stats = h2_df_stats
            if fname_h2_tors is not None:
                res.distribution_torsion_h2(fname_h2_tors)

        return res

    def load_hs(self, fname_hs):

        res = MappingsH2HS2()
        res.hs2_res.load_ind(fname_hs)
        res._distribution_hs2()

        return res

    def new_fig(self, title=None, nrows=2, ncols=1, dpi=150, size=(12.0, 5.0)):

        if self.interactive:
            subplots = plt.subplots
        else:
            subplots = mplutils.subplots

        fig, axes = subplots(nrows=nrows, ncols=ncols, dpi=dpi, figsize=size)
        axes = axes.ravel()
        if title is not None:
            fig.suptitle(title)
        return fig, axes

    def set_axes_label_grid(self, axes, setlegend=False):
        for ax in axes.ravel():
            if setlegend:
                leg = ax.legend(loc='best')
                if leg is not None:
                    leg.get_frame().set_alpha(0.5)