Add sowfa suite (NREL#2)

* fix layout to be same size

* start first examples

* add get points function

* initial example version

* update examples

* update examples

* bug fix

examples/example_input.json

@@ -22,10 +22,8 @@
       "layout_x": [
         0,
         756,
-        1512,
         0,
-        756,
-        1512
+        756
       ],
       "layout_y": [
         0.0,

examples/sowfa_power_suite/example_0000_get_sowfa_inflow_high.py

@@ -0,0 +1,99 @@
+#
+# Copyright 2019 NREL
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not use
+# this file except in compliance with the License. You may obtain a copy of the
+# License at http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software distributed
+# under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
+# CONDITIONS OF ANY KIND, either express or implied. See the License for the
+# specific language governing permissions and limitations under the License.
+#
+
+# See read the https://floris.readthedocs.io for documentation
+
+from floris.utilities import Vec3
+import floris.tools as wfct
+import floris.tools.visualization as vis
+import floris.tools.cut_plane as cp
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import pandas as pd
+
+# # Define a minspeed and maxspeed to use across visualiztions
+minspeed = 4.0
+maxspeed = 8.5
+
+for ti in ['low','hi']:
+
+    # Load the SOWFA case in
+    sowfa_root = '/Users/pfleming/Box Sync/sowfa_library/full_runs/near_wake'
+    inflow_case = '%s_no_turbine' % ti
+
+    si = wfct.sowfa_utilities.SowfaInterface(os.path.join(sowfa_root,inflow_case))
+
+    flow_origin = si.flow_data.origin
+    print("origin of saved flow field = ", flow_origin)
+
+    # Re-origin flow to 0
+    si.flow_data.x = si.flow_data.x + flow_origin.x1
+    si.flow_data.y = si.flow_data.y + flow_origin.x2
+    si.flow_data.origin = Vec3(0,0,0)
+
+
+    # Get the hub-height velocities at turbine locations
+    y_points = np.arange(200,1600,10.) 
+    x_points = np.ones_like(y_points) * 1000. 
+    z_points = np.ones_like(y_points) * 90.
+
+    u_points = si.flow_data.get_points_from_flow_data(x_points,y_points,z_points)
+
+    x_points_2 = np.ones_like(y_points) * 2000.
+    u_points_2 = si.flow_data.get_points_from_flow_data(x_points_2,y_points,z_points)
+
+    x_points_3 = np.ones_like(y_points) * 3000. 
+    u_points_3 = si.flow_data.get_points_from_flow_data(x_points_3,y_points,z_points)
+
+    # Grab a set of points to describe the flow field by
+    x_p = np.array([950,950,950,950,1500,1500,1500,1500,2000,2000,2000,2000])
+    y_p = np.array([300,700,1100,1500,300,700,1100,1500,300,700,1100,1500]) 
+    z_p = z_points = np.ones_like(x_p) * 90.
+    u_p = si.flow_data.get_points_from_flow_data(x_p,y_p,z_p)
+
+    # Save the flow points as a dataframe
+    df = pd.DataFrame({'x':x_p,
+                        'y':y_p,
+                        'z':z_p,
+                        'u':u_p})
+    df.to_pickle('flow_data_%s.p' % ti)
+
+    # Plot the SOWFA flow and turbines using the input information
+    fig, axarr = plt.subplots(2, 1, figsize=(5, 5))
+
+    ax = axarr[0]
+    sowfa_flow_data = si.flow_data
+    hor_plane = si.get_hor_plane(90)
+    wfct.visualization.visualize_cut_plane(
+        hor_plane, ax=ax, minSpeed=minspeed, maxSpeed=maxspeed)
+    ax.plot(x_points,y_points,color='k')
+    ax.plot(x_points_2,y_points,color='g')
+    ax.plot(x_points_3,y_points,color='r')
+    ax.scatter(x_p,y_p,color='m')
+
+    ax.set_title('SOWFA')
+    ax.set_ylabel('y location [m]')
+
+    ax = axarr[1]
+    ax.plot(y_points,u_points,color='k')
+    ax.plot(y_points,u_points_2,color='g')
+    ax.plot(y_points,u_points_3,color='r')
+    vis.reverse_cut_plane_x_axis_in_plot(ax)
+    ax.set_title('Looking downstream')
+    fig.suptitle(ti)
+
+
+
+
+plt.show()

examples/sowfa_power_suite/example_0001_set_up_first_case.py

@@ -0,0 +1,113 @@
+#
+# Copyright 2019 NREL
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not use
+# this file except in compliance with the License. You may obtain a copy of the
+# License at http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software distributed
+# under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
+# CONDITIONS OF ANY KIND, either express or implied. See the License for the
+# specific language governing permissions and limitations under the License.
+#
+
+# See read the https://floris.readthedocs.io for documentation
+
+import floris.tools as wfct
+import floris.tools.visualization as vis
+import floris.tools.cut_plane as cp
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import pandas as pd
+
+
+# Parameters
+ti = 'hi'  #hi/low
+
+# Load in SOWFA flow data and power results
+df_flow = pd.read_pickle('flow_data_%s.p' % ti)
+df_power = pd.read_pickle('data_sowfa.p')
+
+if ti == 'hi':
+    df_power = df_power[df_power.TI > 0.07]
+    het_wind_speed_gain = 1.03
+    hom_wind_speed = 8.25
+else:
+    df_power = df_power[df_power.TI < 0.07]
+    het_wind_speed_gain = 1.02
+    hom_wind_speed = 8.33
+
+df_power = df_power[df_power.yaw_0==0]
+
+# Iniitialize FLORIS
+fi = wfct.floris_interface.FlorisInterface("example_input.json")
+
+# Pick a random case and set FLORIS to match
+random_case = df_power.sample()
+print(random_case)
+x_locations = np.array(random_case.layout_x.values[0])
+y_locations = np.array(random_case.layout_y.values[0])
+fi.reinitialize_flow_field(wind_speed=[hom_wind_speed],layout_array =[x_locations,y_locations])
+
+yaw_array = np.array([random_case.yaw_0.values,random_case.yaw_1.values,random_case.yaw_2.values,random_case.yaw_3.values])
+sowfa_power_array = np.array([random_case.sowfa_power_0.values,random_case.sowfa_power_1.values,random_case.sowfa_power_2.values,random_case.sowfa_power_3.values])
+fi.calculate_wake(yaw_angles=yaw_array)
+
+# Get horizontal plane at default height (hub-height)
+
+
+# Plot and show homogenous flow
+fig, axarr = plt.subplots(3,2,sharex='row',sharey='row',figsize=(8,10))
+
+ax = axarr[0,0]
+hor_plane = fi.get_hor_plane()
+wfct.visualization.visualize_cut_plane(hor_plane, ax=ax)
+ax.set_title('Homogenous')
+
+# Compare homogenous results
+floris_power_array = np.array([p[0]/1000. for p in fi.get_turbine_power()])
+
+het_wind_speed = het_wind_speed_gain * df_flow.u.values
+
+ax = axarr[1,0]
+row_locs = [1,3]
+ax.plot(x_locations[row_locs],sowfa_power_array[row_locs],'.-',color='k',label='SOWFA')
+ax.plot(x_locations[row_locs],floris_power_array[row_locs],'.-',color='r',label='FLORIS')
+ax.set_title('Top Row')
+
+ax = axarr[2,0]
+row_locs = [0,2]
+ax.plot(x_locations[row_locs],sowfa_power_array[row_locs],'.-',color='k',label='SOWFA')
+ax.plot(x_locations[row_locs],floris_power_array[row_locs],'.-',color='r',label='FLORIS')
+ax.set_title('Bottom Row')
+
+# Redo to Het
+
+fi.reinitialize_flow_field( wind_speed=het_wind_speed.tolist(),#wind_direction=270. * np.ones_like(df_flow.u.values),
+                                wind_layout=[df_flow.x.values, df_flow.y.values])
+
+fi.calculate_wake()
+
+# Compare homogenous results
+floris_power_array = np.array([p[0]/1000. for p in fi.get_turbine_power()])
+
+ax = axarr[0,1]
+hor_plane = fi.get_hor_plane()
+wfct.visualization.visualize_cut_plane(hor_plane, ax=ax)
+ax.set_title('Heterogenous')
+
+ax = axarr[1,1]
+row_locs = [1,3]
+ax.plot(x_locations[row_locs],sowfa_power_array[row_locs],'.-',color='k',label='SOWFA')
+ax.plot(x_locations[row_locs],floris_power_array[row_locs],'.-',color='r',label='FLORIS')
+ax.set_title('Top Row')
+
+ax = axarr[2,1]
+row_locs = [0,2]
+ax.plot(x_locations[row_locs],sowfa_power_array[row_locs],'.-',color='k',label='SOWFA')
+ax.plot(x_locations[row_locs],floris_power_array[row_locs],'.-',color='r',label='FLORIS')
+ax.set_title('Bottom Row')
+
+plt.show()
+

examples/sowfa_power_suite/example_0002_compare_to_sowfa.py

@@ -0,0 +1,110 @@
+#
+# Copyright 2019 NREL
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not use
+# this file except in compliance with the License. You may obtain a copy of the
+# License at http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software distributed
+# under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR
+# CONDITIONS OF ANY KIND, either express or implied. See the License for the
+# specific language governing permissions and limitations under the License.
+#
+
+# See read the https://floris.readthedocs.io for documentation
+
+import floris.tools as wfct
+import floris.tools.visualization as vis
+import floris.tools.cut_plane as cp
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+import pandas as pd
+
+
+D = 126
+
+fig, axarr = plt.subplots(2,3,sharex=True,sharey=True)
+
+
+# Load in results
+df_power_full = pd.read_pickle('data_sowfa.p')
+
+# List simulation options
+ti_vals = ['low','hi']
+yaw_values = [0,10,20]
+x_locations_unique = sorted(df_power_full.layout_x.unique())
+
+# These dont change
+y_locations = np.array(df_power_full.layout_y.values[0])
+
+
+# Iniitialize FLORIS
+fi = wfct.floris_interface.FlorisInterface("example_input.json")
+fi_b = wfct.floris_interface.FlorisInterface("example_input.json")
+fi_b.floris.farm.set_wake_model('blondel')
+
+
+
+for ti_idx, ti in enumerate(ti_vals):
+    if ti == 'hi':
+        df_power_ti = df_power_full[df_power_full.TI > 0.07]
+        ti_val = 0.09
+        wind_speed = 8.25
+    else:
+        df_power_ti = df_power_full[df_power_full.TI < 0.07]
+        ti_val = 0.065
+        wind_speed = 8.33
+
+    for yaw_idx, yaw in enumerate(yaw_values):
+
+        df_power_yaw = df_power_ti[df_power_ti.yaw_0==yaw]
+        d_array = []
+        results_sowfa_1 = []
+        results_sowfa_2 = []
+        gauss = []
+        blondel = []
+
+        for x_locations in x_locations_unique:
+
+
+            df_inner = df_power_yaw[df_power_yaw.layout_x==x_locations]
+
+            yaw_array = np.array([df_inner.yaw_0.values,df_inner.yaw_1.values,df_inner.yaw_2.values,df_inner.yaw_3.values])
+            sowfa_power_array = np.array([df_inner.sowfa_power_0.values,df_inner.sowfa_power_1.values,df_inner.sowfa_power_2.values,df_inner.sowfa_power_3.values])
+
+            # Set up FLORIS and get powers
+            fi.reinitialize_flow_field(wind_speed=[wind_speed],turbulence_intensity=[ti_val],layout_array =[x_locations,y_locations])
+            fi.calculate_wake(yaw_angles=yaw_array)
+            floris_power_array = np.array([p[0]/1000. for p in fi.get_turbine_power()])
+
+            # Repeat BLONDEL
+            fi_b.reinitialize_flow_field(wind_speed=[wind_speed],turbulence_intensity=[ti_val],layout_array =[x_locations,y_locations])
+            fi_b.calculate_wake(yaw_angles=yaw_array)
+            floris_b_power_array = np.array([p[0]/1000. for p in fi_b.get_turbine_power()])
+
+
+            # Save all the results
+            d_loc = (x_locations[2] - x_locations[0])/D
+            d_array.append(d_loc)
+            results_sowfa_1.append(sowfa_power_array[2])
+            results_sowfa_2.append(sowfa_power_array[3])
+            gauss.append(floris_power_array[3]) # Identical, just pick one
+            blondel.append(floris_b_power_array[3]) # Identical, just pick one
+
+        ax = axarr[ti_idx, yaw_idx]
+        ax.plot(d_array,results_sowfa_1, color='k',marker='o',ls='None',label='SOWFA 1')
+        ax.plot(d_array,results_sowfa_2, color='k',marker='x',ls='None',label='SOWFA 2')
+        ax.plot(d_array,gauss, color='g',marker='.',label="gauss")
+        ax.plot(d_array,blondel, color='violet',marker='.',label="blondel")
+        ax.set_title('%s TI, yaw=%d' % (ti,yaw))
+        ax.grid(True)
+        if (ti_idx==0) and (yaw_idx==0):
+            ax.legend()
+
+
+
+
+
+plt.show()
+