In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload
import numpy as np
import scipy as sp
import pandas as pd
from scipy import interpolate
import matplotlib.pyplot as plt

import os,sys,inspect
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parentdir = os.path.dirname(os.path.dirname(currentdir))
#sys.path.insert(0,parentdir) 
#sys.path.append(parentdir+'/gclarsen/src/gclarsen')

import fusedwake.WindTurbine as wt
import fusedwake.WindFarm as wf
#import fusedwake.fused as fused_gcl
from fusedwake.gcl import *

In [2]:
#from gclarsen.fusedwasp import  PlantFromWWH
#wwh = PlantFromWWH(filename = parentdir+'/wind-farm-wake-model/gclarsen/src/gclarsen/test/wind_farms/horns_rev/hornsrev1_turbine_nodescription.wwh')
#wwh = PlantFromWWH(filename = 'hornsrev1_turbine_nodescription.wwh')

Verification of the FUSED-Wind wrapper

common inputs


In [9]:
wf.WindFarm?

In [11]:
v80 = wt.WindTurbine('Vestas v80 2MW offshore','V80_2MW_offshore.dat',70,40)
HR1 = wf.WindFarm(name='Horns Rev 1',yml='hornsrev.yml')#,v80)
WD = range(0,360,1)

FUSED-Wind implementation


In [12]:
"""
##Fused inputs
inputs = dict(
    wind_speed=8.0,
    roughness=0.0001,
    TI=0.05,
    NG=4,
    sup='lin',
    wt_layout = fused_gcl.generate_GenericWindFarmTurbineLayout(HR1))

fgcl = fused_gcl.FGCLarsen()
# Setting the inputs
for k,v in inputs.iteritems():
    setattr(fgcl, k, v)

fP_WF = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    fgcl.wind_direction = wd
    fgcl.run()
    fP_WF[iwd] = fgcl.power
"""


Out[12]:
"\n##Fused inputs\ninputs = dict(\n    wind_speed=8.0,\n    roughness=0.0001,\n    TI=0.05,\n    NG=4,\n    sup='lin',\n    wt_layout = fused_gcl.generate_GenericWindFarmTurbineLayout(HR1))\n\nfgcl = fused_gcl.FGCLarsen()\n# Setting the inputs\nfor k,v in inputs.iteritems():\n    setattr(fgcl, k, v)\n\nfP_WF = np.zeros([len(WD)])\nfor iwd, wd in enumerate(WD):\n    fgcl.wind_direction = wd\n    fgcl.run()\n    fP_WF[iwd] = fgcl.power\n"

pure python implementation


In [13]:
P_WF = np.zeros([len(WD)])
P_WF_v0 = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    P_WT,U_WT,CT_WT = GCLarsen(WS=8.0,z0=0.0001,TI=0.05,WD=wd,WF=HR1,NG=4,sup='lin')
    P_WF[iwd] = P_WT.sum()
    P_WT,U_WT,CT_WT = GCLarsen_v0(WS=8.0,z0=0.0001,TI=0.05,WD=wd,WF=HR1,NG=4,sup='lin')
    P_WF_v0[iwd] = P_WT.sum()


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-4a1f1cb162b9> in <module>()
      2 P_WF_v0 = np.zeros([len(WD)])
      3 for iwd, wd in enumerate(WD):
----> 4     P_WT,U_WT,CT_WT = GCLarsen(WS=8.0,z0=0.0001,TI=0.05,WD=wd,WF=HR1,NG=4,sup='lin')
      5     P_WF[iwd] = P_WT.sum()
      6     P_WT,U_WT,CT_WT = GCLarsen_v0(WS=8.0,z0=0.0001,TI=0.05,WD=wd,WF=HR1,NG=4,sup='lin')

NameError: name 'GCLarsen' is not defined

In [6]:
fig, ax = plt.subplots()
ax.plot(WD,P_WF/(HR1.WT[0].get_P(8.0)*HR1.nWT),'-o', label='python')
ax.plot(WD,P_WF_v0/(HR1.WT[0].get_P(8.0)*HR1.nWT),'-d', label='python v0')
ax.set_xlabel('wd [deg]')
ax.set_ylabel('Wind farm efficiency [-]')
ax.set_title(HR1.name)
ax.legend(loc=3)
plt.savefig(HR1.name+'_Power_wd_360.pdf')


Asserting new implementation


In [7]:
WD = 261.05
P_WT,U_WT,CT_WT = GCLarsen_v0(WS=10.,z0=0.0001,TI=0.1,WD=WD,WF=HR1, NG=5, sup='quad')
P_WT_2,U_WT_2,CT_WT_2 = GCLarsen(WS=10.,z0=0.0001,TI=0.1,WD=WD,WF=HR1, NG=5, sup='quad')

In [8]:
np.testing.assert_array_almost_equal(U_WT,U_WT_2)
np.testing.assert_array_almost_equal(P_WT,P_WT_2)


