<commit_before>import numpy as np
from scipy.integrate import quad
from scipy.io import loadmat
import pylab as plt
import time

from openmdao.api import Component, Problem, Group


class GaussianWake(Component):

    def __init__(self, nTurbines, direction_id=0, options=None):
        super(GaussianWake, self).__init__()

        self.deriv_options['type'] = 'fd'
        self.deriv_options['form'] = 'central'
        self.deriv_options['step_size'] = 1.0e-12
        self.deriv_options['step_calc'] = 'relative'

        self.nTurbines = nTurbines
        self.direction_id = direction_id

        if options is None:
            self.radius_multiplier = 1.0
            self.nSamples = nSamples = 0
        else:
            # self.radius_multiplier = options['radius multiplier']
            self.nSamples = nSamples = options['nSamples']

        # unused but required for compatibility
        self.add_param('hubHeight', np.zeros(nTurbines), units='m')
        self.add_param('wakeCentersYT', np.zeros(nTurbines*nTurbines), units='m')
        self.add_param('wakeDiametersT', np.zeros(nTurbines*nTurbines), units='m')
        self.add_param('wakeOverlapTRel', np.zeros(nTurbines*nTurbines))

        # used
        self.add_param('turbineXw', val=np.zeros(nTurbines), units='m')
        self.add_param('turbineYw', val=np.zeros(nTurbines), units='m')
        self.add_param('turbineZ', val=np.zeros(nTurbines), units='m')
        self.add_param('yaw%i' % direction_id, np.zeros(nTurbines), units='deg')
        self.add_param('rotorDiameter', val=np.zeros(nTurbines)+126.4, units='m')
        self.add_param('Ct', np.zeros(nTurbines), desc='Turbine thrust coefficients')
        self.add_param('wind_speed', val=8.0, units='m/s')
        self.add_param('axialInduction', val=np.zeros(nTurbines)+1./3.)

        # params for Bastankhah with yaw
        self.add_param('model_params:ky', val=0.022, pass_by_object=True)
        self.add_param('model_params:kz', val=0.022, pass_by_object=True)
        self.add_param('model_params:alpha', val=2.32, pass_by_object=True)
        self.add_param('model_params:beta', val=0.154, pass_by_object=True)
        self.add_param('model_params:I', val=0.075, pass_by_object=True, desc='turbulence intensity')\

        self.add_output('wtVelocity%i' % direction_id, val=np.zeros(nTurbines), units='m/s')

        if nSamples > 0:
            # visualization input
            self.add_param('wsPositionXw', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='downwind position of desired measurements in wind ref. frame')
            self.add_param('wsPositionYw', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='crosswind position of desired measurements in wind ref. frame')
            self.add_param('wsPositionZ', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='position of desired measurements in wind ref. frame')

            # visualization output
            self.add_output('wsArray%i' % direction_id, np.zeros(nSamples), units='m/s', pass_by_object=True,
                            desc='wind speed at measurement locations')

    def solve_nonlinear(self, params, unknowns, resids):

        nTurbines = self.nTurbines
        direction_id = self.direction_id
        nSamples = self.nSamples

        # params for Bastankhah model with yaw
        ky = params['model_params:ky']
        kz = params['model_params:kz']
        alpha = params['model_params:alpha']
        beta = params['model_params:beta']
        I = params['model_params:I']

        # rename inputs and outputs
        turbineXw = params['turbineXw']
        turbineYw = params['turbineYw']
        turbineZ = params['turbineZ']
        yaw = params['yaw%i' % direction_id]
        rotorDiameter = params['rotorDiameter']
        Ct = params['Ct']
        axialInduction = params['axialInduction']
        wind_speed = params['wind_speed']

        # print(turbineXw, turbineYw)

        # Dw0[0] = Dw0[1]
        # m[0] = m[1]

        # rotorDiameter *= 2

        for i in range(0, nTurbines):
            if (Ct[i] > 0.96): # Glauert condition
                axialInduction[i] = 0.143 + np.sqrt(0.0203-0.6427*(0.889 - Ct[i]))
            else:
                axialInduction[i] = 0.5*(1.0-np.sqrt(1.0-Ct[i]))

        print(yaw)
        # NOTE: Bastankhah and <NAME>gel 2016 defines yaw as positive clockwise, the negative below accounts for this
        yaw *= -np.pi/180.

        if self.nSamples > 0:
            velX = params['wsPositionXw']
            velY = params['wsPositionYw']
            velZ = params['wsPositionZ']
        else:
            velX = np.zeros([0, 0])
            velY = np.zeros([0, 0])
            velZ = np.zeros([0, 0])


        velocitiesTurbines = np.tile(wind_speed, nTurbines)
        ws_array = np.tile(wind_speed, nSamples)

        for turb in range(0, nTurbines):
            x0 = rotorDiameter[turb] * (np.cos(yaw[turb]) * (1.0 + np.sqrt(1.0 - Ct[turb])) /
                                         (np.sqrt(2.0) * (alpha * I + beta * (1.0 - np.sqrt(1.0 - Ct[turb])))))
            theta_c_0 = 0.3 * yaw[turb] * (1.0 - np.sqrt(1.0 - Ct[turb] * np.cos(yaw[turb]))) / np.cos(yaw[turb])

            for loc in range(0, nSamples):  # at velX-locations
                deltax0 = velX[loc] - (turbineXw[turb] + x0)
                if deltax0 > 0.0:
                    sigmay = rotorDiameter[turb] * (ky * deltax0 / rotorDiameter[turb]
                                                    + np.cos(yaw[turb]) / np.sqrt(8.0))
                    sigmaz = rotorDiameter[turb] * (kz * deltax0 / rotorDiameter[turb]
                                                    + 1.0 / np.sqrt(8.0))
                    wake_offset = rotorDiameter[turb] * (
                        theta_c_0 * x0 / rotorDiameter[turb] +
                        (theta_c_0 / 14.7) * np.sqrt(np.cos(yaw[turb]) / (ky * kz * Ct[turb])) *
                        (2.9 + 1.3 * np.sqrt(1.0 - Ct[turb]) - Ct[turb]) *
                        np.log(
                            ((1.6 + np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              - np.sqrt(Ct[turb]))) /
                            ((1.6 - np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              + np.sqrt(Ct[turb])))
                        )
                    )

                    deltay = velY[loc] - (turbineYw[turb] + wake_offset)

                    deltav = wind_speed * (
                        (1.0 - np.sqrt(1.0 - Ct[turb] *
                                       np.cos(yaw[turb]) / (8.0 * sigmay * sigmaz /
                                                            (rotorDiameter[turb] ** 2)))) *
                        np.exp(-0.5 * ((deltay) / sigmay) ** 2) *
                        np.exp(-0.5 * ((velZ[loc] - turbineZ[turb]) / sigmaz) ** 2)
                    )

                    ws_array[loc] -= deltav

            for turbI in range(0, nTurbines):  # at turbineX-locations

                deltax0 = turbineXw[turbI] - (turbineXw[turb] + x0)

                if deltax0 > 0.0:
                    sigmay = rotorDiameter[turb] * (ky * deltax0 / rotorDiameter[turb]
                                                     + np.cos(yaw[turb]) / np.sqrt(8.0))
                    sigmaz = rotorDiameter[turb] * (kz * deltax0 / rotorDiameter[turb]
                                                     + 1.0 / np.sqrt(8.0))
                    wake_offset = rotorDiameter[turb] * (
                        theta_c_0 * x0 / rotorDiameter[turb] +
                        (theta_c_0 / 14.7) * np.sqrt(np.cos(yaw[turb]) / (ky * kz * Ct[turb])) *
                        (2.9 + 1.3 * np.sqrt(1.0 - Ct[turb]) - Ct[turb]) *
                        np.log(
                            ((1.6 + np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              - np.sqrt(Ct[turb]))) /<commit_msg>include full wake - causes undefined point - checking with Austin<commit_after>import numpy as np
from scipy.integrate import quad
from scipy.io import loadmat
import pylab as plt
import time

from openmdao.api import Component, Problem, Group


class GaussianWake(Component):

    def __init__(self, nTurbines, direction_id=0, options=None):
        super(GaussianWake, self).__init__()

        self.deriv_options['type'] = 'fd'
        self.deriv_options['form'] = 'central'
        self.deriv_options['step_size'] = 1.0e-12
        self.deriv_options['step_calc'] = 'relative'

        self.nTurbines = nTurbines
        self.direction_id = direction_id

        if options is None:
            self.radius_multiplier = 1.0
            self.nSamples = nSamples = 0
        else:
            # self.radius_multiplier = options['radius multiplier']
            self.nSamples = nSamples = options['nSamples']

        # unused but required for compatibility
        self.add_param('hubHeight', np.zeros(nTurbines), units='m')
        self.add_param('wakeCentersYT', np.zeros(nTurbines*nTurbines), units='m')
        self.add_param('wakeDiametersT', np.zeros(nTurbines*nTurbines), units='m')
        self.add_param('wakeOverlapTRel', np.zeros(nTurbines*nTurbines))

        # used
        self.add_param('turbineXw', val=np.zeros(nTurbines), units='m')
        self.add_param('turbineYw', val=np.zeros(nTurbines), units='m')
        self.add_param('turbineZ', val=np.zeros(nTurbines), units='m')
        self.add_param('yaw%i' % direction_id, np.zeros(nTurbines), units='deg')
        self.add_param('rotorDiameter', val=np.zeros(nTurbines)+126.4, units='m')
        self.add_param('Ct', np.zeros(nTurbines), desc='Turbine thrust coefficients')
        self.add_param('wind_speed', val=8.0, units='m/s')
        self.add_param('axialInduction', val=np.zeros(nTurbines)+1./3.)

        # params for Bastankhah with yaw
        self.add_param('model_params:ky', val=0.022, pass_by_object=True)
        self.add_param('model_params:kz', val=0.022, pass_by_object=True)
        self.add_param('model_params:alpha', val=2.32, pass_by_object=True)
        self.add_param('model_params:beta', val=0.154, pass_by_object=True)
        self.add_param('model_params:I', val=0.075, pass_by_object=True, desc='turbulence intensity')\

        self.add_output('wtVelocity%i' % direction_id, val=np.zeros(nTurbines), units='m/s')

        if nSamples > 0:
            # visualization input
            self.add_param('wsPositionXw', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='downwind position of desired measurements in wind ref. frame')
            self.add_param('wsPositionYw', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='crosswind position of desired measurements in wind ref. frame')
            self.add_param('wsPositionZ', np.zeros(nSamples), units='m', pass_by_object=True,
                           desc='position of desired measurements in wind ref. frame')

            # visualization output
            self.add_output('wsArray%i' % direction_id, np.zeros(nSamples), units='m/s', pass_by_object=True,
                            desc='wind speed at measurement locations')

    def solve_nonlinear(self, params, unknowns, resids):

        nTurbines = self.nTurbines
        direction_id = self.direction_id
        nSamples = self.nSamples

        # params for Bastankhah model with yaw
        ky = params['model_params:ky']
        kz = params['model_params:kz']
        alpha = params['model_params:alpha']
        beta = params['model_params:beta']
        I = params['model_params:I']

        # rename inputs and outputs
        turbineXw = params['turbineXw']
        turbineYw = params['turbineYw']
        turbineZ = params['turbineZ']
        yaw = params['yaw%i' % direction_id]
        rotorDiameter = params['rotorDiameter']
        Ct = params['Ct']
        axialInduction = params['axialInduction']
        wind_speed = params['wind_speed']

        # print(turbineXw, turbineYw)

        # Dw0[0] = Dw0[1]
        # m[0] = m[1]

        # rotorDiameter *= 2

        for i in range(0, nTurbines):
            if (Ct[i] > 0.96): # Glauert condition
                axialInduction[i] = 0.143 + np.sqrt(0.0203-0.6427*(0.889 - Ct[i]))
            else:
                axialInduction[i] = 0.5*(1.0-np.sqrt(1.0-Ct[i]))

        print(yaw)
        # NOTE: Bastankhah and <NAME> 2016 defines yaw as positive clockwise, the negative below accounts for this
        yaw *= -np.pi/180.

        if self.nSamples > 0:
            velX = params['wsPositionXw']
            velY = params['wsPositionYw']
            velZ = params['wsPositionZ']
        else:
            velX = np.zeros([0, 0])
            velY = np.zeros([0, 0])
            velZ = np.zeros([0, 0])


        velocitiesTurbines = np.tile(wind_speed, nTurbines)
        ws_array = np.tile(wind_speed, nSamples)

        for turb in range(0, nTurbines):
            x0 = rotorDiameter[turb] * (np.cos(yaw[turb]) * (1.0 + np.sqrt(1.0 - Ct[turb])) /
                                         (np.sqrt(2.0) * (alpha * I + beta * (1.0 - np.sqrt(1.0 - Ct[turb])))))
            theta_c_0 = 0.3 * yaw[turb] * (1.0 - np.sqrt(1.0 - Ct[turb] * np.cos(yaw[turb]))) / np.cos(yaw[turb])

            for loc in range(0, nSamples):  # at velX-locations
                deltax0 = velX[loc] - (turbineXw[turb] + x0)
                if deltax0 + x0 > 0.0:
                    sigmay = rotorDiameter[turb] * (ky * deltax0 / rotorDiameter[turb]
                                                    + np.cos(yaw[turb]) / np.sqrt(8.0))
                    sigmaz = rotorDiameter[turb] * (kz * deltax0 / rotorDiameter[turb]
                                                    + 1.0 / np.sqrt(8.0))
                    wake_offset = rotorDiameter[turb] * (
                        theta_c_0 * x0 / rotorDiameter[turb] +
                        (theta_c_0 / 14.7) * np.sqrt(np.cos(yaw[turb]) / (ky * kz * Ct[turb])) *
                        (2.9 + 1.3 * np.sqrt(1.0 - Ct[turb]) - Ct[turb]) *
                        np.log(
                            ((1.6 + np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              - np.sqrt(Ct[turb]))) /
                            ((1.6 - np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              + np.sqrt(Ct[turb])))
                        )
                    )

                    deltay = velY[loc] - (turbineYw[turb] + wake_offset)

                    deltav = wind_speed * (
                        (1.0 - np.sqrt(1.0 - Ct[turb] *
                                       np.cos(yaw[turb]) / (8.0 * sigmay * sigmaz /
                                                            (rotorDiameter[turb] ** 2)))) *
                        np.exp(-0.5 * ((deltay) / sigmay) ** 2) *
                        np.exp(-0.5 * ((velZ[loc] - turbineZ[turb]) / sigmaz) ** 2)
                    )

                    ws_array[loc] -= deltav

            for turbI in range(0, nTurbines):  # at turbineX-locations

                deltax0 = turbineXw[turbI] - (turbineXw[turb] + x0)

                if deltax0 + x0 > 0.0:
                    sigmay = rotorDiameter[turb] * (ky * deltax0 / rotorDiameter[turb]
                                                     + np.cos(yaw[turb]) / np.sqrt(8.0))
                    sigmaz = rotorDiameter[turb] * (kz * deltax0 / rotorDiameter[turb]
                                                     + 1.0 / np.sqrt(8.0))
                    wake_offset = rotorDiameter[turb] * (
                        theta_c_0 * x0 / rotorDiameter[turb] +
                        (theta_c_0 / 14.7) * np.sqrt(np.cos(yaw[turb]) / (ky * kz * Ct[turb])) *
                        (2.9 + 1.3 * np.sqrt(1.0 - Ct[turb]) - Ct[turb]) *
                        np.log(
                            ((1.6 + np.sqrt(Ct[turb])) *
                             (1.6 * np.sqrt(8.0 * sigmay * sigmaz /
                                            (np.cos(yaw[turb]) * rotorDiameter[turb] ** 2))
                              - np.sqrt(Ct[turb]))) /
{'fraction_non_alphanumeric': 0.13447680304925258, 'fraction_numerical': 0.022750282889643262, 'mean_word_length': 1.6481627503548337, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4870465', 'n_tokens_mistral': 5834, 'n_tokens_neox': 5395, 'n_words': 958}