---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-8-987fd8803fdf> in <module>()
----> 1 np.testing.assert_array_almost_equal(U_WT,U_WT_2)
      2 np.testing.assert_array_almost_equal(P_WT,P_WT_2)

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/numpy/testing/utils.pyc in assert_array_almost_equal(x, y, decimal, err_msg, verbose)
    890     assert_array_compare(compare, x, y, err_msg=err_msg, verbose=verbose,
    891              header=('Arrays are not almost equal to %d decimals' % decimal),
--> 892              precision=decimal)
    893 
    894 

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/numpy/testing/utils.pyc in assert_array_compare(comparison, x, y, err_msg, verbose, header, precision)
    711                                 names=('x', 'y'), precision=precision)
    712             if not cond:
--> 713                 raise AssertionError(msg)
    714     except ValueError:
    715         import traceback

AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 90.0%)
 x: array([ 9.962871,  9.962871,  9.962871,  9.962871,  9.962871,  9.962871,
        9.962871,  9.962871,  9.641319,  9.641319,  9.641319,  9.641319,
        9.641319,  9.642722,  9.642722,  9.642722,  9.648799,  9.648799,...
 y: array([ 9.962871,  9.962871,  9.962871,  9.962871,  9.962871,  9.962871,
        9.962871,  9.962871,  9.638172,  9.638172,  9.638172,  9.638172,
        9.638172,  9.639587,  9.639587,  9.639587,  9.64599 ,  9.64599 ,...

There was a bug corrected in the new implementation of the GCL model

Time comparison

New implementation is wrapped inside fusedwind


In [9]:
WD = range(0,360,1)

In [10]:
%%timeit
fP_WF = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    fgcl.wind_direction = wd
    fgcl.run()
    fP_WF[iwd] = fgcl.power


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-618732d476cb> in <module>()
----> 1 get_ipython().run_cell_magic(u'timeit', u'', u'fP_WF = np.zeros([len(WD)])\nfor iwd, wd in enumerate(WD):\n    fgcl.wind_direction = wd\n    fgcl.run()\n    fP_WF[iwd] = fgcl.power')

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2291             magic_arg_s = self.var_expand(line, stack_depth)
   2292             with self.builtin_trap:
-> 2293                 result = fn(magic_arg_s, cell)
   2294             return result
   2295 

<decorator-gen-59> in timeit(self, line, cell)

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    191     # but it's overkill for just that one bit of state.
    192     def magic_deco(arg):
--> 193         call = lambda f, *a, **k: f(*a, **k)
    194 
    195         if callable(arg):

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/IPython/core/magics/execution.pyc in timeit(self, line, cell)
   1035             number = 1
   1036             for _ in range(1, 10):
-> 1037                 time_number = timer.timeit(number)
   1038                 worst_tuning = max(worst_tuning, time_number / number)
   1039                 if time_number >= 0.2:

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/IPython/core/magics/execution.pyc in timeit(self, number)
    131         gc.disable()
    132         try:
--> 133             timing = self.inner(it, self.timer)
    134         finally:
    135             if gcold:

<magic-timeit> in inner(_it, _timer)

NameError: global name 'fgcl' is not defined

In [11]:
%%timeit
#%%prun -s cumulative #profiling
P_WF = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    P_WT,U_WT,CT_WT = GCLarsen(WS=8.0,z0=0.0001,TI=0.05,WD=wd,WF=HR1,NG=4,sup='lin')
    P_WF[iwd] = P_WT.sum()


1 loops, best of 3: 10.7 s per loop

Pandas


In [12]:
df=pd.DataFrame(data=P_WF, index=WD, columns=['P_WF'])

In [13]:
df.plot()


Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x10957f710>

WD uncertainty

Normally distributed wind direction uncertainty (reference wind direction, not for individual turbines).


In [14]:
P_WF_GAv8 = np.zeros([len(WD)])
P_WF_GAv16 = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    P_WT_GAv,U_WT,CT_WT = GCL_P_GaussQ_Norm_U_WD(meanWD=wd,stdWD=2.5,NG_P=8, WS=8.0,z0=0.0001,TI=0.05,WF=HR1,NG=4,sup='lin')
    P_WF_GAv8[iwd] = P_WT_GAv.sum()
    P_WT_GAv,U_WT,CT_WT = GCL_P_GaussQ_Norm_U_WD(meanWD=wd,stdWD=2.5,NG_P=16, WS=8.0,z0=0.0001,TI=0.05,WF=HR1,NG=4,sup='lin')
    P_WF_GAv16[iwd] = P_WT_GAv.sum()


---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-14-72fcceecdf2d> in <module>()
      4     P_WT_GAv,U_WT,CT_WT = GCL_P_GaussQ_Norm_U_WD(meanWD=wd,stdWD=2.5,NG_P=8, WS=8.0,z0=0.0001,TI=0.05,WF=HR1,NG=4,sup='lin')
      5     P_WF_GAv8[iwd] = P_WT_GAv.sum()
----> 6     P_WT_GAv,U_WT,CT_WT = GCL_P_GaussQ_Norm_U_WD(meanWD=wd,stdWD=2.5,NG_P=16, WS=8.0,z0=0.0001,TI=0.05,WF=HR1,NG=4,sup='lin')
      7     P_WF_GAv16[iwd] = P_WT_GAv.sum()

/Users/pe/git/FUSED-Wake/fusedwake/gcl/python/gcl.py in GCL_P_GaussQ_Norm_U_WD(WF, WS, meanWD, stdWD, NG_P, TI, z0, alpha, inflow, NG, sup, pars)
    603                             NG=NG,
    604                             sup=sup,
--> 605                             pars=pars)
    606         meanP_WT += wi[i]*P_WT*(1./np.sqrt(np.pi))
    607         meanU_WT += wi[i]*U_WT*(1./np.sqrt(np.pi))

/Users/pe/git/FUSED-Wake/fusedwake/gcl/python/gcl.py in GCLarsen(WF, WS, WD, TI, z0, alpha, inflow, NG, sup, pars)
    477         cU = U_WT[cWT]
    478         if cU>WF.WT[i].u_cutin:
--> 479             Ct[cWT] = WF.WT[i].get_CT(U_WT[cWT])
    480             P_WT[cWT] = WF.WT[i].get_P(U_WT[cWT])
    481         else:

/Users/pe/git/FUSED-Wake/fusedwake/WindTurbine.pyc in get_CT(self, u)
    126         """
    127         #return np.interp(u, self.ref_u, self.ref_CT)
--> 128         return ((u>=self.u_cutin)&(u<=self.u_cutout))*self.CTCI(u) + \
    129                 ((u<self.u_cutin)|(u>self.u_cutout))*self.CT_idle
    130 

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/scipy/interpolate/polyint.pyc in __call__(self, x)
     77         """
     78         x, x_shape = self._prepare_x(x)
---> 79         y = self._evaluate(x)
     80         return self._finish_y(y, x_shape)
     81 

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/scipy/interpolate/interpolate.pyc in _evaluate(self, x_new)
    495         #    The behavior is set by the bounds_error variable.
    496         x_new = asarray(x_new)
--> 497         out_of_bounds = self._check_bounds(x_new)
    498         y_new = self._call(self, x_new)
    499         if len(y_new) > 0:

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/scipy/interpolate/interpolate.pyc in _check_bounds(self, x_new)
    521 
    522         # !! Could provide more information about which values are out of bounds
--> 523         if self.bounds_error and below_bounds.any():
    524             raise ValueError("A value in x_new is below the interpolation "
    525                 "range.")

/Users/pe/.virtualenvs/fusedwake/lib/python2.7/site-packages/numpy/core/_methods.pyc in _any(a, axis, dtype, out, keepdims)
     36 
     37 def _any(a, axis=None, dtype=None, out=None, keepdims=False):
---> 38     return umr_any(a, axis, dtype, out, keepdims)
     39 
     40 def _all(a, axis=None, dtype=None, out=None, keepdims=False):

KeyboardInterrupt: 

In [ ]:
fig, ax = plt.subplots()
fig.set_size_inches([12,6])
ax.plot(WD,P_WF/(HR1.WT.get_P(8.0)*HR1.nWT),'-o', label='Pure python')
ax.plot(WD,fP_WF/(HR1.WT.get_P(8.0)*HR1.nWT),'-d', label='FUSED wrapper')
ax.plot(WD,P_WF_GAv16/(HR1.WT.get_P(8.0)*HR1.nWT),'-', label='Gauss Avg. FUSED wrapper, NG_P = 16')
ax.plot(WD,P_WF_GAv8/(HR1.WT.get_P(8.0)*HR1.nWT),'-', label='Gauss Avg. FUSED wrapper, NG_P = 8')
ax.set_xlabel('wd [deg]')
ax.set_ylabel('Wind farm efficiency [-]')
ax.set_title(HR1.name)
ax.legend(loc=3)
plt.savefig(HR1.name+'_Power_wd_360.pdf')

Uniformly distributed wind direction uncertainty (bin/sectors definition)


In [ ]:
P_WF_GA_u8 = np.zeros([len(WD)])
for iwd, wd in enumerate(WD):
    P_WT_GAv,U_WT,CT_WT = GCL_P_GaussQ_Uni_U_WD(meanWD=wd,U_WD=2.5,NG_P=8, WS=8.0,z0=0.0001,TI=0.05,WF=HR1,NG=4,sup='lin')
    P_WF_GA_u8[iwd] = P_WT_GAv.sum()

In [ ]:
fig, ax = plt.subplots()
fig.set_size_inches([12,6])
ax.plot(WD,fP_WF/(HR1.WT.get_P(8.0)*HR1.nWT),'-d', label='FUSED wrapper')
ax.plot(WD,P_WF_GAv8/(HR1.WT.get_P(8.0)*HR1.nWT),'-', label='Gauss Quad. Normal, NG_P = 8')
ax.plot(WD,P_WF_GA_u8/(HR1.WT.get_P(8.0)*HR1.nWT),'-', label='Gauss Quad. Uniform, NG_P = 8')
ax.set_xlabel('wd [deg]')
ax.set_ylabel('Wind farm efficiency [-]')
ax.set_title(HR1.name)
ax.legend(loc=3)
plt.savefig(HR1.name+'_Power_wd_360.pdf')

In [ ]:


In [ ]:


In [ ]:


In [ ]: