diff --git a/demo/documented/studies/powercurve_validation/CT_calibration.py b/demo/documented/studies/powercurve_validation/CT_calibration.py new file mode 100644 index 00000000..31d5e8a1 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/CT_calibration.py @@ -0,0 +1,146 @@ +import yaml + +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +import seaborn + +################################################################################ +# this is a brief script for calibration of CT. +# +# data comes from a run at 7.0 m/s (region II) with all calibration factors set +# to zero. a convergence study is run with single_turbine.yaml, and the results +# are stored and saved in CTcalibration_data.csv. +# +# we target matching the reference CT at 7.0 m/s with the asymptotic simulation. +# we assume geometric convergence of thrust w.r.t. dx. this allows the hand-fit +# of the convergence data for an estimate of the asymptotic thrust, shown below. +################################################################################ + +# simulation properties +rho_fluid = 1.225 # kg/m^3 +diameter_rotor = 130.0 # m +area_rotor = np.pi/4*diameter_rotor**2 # m^2 +V_inf = 8.0 # m/s +Tf_fluid = 0.5*rho_fluid*area_rotor*V_inf**2/1e3 +P_fluid = 0.5*rho_fluid*area_rotor*V_inf**3/1e6 +volume_domain = (780.0 - -390.0)*(390.0 - -390.0)*(520.0 - 0.04) + +# load, amend calibration simulation data +fn_CT_data = "CTcalibration_data.csv" +df_data = pd.read_csv(fn_CT_data) +df_data["dx_elem"] = volume_domain/df_data.Nelem + +# fit loge convergence with linear model +mb_Tf = np.polyfit( + np.log10(df_data.dx_elem.to_numpy()[:-1]), + np.log10(np.abs(rho_fluid*df_data.Tf.to_numpy()[:-1] - rho_fluid*df_data.Tf.to_numpy()[-1])), + deg=1, +) +mb_P = np.polyfit( + np.log10(df_data.dx_elem.to_numpy()[:-1]), + np.log10(np.abs(rho_fluid*df_data.P.to_numpy()[:-1] - rho_fluid*df_data.P.to_numpy()[-1])), + deg=1, +) +# i.e.: +# log10(|J-Jinf|) = m*log(dx) + b +# |J-Jinf| = dx**m * 10**b +# Jinf = J -/+ dx**m * 10**b + +# estimates of the true thrust based on each sample +Tf_inf_est = rho_fluid*df_data.Tf + df_data.dx_elem**mb_Tf[0] * 10**mb_Tf[1] +signpattern_P = np.ones_like(df_data.dx_elem) # np.array([1.0, 1.0, 1.0, 1.0]) +P_inf_est = rho_fluid*df_data.P + signpattern_P*df_data.dx_elem**mb_P[0] * 10**mb_P[1] + +# plot the convergence behavior for validation +fig, axes = plt.subplots(2, sharex=True) + +# thrust vs. discretization length +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + rho_fluid*df_data.Tf.to_numpy(), +) +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + Tf_inf_est.to_numpy(), + "--", label="corrected value", +) +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + Tf_inf_est.to_numpy()[-1]*np.ones_like(df_data.dx_elem.to_numpy()), + "--", label="estimated asymptote", +) +axes[0].set_xlabel("discretization length, $\\Delta x$ (m)") +axes[0].set_ylabel("thrust force, $T$ (kN)") +axes[0].legend() + +# apparent error vs. discretization length +axes[1].loglog( + df_data.dx_elem.to_numpy()[:-1], + np.abs(rho_fluid*df_data.Tf.to_numpy()[:-1] - rho_fluid*df_data.Tf.to_numpy()[-1]), +) +axes[1].loglog( + df_data.dx_elem.to_numpy()[:-1], + 10**(mb_Tf[0]*np.log10(df_data.dx_elem.to_numpy()[:-1]) + mb_Tf[1]), + "--", label="line of best fit", +) +axes[1].set_xlabel("discretization length, $\\Delta x$ (m)") +axes[1].set_ylabel("apparent thrust error, $e_T$ (kN)") +axes[1].legend() + +# plot the convergence behavior for validation +fig, axes = plt.subplots(2, sharex=True) + +# power vs. discretization length +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + rho_fluid*df_data.P.to_numpy(), +) +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + P_inf_est.to_numpy(), + "--", label="corrected value", +) +axes[0].semilogx( + df_data.dx_elem.to_numpy(), + P_inf_est.to_numpy()[-1]*np.ones_like(df_data.dx_elem.to_numpy()), + "--", label="estimated asymptote", +) +axes[0].set_xlabel("discretization length, $\\Delta x$ (m)") +axes[0].set_ylabel("power, $P$ (MW)") +axes[0].legend() + +# apparent error vs. discretization length +axes[1].loglog( + df_data.dx_elem.to_numpy()[:-1], + np.abs(rho_fluid*df_data.P.to_numpy()[:-1] - rho_fluid*df_data.P.to_numpy()[-1]), +) +axes[1].loglog( + df_data.dx_elem.to_numpy()[:-1], + 10**(mb_P[0]*np.log10(df_data.dx_elem.to_numpy()[:-1]) + mb_P[1]), + "--", label="line of best fit", +) +axes[1].set_xlabel("discretization length, $\\Delta x$ (m)") +axes[1].set_ylabel("power error, $e_P$ (MW)") +axes[1].legend() + +# get the reference power curve we want to match +with open("writeup/FLORIS_IEA-3p4-130-RWT.yaml", "r") as f_yaml: + data_refPC = yaml.safe_load(f_yaml) +V_refPC = np.array(data_refPC["power_thrust_table"]["wind_speed"]) +CT_refPC = np.array(data_refPC["power_thrust_table"]["thrust_coefficient"]) +P_refPC = np.array(data_refPC["power_thrust_table"]["power"])/1e3 + +# print out the CT values we see +print(f"estimated asymptotic CT: {Tf_inf_est.to_numpy()[-1]/Tf_fluid}") +print(f"reference powercurve CT: {np.interp(V_inf, V_refPC, CT_refPC)}") +print(f"correction: {np.interp(V_inf, V_refPC, CT_refPC)/(Tf_inf_est.to_numpy()[-1]/Tf_fluid)}\n") + +print(f"estimated asymptotic CP: {P_inf_est.to_numpy()[-1]/P_fluid}") +print(f"reference powercurve CP: {np.interp(V_inf, V_refPC, P_refPC)/P_fluid}") +print(f"correction: {np.interp(V_inf, V_refPC, P_refPC)/(P_inf_est.to_numpy()[-1])}") + +# show and exit +plt.show() + +### diff --git a/demo/documented/studies/powercurve_validation/CTcalibration_data.csv b/demo/documented/studies/powercurve_validation/CTcalibration_data.csv new file mode 100644 index 00000000..4044b937 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/CTcalibration_data.csv @@ -0,0 +1,5 @@ +Nelem,Ndof,V,P,Tf +63744,46048,8.0,1.404740050380867578e+00,3.039104256841761185e+02 +158074,112160,8.0,1.500207951806328710e+00,3.177386680432255730e+02 +451005,314080,8.0,1.553362153286078184e+00,3.254239347945710961e+02 +1110316,765984,8.0,1.586169749666255901e+00,3.300853155380664816e+02 \ No newline at end of file diff --git a/demo/documented/studies/powercurve_validation/multi_turbine.yaml b/demo/documented/studies/powercurve_validation/multi_turbine.yaml new file mode 100644 index 00000000..37057777 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/multi_turbine.yaml @@ -0,0 +1,79 @@ +general: + name: 3D_Wind_Farm # Name of the output folder + output: ["mesh","initial_guess","height","turbine_force","solution"] + +wind_farm: + + ######################### Grid Wind Farm ######################### + type: grid # | + grid_rows: 1 # Number of rows | - + grid_cols: 3 # Number of columns | - + x_spacing: 390 # horizontal spacing | m 390, 520, 650, 780, 910 + y_spacing: 0 # vertical spacing | m + ################################################################### + +turbines: + + ######################### Turbine Properties ######################## + type: power_disk # | + RD: 130.0 # Turbine Diameter | m + thickness: 13.0 # Effective Thickness | m + CPprime0: 1.022 # Region 2 Power Coeff | - + # CPprime0: 1.043 # Region 2 Power Coeff | - + # CPprime0: 1.360 # Region 2 Power Coeff | - + CTprime0: 1.319 # Region 2 Thrust Coeff | - + # CTprime0: 1.350 # Region 2 Thrust Coeff | - + # CTprime0: 1.431 # Region 2 Thrust Coeff | - + # CTprime0: 1.639 # Region 2 Thrust Coeff | - + Prated: 2570282.74287 # Specific Rated Power | W + # Prated: 2532257.89912 # Specific Rated Power | W + # Prated: 3077551.02040 # Specific Rated Power | W + Trated: 368163.26530 # Specific Thrust | N + force: sine # radial force distribution | - + HH: 110 # Hub Height | m + yaw: 0.0 # Yaw | rads + ################################################################### + +domain: + + ########################### Box Domain ############################ + type: box # | + x_range: [-390, 1950] # x-range of the domain | m 1950, 2340, 2730, 3120, 3510 + y_range: [-390, 390] # y-range of the domain | m + z_range: [0.04, 520] # z-range of the domain | m + nx: 54 # Number of x-nodes | - 54, 63, 72, 81, 90 + ny: 18 # Number of y-nodes | - + nz: 9 # Number of z-nodes | - + ################################################################### + +refine: + warp_type: split + warp_percent: 0.85 # percent of cells moved | - + warp_height: 240 # move cell below this value | m + refine_custom: + 1: + type: box + x_range: [-130, 1300] # 1300, 1690, 2080, 2470, 2860 + y_range: [-130, 130] + z_range: [0, 240] + turbine_num: 1 # number of turbine refinements| - + turbine_factor: 1.25 # turbine radius multiplier | - + +function_space: + type: linear + +boundary_conditions: + vel_profile: log + HH_vel: 8.0 + k: 0.4 + inflow_angle: 270.0 + +problem: + type: stabilized + viscosity: 5 + lmax: 50 + +solver: + type: steady + save_power: true + save_thrust: true diff --git a/demo/documented/studies/powercurve_validation/multi_turbine_fmt.yaml b/demo/documented/studies/powercurve_validation/multi_turbine_fmt.yaml new file mode 100644 index 00000000..8fe56f96 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/multi_turbine_fmt.yaml @@ -0,0 +1,79 @@ +general: + name: CASENAME_TO_REPLACE # Name of the output folder + output: ["mesh","initial_guess","height","turbine_force","solution"] + +wind_farm: + + ######################### Grid Wind Farm ######################### + type: grid # | + grid_rows: 1 # Number of rows | - + grid_cols: 3 # Number of columns | - + x_spacing: 910 # horizontal spacing | m 390, 520, 650, 780, 910 + y_spacing: 0 # vertical spacing | m + ################################################################### + +turbines: + + ######################### Turbine Properties ######################## + type: power_disk # | + RD: 130.0 # Turbine Diameter | m + thickness: 13.0 # Effective Thickness | m + CPprime0: 1.022 # Region 2 Power Coeff | - + # CPprime0: 1.043 # Region 2 Power Coeff | - + # CPprime0: 1.360 # Region 2 Power Coeff | - + CTprime0: 1.319 # Region 2 Thrust Coeff | - + # CTprime0: 1.350 # Region 2 Thrust Coeff | - + # CTprime0: 1.431 # Region 2 Thrust Coeff | - + # CTprime0: 1.639 # Region 2 Thrust Coeff | - + Prated: 2570282.74287 # Specific Rated Power | W + # Prated: 2532257.89912 # Specific Rated Power | W + # Prated: 3077551.02040 # Specific Rated Power | W + Trated: 368163.26530 # Specific Thrust | N + force: sine # radial force distribution | - + HH: 110.0 # Hub Height | m + yaw: 0.0 # Yaw | rads + ################################################################### + +domain: + + ########################### Box Domain ############################ + type: box # | + x_range: [-390, 3510] # x-range of the domain | m 1950, 2340, 2730, 3120, 3510 + y_range: [-390, 390] # y-range of the domain | m + z_range: [0.04, 520] # z-range of the domain | m + nx: 90 # Number of x-nodes | - 54, 63, 72, 81, 90 + ny: 18 # Number of y-nodes | - + nz: 9 # Number of z-nodes | - + ################################################################### + +refine: + warp_type: split + warp_percent: 0.85 # percent of cells moved | - + warp_height: 240 # move cell below this value | m + refine_custom: + 1: + type: box + x_range: [-130, 2860] # 1300, 1690, 2080, 2470, 2860 + y_range: [-130, 130] + z_range: [0, 240] + turbine_num: 1 # number of turbine refinements| - + turbine_factor: 1.25 # turbine radius multiplier | - + +function_space: + type: linear + +boundary_conditions: + vel_profile: log + HH_vel: VINF_TO_REPLACE + k: 0.4 + inflow_angle: 270.0 + +problem: + type: stabilized + viscosity: 5 + lmax: 50 + +solver: + type: steady + save_power: true + save_thrust: true diff --git a/demo/documented/studies/powercurve_validation/single_turbine.yaml b/demo/documented/studies/powercurve_validation/single_turbine.yaml new file mode 100644 index 00000000..a090b3f3 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/single_turbine.yaml @@ -0,0 +1,79 @@ +general: + name: 3D_Wind_Farm # Name of the output folder + output: ["mesh","initial_guess","height","turbine_force","solution"] + +wind_farm: + + ######################### Grid Wind Farm ######################### + type: grid # | + grid_rows: 1 # Number of rows | - + grid_cols: 1 # Number of columns | - + x_spacing: 0 # horizontal spacing | m + y_spacing: 0 # vertical spacing | m + ################################################################### + +turbines: + + ######################### Turbine Properties ######################## + type: power_disk # | + RD: 130.0 # Turbine Diameter | m + thickness: 13.0 # Effective Thickness | m + CPprime0: 1.022 # Region 2 Power Coeff | - + # CPprime0: 1.043 # Region 2 Power Coeff | - + # CPprime0: 1.360 # Region 2 Power Coeff | - + CTprime0: 1.319 # Region 2 Thrust Coeff | - + # CTprime0: 1.350 # Region 2 Thrust Coeff | - + # CTprime0: 1.431 # Region 2 Thrust Coeff | - + # CTprime0: 1.639 # Region 2 Thrust Coeff | - + Prated: 2570282.74287 # Specific Rated Power | W + # Prated: 2532257.89912 # Specific Rated Power | W + # Prated: 3077551.02040 # Specific Rated Power | W + Trated: 368163.26530 # Specific Thrust | N + force: sine # radial force distribution | - + HH: 110 # Hub Height | m + yaw: 0.0 # Yaw | rads + ################################################################### + +domain: + + ########################### Box Domain ############################ + type: box # | + x_range: [-390, 780] # x-range of the domain | m + y_range: [-390, 390] # y-range of the domain | m + z_range: [0.04, 520] # z-range of the domain | m + nx: 54 # 19 # 27 # 38 # 54 # Number of x-nodes | - + ny: 36 # 13 # 18 # 25 # 36 # Number of y-nodes | - + nz: 18 # 6 # 9 # 13 # 18 # Number of z-nodes | - + ################################################################### + +refine: + warp_type: split + warp_percent: 0.85 # percent of cells moved | - + warp_height: 240 # move cell below this value | m + refine_custom: + 1: + type: box + x_range: [-130, 130] + y_range: [-130, 130] + z_range: [0, 240] + turbine_num: 1 # number of turbine refinements| - + turbine_factor: 1.25 # turbine radius multiplier | - + +function_space: + type: linear + +boundary_conditions: + vel_profile: log + HH_vel: 15.0 # 8.0 + k: 0.4 + inflow_angle: 270.0 + +problem: + type: stabilized + viscosity: 5 + lmax: 50 + +solver: + type: steady + save_power: true + save_thrust: true diff --git a/demo/documented/studies/powercurve_validation/single_turbine_fmt.yaml b/demo/documented/studies/powercurve_validation/single_turbine_fmt.yaml new file mode 100644 index 00000000..8e5f18d8 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/single_turbine_fmt.yaml @@ -0,0 +1,79 @@ +general: + name: CASENAME_TO_REPLACE # Name of the output folder + output: ["mesh","initial_guess","height","turbine_force","solution"] + +wind_farm: + + ######################### Grid Wind Farm ######################### + type: grid # | + grid_rows: 1 # Number of rows | - + grid_cols: 1 # Number of columns | - + x_spacing: 0 # horizontal spacing | m + y_spacing: 0 # vertical spacing | m + ################################################################### + +turbines: + + ######################### Turbine Properties ######################## + type: power_disk # | + RD: 130.0 # Turbine Diameter | m + thickness: 13.0 # Effective Thickness | m + CPprime0: 1.022 # Region 2 Power Coeff | - + # CPprime0: 1.043 # Region 2 Power Coeff | - + # CPprime0: 1.360 # Region 2 Power Coeff | - + CTprime0: 1.319 # Region 2 Thrust Coeff | - + # CTprime0: 1.350 # Region 2 Thrust Coeff | - + # CTprime0: 1.431 # Region 2 Thrust Coeff | - + # CTprime0: 1.639 # Region 2 Thrust Coeff | - + Prated: 2570282.74287 # Specific Rated Power | W + # Prated: 2532257.89912 # Specific Rated Power | W + # Prated: 3077551.02040 # Specific Rated Power | W + Trated: 368163.26530 # Specific Thrust | N + force: sine # radial force distribution | - + HH: 110.0 # Hub Height | m + yaw: 0.0 # Yaw | rads + ################################################################### + +domain: + + ########################### Box Domain ############################ + type: box # | + x_range: [-390, 780] # x-range of the domain | m + y_range: [-390, 390] # y-range of the domain | m + z_range: [0.04, 520] # z-range of the domain | m + nx: 27 # 19 # 27 # 38 # Number of x-nodes | - + ny: 18 # 13 # 18 # 25 # Number of y-nodes | - + nz: 9 # 6 # 9 # 13 # Number of z-nodes | - + ################################################################### + +refine: + warp_type: split + warp_percent: 0.85 # percent of cells moved | - + warp_height: 240 # move cell below this value | m + refine_custom: + 1: + type: box + x_range: [-130, 130] + y_range: [-130, 130] + z_range: [0, 240] + turbine_num: 1 # number of turbine refinements| - + turbine_factor: 1.25 # turbine radius multiplier | - + +function_space: + type: linear + +boundary_conditions: + vel_profile: log + HH_vel: VINF_TO_REPLACE + k: 0.4 + inflow_angle: 270.0 + +problem: + type: stabilized + viscosity: 5 + lmax: 50 + +solver: + type: steady + save_power: true + save_thrust: true diff --git a/demo/documented/studies/powercurve_validation/writeup/FLORIS_IEA-3p4-130-RWT.yaml b/demo/documented/studies/powercurve_validation/writeup/FLORIS_IEA-3p4-130-RWT.yaml new file mode 100644 index 00000000..a2328e68 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/FLORIS_IEA-3p4-130-RWT.yaml @@ -0,0 +1,174 @@ +TSR: 8.01754386 +hub_height: 110.0 +power_thrust_table: + cosine_loss_exponent_tilt: 1.88 + cosine_loss_exponent_yaw: 1.88 + power: + - 0.0 + - 0.0 + - 55.48565038827394 + - 132.60594949759707 + - 229.20900366592642 + - 342.71254323810695 + - 468.6008480810116 + - 602.8982926864225 + - 740.4915214111589 + - 880.4618007632128 + - 1020.8388498653877 + - 1156.4922707088438 + - 1283.0654945862623 + - 1396.4539297233146 + - 1492.9570847061225 + - 1569.412461109748 + - 1623.3079057292161 + - 1652.8696585473415 + - 1657.123938275946 + - 1678.4995701413166 + - 1726.1238569674936 + - 1801.3513158108537 + - 1906.3000805443924 + - 2043.8824757534885 + - 2217.8441255643593 + - 2432.8096396619285 + - 2694.3325007912877 + - 3008.946421349339 + - 3384.21515097733 + - 3622.320140431238 + - 3622.3451711828 + - 3622.3459040961607 + - 3622.3452257887666 + - 3622.3451528700944 + - 3622.357444012131 + - 3622.3451658123345 + - 3622.346991064941 + - 3622.3544214545163 + - 3622.3813238557404 + - 3622.3451993631197 + - 3622.3453899482156 + - 3622.3467481326484 + - 3622.3475934972635 + - 3622.347589398676 + - 3622.3647319812258 + - 3622.3816921026446 + - 3622.350066618316 + - 3622.345086853103 + - 3622.3458511887966 + - 3622.3493778845755 + - 3622.3452304136345 + - 3622.4579153703894 + - 0.0 + - 0.0 + ref_air_density: 1.225 + ref_tilt: 5 + thrust_coefficient: + - 0.0 + - 0.0 + - 0.8140283854894109 + - 0.8041771483638767 + - 0.7982113770373608 + - 0.7945189424815273 + - 0.791309102845844 + - 0.786112395855186 + - 0.7760330855568125 + - 0.766405559051694 + - 0.766405559051694 + - 0.7664055590516942 + - 0.7664055590516942 + - 0.766405559051694 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516942 + - 0.766405559051694 + - 0.7664055590516943 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516939 + - 0.766405559051694 + - 0.7664055590516943 + - 0.7664055590516939 + - 0.7664055590516944 + - 0.766405559051694 + - 0.7664055590516939 + - 0.7664055590516942 + - 0.8067632947740724 + - 0.5306091101592825 + - 0.4466865448824549 + - 0.3801106706360693 + - 0.3254312907353938 + - 0.2798429968252726 + - 0.2414778385628342 + - 0.20904693927376936 + - 0.1815926995732083 + - 0.1582789238889643 + - 0.13843565625820317 + - 0.12141358275643703 + - 0.10691750603387402 + - 0.09454669711591918 + - 0.08392954309933542 + - 0.07479548070995087 + - 0.06691762496259271 + - 0.060104597639619645 + - 0.05418408689005904 + - 0.04906125901748401 + - 0.044599887927284615 + - 0.04063092405931261 + - 0.03722282641218611 + - 0.0 + - 0.0 + wind_speed: + - 0.0 + - 2.9999 + - 3.0 + - 3.539159827293725 + - 4.048452625194753 + - 4.526747636329809 + - 4.9729829261640415 + - 5.386167740761253 + - 5.76538470650541 + - 6.109791866899474 + - 6.418624551919384 + - 6.691197075772739 + - 6.9269042592927175 + - 7.125222773587183 + - 7.285712301959686 + - 7.408016517522657 + - 7.491863874332182 + - 7.537068210287934 + - 7.543529160459564 + - 7.575825940997163 + - 7.646808745314545 + - 7.756319973835029 + - 7.9041164842451614 + - 8.089870131331459 + - 8.313168495544858 + - 8.57351579867522 + - 8.870334004602977 + - 9.202964102683834 + - 9.570667570917234 + - 9.812675420388173 + - 10.407952984173718 + - 10.875675946194342 + - 11.374758439769723 + - 11.904092376955541 + - 12.462502504037497 + - 13.048749010888468 + - 13.661530283657028 + - 14.299485794675704 + - 14.96119912317263 + - 15.645201100079884 + - 16.349973069956214 + - 17.0739502627819 + - 17.815525268139492 + - 18.57305160406695 + - 19.344847372659316 + - 20.12919899430266 + - 20.924365012249385 + - 21.728579959087647 + - 22.54005827652058 + - 23.356998279752162 + - 24.177586157678096 + - 25.0 + - 25.01 + - 50.0 +rotor_diameter: 130.0 +turbine_type: IEA-3.4-130-RWT diff --git a/demo/documented/studies/powercurve_validation/writeup/FLORIS_case.yaml b/demo/documented/studies/powercurve_validation/writeup/FLORIS_case.yaml new file mode 100644 index 00000000..3fa8d73a --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/FLORIS_case.yaml @@ -0,0 +1,90 @@ + +name: WindSE 3-turbine validation +description: Three turbines for WindSE power curve sanity checking +floris_version: v4 + +logging: + console: + enable: true + level: WARNING # Can be one of "CRITICAL", "ERROR", "WARNING", "INFO", "DEBUG". + file: + enable: false + level: WARNING # Can be one of "CRITICAL", "ERROR", "WARNING", "INFO", "DEBUG". + +solver: + type: turbine_grid # Can be one of: "turbine_grid", "flow_field_grid", "flow_field_planar_grid". + turbine_grid_points: 3 + +farm: + layout_x: + - 0.0 + - 975.0 + - 1950.0 + layout_y: + - 0.0 + - 0.0 + - 0.0 + turbine_type: + - !include FLORIS_IEA-3p4-130-RWT.yaml + +flow_field: + air_density: 1.225 + reference_wind_height: -1 + turbulence_intensities: + - 0.06 + wind_directions: + - 270.0 + wind_shear: 0.12 + wind_speeds: + - 8.0 + wind_veer: 0.0 + multidim_conditions: + Tp: 2.5 + Hs: 3.01 + +wake: + model_strings: + combination_model: sosfs + deflection_model: gauss + turbulence_model: crespo_hernandez + velocity_model: gauss + enable_secondary_steering: true + enable_yaw_added_recovery: true + enable_active_wake_mixing: false + enable_transverse_velocities: true + wake_deflection_parameters: + gauss: + ad: 0.0 + alpha: 0.58 + bd: 0.0 + beta: 0.077 + dm: 1.0 + ka: 0.38 + kb: 0.004 + jimenez: + ad: 0.0 + bd: 0.0 + kd: 0.05 + wake_velocity_parameters: + cc: + a_s: 0.179367259 + b_s: 0.0118889215 + c_s1: 0.0563691592 + c_s2: 0.13290157 + a_f: 3.11 + b_f: -0.68 + c_f: 2.41 + alpha_mod: 1.0 + gauss: + alpha: 0.58 + beta: 0.077 + ka: 0.38 + kb: 0.004 + jensen: + we: 0.05 + wake_turbulence_parameters: + crespo_hernandez: + initial: 0.1 + constant: 0.5 + ai: 0.8 + downstream: -0.32 \ No newline at end of file diff --git a/demo/documented/studies/powercurve_validation/writeup/floris.csv b/demo/documented/studies/powercurve_validation/writeup/floris.csv new file mode 100644 index 00000000..7523b10f --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/floris.csv @@ -0,0 +1,122 @@ +V,P1,T1,P2,T2,P3,T3,Ptot,Ttot +0.25,0.0,5.081157961792399e-05,0.0,5.081157961792399e-05,0.0,5.081157961792399e-05,0.0,5.081157961792399e-05 +0.5,0.0,0.00020324631847169595,0.0,0.00020324631847169595,0.0,0.00020324631847169595,0.0,0.00020324631847169595 +0.75,0.0,0.0004573042165613159,0.0,0.0004573042165613159,0.0,0.0004573042165613159,0.0,0.0004573042165613159 +1.0,0.0,0.0008129852738867838,0.0,0.0008129852738867838,0.0,0.0008129852738867838,0.0,0.0008129852738867838 +1.25,0.0,0.0012702894904480996,0.0,0.0012702894904480996,0.0,0.0012702894904480996,0.0,0.0012702894904480996 +1.5,0.0,0.0018292168662452635,0.0,0.0018292168662452635,0.0,0.0018292168662452635,0.0,0.0018292168662452635 +1.75,0.0,0.0024897674012782757,0.0,0.0024897674012782757,0.0,0.0024897674012782757,0.0,0.0024897674012782757 +2.0,0.0,0.0032519410955471353,0.0,0.0032519410955471353,0.0,0.0032519410955471353,0.0,0.0032519410955471353 +2.25,0.0,0.004115737949051843,0.0,0.004115737949051843,0.0,0.004115737949051843,0.0,0.004115737949051843 +2.5,0.0,0.0050811579617923985,0.0,0.0050811579617923985,0.0,0.0050811579617923985,0.0,0.0050811579617923985 +2.75,0.0,0.006148201133768802,0.0,0.006148201133768802,0.0,0.006148201133768802,0.0,0.006148201133768803 +3.0,0.0,0.007316867464981054,0.0,0.007316867464981054,0.0,0.007316867464981054,0.0,0.007316867464981054 +3.25,0.0901693575205025,69.52144565758395,0.0,0.008587156955429154,0.06875411039666096,69.75635190057609,0.05297448930572115,46.428794905038494 +3.5,0.12584608527593116,80.17455806429254,0.060593437228258405,81.00467442997646,0.06563378938592859,80.94055321895931,0.08402443729670606,80.70659523774277 +3.75,0.1709522373556338,91.66766527464611,0.0886533603647008,92.58027710071352,0.09453002756410211,92.49445515052608,0.11804520842814559,92.24746584196191 +4.0,0.21826269586179753,103.91738701701684,0.11655278414611135,104.87220699263722,0.12331891303484634,104.75978144913525,0.15271146434758506,104.5164584862631 +4.25,0.2747040048409826,116.9963869317817,0.14808177318305457,117.94942459654288,0.1583240766270319,117.85654209087265,0.19370328488368968,117.6007845397324 +4.5,0.33389387518793223,130.84850009464117,0.1846717654894301,131.86195456291392,0.19619686286095211,131.7447814264561,0.23825416784610484,131.4850786946704 +4.75,0.40259370750007256,145.45838812520546,0.22111615948745533,146.50742931517792,0.23516514148930237,146.3802134928729,0.28629166949227675,146.11534364441877 +5.0,0.47362141525102747,160.79117658029193,0.2645419992852367,161.99990561207412,0.2822399613242924,161.88288850732934,0.3401344586201855,161.55799023323178 +5.25,0.5546908661961891,176.56933097905352,0.3100716050617642,178.27300256695568,0.3291174013778394,178.13416601384597,0.39795995754526414,177.65883318661838 +5.5,0.6395825913201537,192.66620581712206,0.3582219490106552,195.29723793078637,0.38204669692449317,195.1478435295249,0.459950412418434,194.3704290924778 +5.75,0.730081457347835,208.79743715413863,0.4129445146870713,213.07999450524133,0.437249783793957,212.91341684707635,0.5267585852762877,211.5969495021521 +6.0,0.8301984839117598,225.31936386004537,0.467780168543108,231.60227747817024,0.49589743853302753,231.2870139840439,0.597958696995965,229.40288510741985 +6.25,0.9376179784698675,243.38923177109183,0.5299232621499016,250.5443412769742,0.5596964608782755,250.17846968398487,0.6757459004993481,248.0373475773503 +6.5,1.0538516754549905,263.24979308361293,0.5911280054788979,270.17526440026194,0.6261735818734251,269.4331709994467,0.7570510876024378,267.6194094944405 +6.75,1.179681054388145,283.8891999378015,0.6579446297589594,289.6952735409024,0.698173190674643,288.6036862383063,0.8452662916072492,287.3960532390034 +7.0,1.3155962481407875,305.30745233365764,0.7260040195544638,309.5654724033882,0.7739094140858843,308.22703538997825,0.9385032272603785,307.69998670900804 +7.25,1.461394783774867,327.50455027118113,0.8003673193737004,329.85873305516077,0.8548361282890699,328.25775468576677,1.038866077145879,328.54034600403617 +7.5,1.6173816158025585,350.4804937503723,0.8763512035493293,350.60979092143623,0.9423986740631011,350.48049375037226,1.145377164471663,350.52359280739364 +7.75,1.784690037695071,374.23528277123074,0.9607242094572293,374.2352827712308,1.0321323375021068,374.2352827712308,1.2591821948848023,374.23528277123086 +8.0,1.9636062319651761,398.7689173337569,1.0478082087280762,398.76891733375686,1.129105942561379,398.7689173337569,1.380173461084877,398.7689173337569 +8.25,2.1537590950872496,424.08139743795044,1.1404823799115325,424.0813974379505,1.231487932044032,424.0813974379505,1.5085764690142716,424.0813974379505 +8.5,2.3558672754952124,450.1727230838116,1.2390979142452607,450.17272308381155,1.339373842097445,450.17272308381155,1.644779677279306,450.1727230838116 +8.75,2.570466947210423,477.04289427134006,1.3424834224693887,477.04289427134,1.4533111179762472,477.04289427134,1.7887538292186862,477.04289427134 +9.0,2.797276397680262,504.691911000536,1.4513050982342601,504.691911000536,1.5732821301415056,504.6919110005359,1.940621208685343,504.691911000536 +9.25,3.035104037397528,533.1197732713995,1.5655804723654327,533.1197732713995,1.7002003605595846,533.1197732713996,2.1002949567741815,533.1197732713995 +9.5,3.2896571931228045,562.3264810839306,1.686110338733192,562.3264810839306,1.8339219687384225,562.3264810839305,2.26989650019814,562.3264810839306 +9.75,3.538456686682414,612.5166944582794,1.7858461172904017,592.3120344381291,1.9851679268296722,592.3120344381292,2.4364902436008298,599.0469211115125 +10.0,3.622327044157654,593.9648185459672,1.9889544342509302,623.0764333339952,2.1069179639718905,623.0764333339952,2.5727331474601582,613.3725617379858 +10.25,3.6223375320497415,525.2025663074712,2.30993145503877,654.6196777715287,2.2140034190286024,654.6196777715286,2.7154241353723716,611.4806406168428 +10.5,3.62234527734387,464.6979372598775,2.666073800941529,686.9417677507297,2.337803551415242,686.9417677507297,2.8754075432335466,612.8604909204456 +10.75,3.6223456681827697,445.0439593608759,2.9606100870746572,720.0427032715979,2.482742209159499,720.0427032715983,3.021899321472309,628.3764553013574 +11.0,3.622345769722794,426.4370647559222,3.27624078018154,753.9224843341342,2.6353575855292637,753.9224843341343,3.1779813784778654,644.7606778080635 +11.25,3.6223454307318868,411.80616837254763,3.578423794301074,822.4510624064279,2.7480868985464744,788.5811109383376,3.3162853745264784,674.2794472391043 +11.5,3.622345212202052,397.73061588587467,3.622330984627368,738.776816790086,3.108656156043992,824.0185830842086,3.451110784291137,653.5086719200564 +11.75,3.6223451778428677,386.2920271761649,3.6223434032677266,617.4630233255486,3.5992860615206186,901.1513918181668,3.6146582142104045,634.9688141066267 +12.0,3.6223466526568617,374.469698853721,3.622345570382401,567.670188946981,3.622334424922358,759.9805070541832,3.6223422159872065,567.3734649516285 +12.25,3.6223521426624448,365.3930081133308,3.622345801110452,532.6196659296982,3.6223452993031855,629.4379184061999,3.6223477476920274,509.1501974830762 +12.5,3.622357264503283,354.76910800661454,3.622345402349194,504.8639535119541,3.622345807703222,581.4420636125067,3.622349491518566,480.35837504369175 +12.75,3.622352040682574,347.52962114046784,3.6223452046832794,481.4415306421037,3.622345557452205,545.3797914285076,3.622347600939353,458.1169810703597 +13.0,3.622346816861865,338.86542794638535,3.622345165112278,459.73757800659337,3.6223452160904928,512.2592853697919,3.622345732688212,436.95409710759014 +13.25,3.6223456739346482,331.7746425985911,3.622349410339327,441.9485871031459,3.622345173969946,487.070039423884,3.622346752747974,420.26442304187367 +13.5,3.6223464168705184,324.85396105814306,3.622355607086224,424.7286808132697,3.622348499943041,463.7863554345775,3.622350174633261,404.45633243533007 +13.75,3.6223476508849335,317.5684183562588,3.622353278523889,410.12694085292713,3.6223550866282928,443.5721451851565,3.622352005345705,390.42250146478074 +14.0,3.6223505559443914,312.1175813992215,3.622347486961833,396.340111042303,3.622353464143539,426.10042488817834,3.622350502349921,378.18603910990095 +14.25,3.6223534610038493,305.64410182658366,3.6223456675677923,383.93037026125967,3.6223473058777307,409.6873954059334,3.6223488114831244,366.42062249792554 +14.5,3.622361209304207,300.34190314783376,3.6223464905927485,372.52394206023916,3.6223457288629564,395.65820351555556,3.6223511429199706,356.1746829078761 +14.75,3.6223713497021293,295.244499680246,3.6223482435252095,361.5668156972668,3.622346592583677,382.27506612476355,3.622355395270338,346.3621271674254 +15.0,3.622381107897377,289.3094836819298,3.6223514746748067,352.0888226404568,3.622348767421511,370.38632976167185,3.6223604499978985,337.26154536135283 +15.25,3.6223679350932647,285.352679875891,3.622355153010334,342.13846598442024,3.6223521073292693,359.5032803609661,3.6223583984776226,328.9981420737591 +15.5,3.6223547622891514,280.65214750492476,3.6223662387187296,334.6838236333273,3.622357768974054,349.0201527333878,3.6223595899939784,321.4520412905466 +15.75,3.6223452178468003,275.8554713842255,3.622377315436947,326.20843001549594,3.6223693100521546,340.1993919560072,3.6223639477786342,314.08776445190955 +16.0,3.6223452852955873,272.14451990805355,3.6223723306154363,319.13543863945705,3.622380662773501,330.60913064794033,3.6223660928948416,307.29636306515033 +16.25,3.6223453527443747,267.78284083947005,3.6223580544193097,312.351326167589,3.6223674886124706,323.47646930823754,3.6223569652587186,301.2035454384322 +16.5,3.6223455997679705,263.77439099358617,3.6223452062593764,305.04360546855776,3.6223527565286413,315.59489208047097,3.6223478541853296,294.80429618087163 +16.75,3.622346067683757,260.4368495245578,3.622345278369686,299.6669940502601,3.622345234579533,308.5879685738418,3.622345526877659,289.56393738288654 +17.0,3.6223465355995437,256.53526623333073,3.622345350809193,293.4777751736312,3.6223453089944946,302.2524902508098,3.622345731801077,284.0885105525906 +17.25,3.6223469033165125,253.15468514952659,3.6223456176378974,287.83747183039765,3.6223453828734202,295.2449858786952,3.6223459679426093,278.74571428620646 +17.5,3.6223471876465987,250.18625027958421,3.6223461141635944,283.04656559992424,3.622345848437908,290.1078555478625,3.6223463834160334,274.4468904757903 +17.75,3.6223474719766844,246.72796513692094,3.6223466126259307,277.56427398262866,3.6223463589144496,284.5003085072003,3.622346814505688,269.59751587558327 +18.0,3.6223475927245348,243.77018586700058,3.6223469668236885,273.1987281519374,3.62234681808986,278.93185407131205,3.622347125879361,265.30025603008335 +18.25,3.6223475913750427,241.12290574356504,3.62234726661626,268.96118812935936,3.6223471237128897,274.6236451878622,3.622347327234731,261.5692463535956 +18.5,3.62234759002555,238.0474518670645,3.6223475674121866,264.13282728566924,3.6223474308101764,269.69490765881375,3.6223475294159706,257.29172893718254 +18.75,3.6223505559090903,235.36562226359507,3.6223475922003496,260.6265803910526,3.6223475928313693,265.298574042765,3.6223485803136035,253.76359223247087 +19.0,3.6223560958826355,233.02054289134293,3.622347590785588,256.8670493997568,3.6223475913922174,261.4790253808091,3.6223504260201467,250.45553922396962 +19.25,3.6223616358561803,230.30018410562772,3.622347720436771,252.63735523204332,3.622347589947677,257.1321081166413,3.622352315413543,246.6898824847708 +19.5,3.622367111113461,227.80474241221805,3.622353527484353,249.67675386674856,3.6223511924067098,253.52303435481736,3.6223572770015084,243.66817687792798 +19.75,3.622372504381353,225.73917370416368,3.622359305873421,246.35612859285837,3.6223570588508642,250.1528788863831,3.6223629563685464,240.7493937278017 +20.0,3.622377897649244,223.34370979340312,3.6223650888797794,242.691400965745,3.622362940768238,246.33418238765256,3.62236864243242,237.45643104893352 +20.25,3.622378751344619,220.97447576191917,3.6223707226850097,240.07264724085917,3.622368657787592,243.2701452469404,3.6223727106057404,234.77242274990624 +20.5,3.6223688313095384,219.16290871351248,3.6223763227349686,237.14983113249136,3.62237433308903,240.3073450100827,3.6223731623778455,232.20669495202887 +20.75,3.6223589112744583,217.06036372639136,3.622381245728034,233.90249927073725,3.6223800223773535,236.95393693707075,3.622373393126615,229.30559997806645 +21.0,3.6223498991938663,214.77725521343012,3.622370912355717,231.59149207536626,3.6223743866347027,234.27491362657628,3.6223650660614286,226.88122030512423 +21.25,3.6223483547555735,213.18048352650325,3.622360651570205,229.02347976311057,3.6223640053420145,231.67586744961363,3.6223576705559313,224.62661024640917 +21.5,3.62234681031728,211.32550233520902,3.6223503716747714,226.12151537994174,3.622353597733897,228.73328500592248,3.6223502599086497,222.06010090702443 +21.75,3.622345265878987,209.20671361479006,3.622348510651997,224.04336285618683,3.622349007565036,226.31548618780528,3.6223475946986734,219.8551875529274 +22.0,3.6223452945520713,207.72861440003223,3.622346917495259,221.770286765771,3.6223473987001054,224.02145701238015,3.622346536915812,217.84011939272781 +22.25,3.622345529483466,206.1391689995733,3.6223453234911682,219.21165417367249,3.6223457867244733,221.42828206062507,3.6223455465663696,215.59303507795696 +22.5,3.62234576441486,204.31694654232373,3.6223452937209326,217.30104897217788,3.622345224237542,219.2177498470157,3.6223454274577787,213.61191512050576 +22.75,3.622346530226876,202.82096730012978,3.6223455357517635,215.33117335057696,3.622345467947866,217.24334010530134,3.622345844642169,211.79849358533602 +23.0,3.6223476069689164,201.4450497591952,3.622345777413467,213.12391674685207,3.6223457119540234,215.01076199749204,3.6223463654454693,209.85990950117977 +23.25,3.6223486837109564,199.86203641422932,3.6223466222140046,211.32274952021768,3.622346327483907,212.96127587895504,3.6223472111362898,208.048687271134 +23.5,3.6223489299758347,198.3161869344144,3.6223477305136,209.59704504938503,3.622347440643897,211.24339282126854,3.6223480337111105,206.3855416016893 +23.75,3.622347669332914,197.0259509164487,3.6223488358576597,207.66806742106073,3.622348554649284,209.2993849713005,3.622348353279953,204.6644677696033 +24.0,3.6223464086899932,195.54644353926327,3.6223487167236206,205.88931405716957,3.6223490392476294,207.3346240507901,3.622348054887081,202.92346054907432 +24.25,3.6223474633081834,193.92795698477354,3.6223474199246186,204.2680387216069,3.62234773870111,205.72647092231404,3.6223475406446375,201.30748887623153 +24.5,3.6223816383744594,192.90313794943168,3.6223461302637183,202.4788585727089,3.6223464381759505,203.91678141799719,3.622358068938043,199.76625931337927 +24.75,3.6224158134407345,191.71261504352606,3.622355835904942,200.7463592922794,3.622347723152897,201.9682922403442,3.622373124166191,198.1424221920499 +25.0,3.6224499885070096,190.3532371147166,3.622390873835434,199.43780060723174,3.6223828363394532,200.67297861998046,3.6224078995606326,196.82133878064292 +25.25,0.0,0.5183289236824425,0.0,0.5183289236824425,0.0,0.5183289236824425,0.0,0.5183289236824425 +25.5,0.0,0.5286436743448812,0.0,0.5286436743448812,0.0,0.5286436743448812,0.0,0.5286436743448812 +25.75,0.0,0.5390600481665556,0.0,0.5390600481665556,0.0,0.5390600481665556,0.0,0.5390600481665556 +26.0,0.0,0.5495780451474659,0.0,0.5495780451474659,0.0,0.5495780451474659,0.0,0.5495780451474659 +26.25,0.0,0.560197665287612,0.0,0.560197665287612,0.0,0.560197665287612,0.0,0.560197665287612 +26.5,0.0,0.5709189085869939,0.0,0.5709189085869939,0.0,0.5709189085869939,0.0,0.5709189085869939 +26.75,0.0,0.5817417750456116,0.0,0.5817417750456116,0.0,0.5817417750456116,0.0,0.5817417750456116 +27.0,0.0,0.5926662646634653,0.0,0.5926662646634653,0.0,0.5926662646634653,0.0,0.5926662646634653 +27.25,0.0,0.6036923774405549,0.0,0.6036923774405549,0.0,0.6036923774405549,0.0,0.6036923774405549 +27.5,0.0,0.6148201133768803,0.0,0.6148201133768803,0.0,0.6148201133768803,0.0,0.6148201133768803 +27.75,0.0,0.6260494724724414,0.0,0.6260494724724414,0.0,0.6260494724724414,0.0,0.6260494724724414 +28.0,0.0,0.6373804547272386,0.0,0.6373804547272386,0.0,0.6373804547272386,0.0,0.6373804547272386 +28.25,0.0,0.6488130601412714,0.0,0.6488130601412714,0.0,0.6488130601412714,0.0,0.6488130601412714 +28.5,0.0,0.6603472887145402,0.0,0.6603472887145402,0.0,0.6603472887145402,0.0,0.6603472887145402 +28.75,0.0,0.6719831404470448,0.0,0.6719831404470448,0.0,0.6719831404470448,0.0,0.6719831404470448 +29.0,0.0,0.6837206153387851,0.0,0.6837206153387851,0.0,0.6837206153387851,0.0,0.683720615338785 +29.25,0.0,0.6955597133897615,0.0,0.6955597133897615,0.0,0.6955597133897615,0.0,0.6955597133897614 +29.5,0.0,0.7075004345999736,0.0,0.7075004345999736,0.0,0.7075004345999736,0.0,0.7075004345999737 +29.75,0.0,0.7195427789694215,0.0,0.7195427789694215,0.0,0.7195427789694215,0.0,0.7195427789694215 +30.0,0.0,0.7316867464981054,0.0,0.7316867464981054,0.0,0.7316867464981054,0.0,0.7316867464981054 +30.25,0.0,0.7439323371860251,0.0,0.7439323371860251,0.0,0.7439323371860251,0.0,0.7439323371860252 diff --git a/demo/documented/studies/powercurve_validation/writeup/floris_run.py b/demo/documented/studies/powercurve_validation/writeup/floris_run.py new file mode 100644 index 00000000..236f3916 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/floris_run.py @@ -0,0 +1,66 @@ + +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd + +from floris import FlorisModel, TimeSeries + +plt.style.use([ + "dark_background", +]) + +fmodel = FlorisModel("FLORIS_case.yaml") + +air_density = 1.225 # kg/m^3 +wind_speeds = np.arange(0, 30*1.01, 0.25)[1:] +wind_directions = 270.0*np.ones_like(wind_speeds) +wind_data = TimeSeries( + wind_directions=wind_directions, + wind_speeds=wind_speeds, + turbulence_intensities=0.06, +) +wind_data.assign_ti_using_IEC_method() +fmodel.set( + air_density=air_density, + wind_data=wind_data, +) + +diameter_rotor = 130.0 # m +area_rotor = np.pi/4*diameter_rotor**2 # m^2 +fluid_force_fun = lambda V : 0.5*air_density*area_rotor*V**2 +fluid_power_fun = lambda V : 0.5*air_density*area_rotor*V**3 + +fmodel.run() + +P_T = fmodel.get_turbine_powers()/1e6 # MW +P_mean = fmodel.get_farm_power()/1e6/P_T.shape[1] # MW + +CTf_T = fmodel.get_turbine_thrust_coefficients() # kN +Tf_T = CTf_T*np.expand_dims(fluid_force_fun(wind_speeds), -1)/1e3 # kN + +fig, ax = plt.subplots() +ax.plot(wind_speeds, P_T, "-") +ax.plot(wind_speeds, P_mean, "--") + +fig, ax = plt.subplots() +# ax.plot(wind_speeds, P_mean) +ax.plot(wind_speeds, Tf_T, "--") + +print(wind_directions) +print(wind_speeds) +print(P_mean) + +plt.show() + +pd.DataFrame({ + "V": wind_speeds.tolist(), + "P1": P_T[:,0].tolist(), + "T1": Tf_T[:,0].tolist(), + "P2": P_T[:,1].tolist(), + "T2": Tf_T[:,1].tolist(), + "P3": P_T[:,2].tolist(), + "T3": Tf_T[:,2].tolist(), + "Ptot": P_mean.tolist(), + "Ttot": np.mean(Tf_T,axis=1).tolist(), +}).to_csv("floris.csv",index=False) + diff --git a/demo/documented/studies/powercurve_validation/writeup/plot.py b/demo/documented/studies/powercurve_validation/writeup/plot.py new file mode 100644 index 00000000..888fd304 --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/plot.py @@ -0,0 +1,274 @@ +import pprint +import yaml + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns + +plt.style.use( + [ + "dark_background", + "https://raw.githubusercontent.com/cfrontin/tools_cvf/main/tools_cvf/stylesheet_cvf.mplstyle", + "https://raw.githubusercontent.com/cfrontin/tools_cvf/main/tools_cvf/stylesheet_nrel.mplstyle", + ] +) + +with open("FLORIS_IEA-3p4-130-RWT.yaml", "r") as f_floris: + data_FLORIS = yaml.safe_load(f_floris) + +do_single_turbine=False +do_multi_turbine=True + +### SINGLE TURBINE CASE + +if do_single_turbine: + df_curves_collection = [] + label_collection = [] + if True: + # extract three levels of discretization + casename = "kestrel_nx19ny13nz06" # "nx19ny13nz06" + df_curves_lo = pd.read_csv(f"../curves_{casename}.csv", names=["V","P","Tf"]) + df_curves_collection.append(df_curves_lo) + label_collection.append("$N_x=19$, $N_y=13$, $N_z=6$") + if True: + casename = "kestrel_nx27ny18nz09" # "nx27ny18nz09" + df_curves_mid = pd.read_csv(f"../curves_{casename}.csv", names=["V","P","Tf"]) + df_curves_collection.append(df_curves_mid) + label_collection.append("$N_x=27$, $N_y=18$, $N_z=9$") + if True: + casename = "kestrel_nx38ny25nz13" # "nx38ny25nz13" + df_curves_hi = pd.read_csv(f"../curves_{casename}.csv", names=["V","P","Tf"]) + df_curves_collection = [df_curves_lo, df_curves_mid, df_curves_hi] + label_collection.append("$N_x=38$, $N_y=25$, $N_z=13$") + if True: + casename = "kestrel_nx54ny36nz18" # "nx54ny36nz18" + df_curves_hi = pd.read_csv(f"../curves_{casename}.csv", names=["V","P","Tf"]) + df_curves_collection = [df_curves_lo, df_curves_mid, df_curves_hi] + label_collection.append("$N_x=54$, $N_y=36$, $N_z=18$") + + # get FLORIS reference data + rho_fluid = 1.225 + D_rotor=130.0 + A_rotor = np.pi/4*D_rotor**2 + Ts_fluid = lambda V: 0.5*rho_fluid*A_rotor*V**2/1e3 # kN + Ps_fluid = lambda V: 0.5*rho_fluid*A_rotor*V**3/1e6 # MW + + maxV = np.max([np.max(df.V) for df in df_curves_collection]) + V_floris = np.array(data_FLORIS["power_thrust_table"]["wind_speed"])[data_FLORIS["power_thrust_table"]["wind_speed"] <= maxV] + P_floris = np.array(data_FLORIS["power_thrust_table"]["power"])[data_FLORIS["power_thrust_table"]["wind_speed"] <= maxV] + Vrated_floris = 9.812675420388173 + Prated_floris = 3622.3451711828 + CPrated_floris = Prated_floris/1e3/Ps_fluid(Vrated_floris) + CT_floris = np.array(data_FLORIS["power_thrust_table"]["thrust_coefficient"])[data_FLORIS["power_thrust_table"]["wind_speed"] <= maxV] + + ## power plot + + fig, ax = plt.subplots() + + for df_power, label_case in zip(df_curves_collection, label_collection): + ax.plot(df_power.V, rho_fluid*df_power.P, label=label_case) + # ax.plot(df_power_lo.V, rho_fluid*df_power_lo.P, label="$N_x=19$, $N_y=13$, $N_z=6$") + # ax.plot(df_power_mid.V, rho_fluid*df_power_mid.P, label="$N_x=27$, $N_y=18$, $N_z=9$") + # ax.plot(df_power_hi.V, rho_fluid*df_power_hi.P, label="$N_x=38$, $N_y=25$, $N_z=13$") + ax.plot( + V_floris, + P_floris/1e3, + "w--", + label="IEA-3.4MW-130m reference", + ) + ax.plot(Vrated_floris, Prated_floris/1e3, "wo") + ax.plot( + V_floris, + np.minimum( + Prated_floris/1e3*np.ones_like(V_floris), + CPrated_floris*Ps_fluid(V_floris)), + "w:", + label="IEA-3.4MW-130m: $\\min(\\frac{1}{2} C_P A V^3, P_{\\mathrm{rated}}/\\rho)$" + ) + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("turbine power, $P(V)$ (MW)") + ax.legend() + fig.savefig("single_turb_Pconv.png", dpi=300, bbox_inches="tight") + + ## CP plot + + fig, ax = plt.subplots() + + for df_power, label_case in zip(df_curves_collection, label_collection): + ax.plot(df_power.V, rho_fluid*df_power.P/Ps_fluid(df_power.V), label=label_case) + # ax.plot(df_power_lo.V, rho_fluid*df_power_lo.P/Ps_fluid(df_power_lo.V), label="$N_x=19$, $N_y=13$, $N_z=6$") + # ax.plot(df_power_mid.V, rho_fluid*df_power_mid.P/Ps_fluid(df_power_mid.V), label="$N_x=27$, $N_y=18$, $N_z=9$") + # ax.plot(df_power_hi.V, rho_fluid*df_power_hi.P/Ps_fluid(df_power_hi.V), label="$N_x=38$, $N_y=25$, $N_z=13$") + ax.plot( + V_floris[1:], + P_floris[1:]/1e3/Ps_fluid(V_floris[1:]), + "w--", + label="IEA-3.4MW-130m reference", + ) + ax.plot(Vrated_floris, CPrated_floris, "wo") + ax.plot( + V_floris[1:], + np.minimum( + Prated_floris/1e3/Ps_fluid(V_floris[1:]), + CPrated_floris), + "w:", + label="IEA-3.4MW-130m: $\\min(C_P, P_{\\mathrm{rated}}/P_{\\mathrm{fluid}})$" + ) + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("power_coefficient, $C_P(V)$ (-)") + ax.legend() + fig.savefig("single_turb_CPconv.png", dpi=300, bbox_inches="tight") + + # fig, ax = plt.subplots() + # + # for df_power, label_case in zip(df_curves_collection, label_collection): + # ax.plot( + # df_power.V, + # ( + # rho_fluid*df_power.P/Ps_fluid(df_power.V) + # )/np.interp( + # df_power.V, + # V_floris[1:], + # P_floris[1:]/1e3/Ps_fluid(V_floris[1:]), + # ), + # label=label_case, + # ) + + ## thrust plot + + fig, ax = plt.subplots() + + for df_thrust, label_case in zip(df_curves_collection, label_collection): + ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf, label=label_case) + # ax.plot(df_thrust_lo.V, rho_fluid*df_thrust_lo.Tf, label="$N_x=19$, $N_y=13$, $N_z=6$") + # ax.plot(df_thrust_mid.V, rho_fluid*df_thrust_mid.Tf, label="$N_x=27$, $N_y=18$, $N_z=9$") + # ax.plot(df_thrust_hi.V, rho_fluid*df_thrust_hi.Tf, label="$N_x=38$, $N_y=25$, $N_z=13$") + ax.plot( + V_floris, + CT_floris*Ts_fluid(V_floris), + "w--", + label="IEA-3.4MW-130m reference", + ) + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("turbine thrust, $T(V)$ (kN)") + ax.legend() + fig.savefig("single_turb_Tconv.png", dpi=300, bbox_inches="tight") + + ## CT plot + + fig, ax = plt.subplots() + + for df_thrust, label_case in zip(df_curves_collection, label_collection): + ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf/Ts_fluid(df_thrust.V), label=label_case) + # ax.plot(df_thrust_lo.V, rho_fluid*df_thrust_lo.Tf/Ts_fluid(df_thrust_lo.V), label="$N_x=19$, $N_y=13$, $N_z=6$") + # ax.plot(df_thrust_mid.V, rho_fluid*df_thrust_mid.Tf/Ts_fluid(df_thrust_mid.V), label="$N_x=27$, $N_y=18$, $N_z=9$") + # ax.plot(df_thrust_hi.V, rho_fluid*df_thrust_hi.Tf/Ts_fluid(df_thrust_hi.V), label="$N_x=38$, $N_y=25$, $N_z=13$") + ax.plot( + V_floris, + CT_floris, + "w--", + label="IEA-3.4MW-130m reference", + ) + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("thrust coefficient, $C_T(V)$ (-)") + ax.legend() + fig.savefig("single_turb_CTconv.png", dpi=300, bbox_inches="tight") + + + +### MULTI TURBINE CASE + +if do_multi_turbine: + df_curves_collection = [] + label_collection = [] + if True: + # extract three levels of discretization + casename = "kestrel_nx64ny13nz06" # "nx64ny13nz06" + df_curves_lo = pd.read_csv(f"../curves_multi_{casename}.csv", names=["V","P1","Tf1","P2","Tf2","P3","Tf3","P_tot","Tf_tot"]) + df_curves_collection.append(df_curves_lo) + label_collection.append("$N_x=64$, $N_y=13$, $N_z=6$") + if True: + # extract three levels of discretization + casename = "kestrel_nx90ny18nz09" # "nx90ny18nz09" + df_curves_lo = pd.read_csv(f"../curves_multi_{casename}.csv", names=["V","P1","Tf1","P2","Tf2","P3","Tf3","P_tot","Tf_tot"]) + df_curves_collection.append(df_curves_lo) + label_collection.append("$N_x=90$, $N_y=18$, $N_z=9$") + if True: + # extract three levels of discretization + casename = "kestrel_nx127ny25nz13" # "nx127ny25nz13" + df_curves_lo = pd.read_csv(f"../curves_multi_{casename}.csv", names=["V","P1","Tf1","P2","Tf2","P3","Tf3","P_tot","Tf_tot"]) + df_curves_collection.append(df_curves_lo) + label_collection.append("$N_x=127$, $N_y=25$, $N_z=13$") + if False: + # extract three levels of discretization + casename = "kestrel_nx180ny36nz18" # "nx180ny36nz18" + df_curves_lo = pd.read_csv(f"../curves_multi_{casename}.csv", names=["V","P1","Tf1","P2","Tf2","P3","Tf3","P_tot","Tf_tot"]) + df_curves_collection.append(df_curves_lo) + label_collection.append("$N_x=180$, $N_y=36$, $N_z=18$") + + # get FLORIS reference data + rho_fluid = 1.225 + D_rotor=130.0 + A_rotor = np.pi/4*D_rotor**2 + Ts_fluid = lambda V: 0.5*rho_fluid*A_rotor*V**2/1e3 # kN + Ps_fluid = lambda V: 0.5*rho_fluid*A_rotor*V**3/1e6 # MW + + df_floris_multi = pd.read_csv("floris.csv") + V_floris = df_floris_multi.V + P1_floris = df_floris_multi.P1 + P2_floris = df_floris_multi.P2 + P3_floris = df_floris_multi.P3 + Tf1_floris = df_floris_multi.T1 + Tf2_floris = df_floris_multi.T2 + Tf3_floris = df_floris_multi.T3 + + ## power plot + + fig, ax = plt.subplots() + + ax.plot([], [], "w--", label="T1") + ax.plot([], [], "w:", label="T2") + ax.plot([], [], "w-.", label="T3") + + for df_power, label_case in zip(df_curves_collection, label_collection): + pt0 = ax.plot([], [], "-", label=label_case) + ax.plot(df_power.V, rho_fluid*df_power.P1, "--", c=pt0[-1].get_color(), label="__") + ax.plot(df_power.V, rho_fluid*df_power.P2, ":", c=pt0[-1].get_color(), label="__") + ax.plot(df_power.V, rho_fluid*df_power.P3, "-.", c=pt0[-1].get_color(), label="__") + # ax.plot(df_power.V, rho_fluid*df_power.P_tot/3, "-", c=pt0[-1].get_color(), label=label_case) + ax.plot([], [], "w-", label="FLORIS") + ax.plot(V_floris, P1_floris, "w--", label="_FLORIS T1_") + ax.plot(V_floris, P2_floris, "w:", label="_FLORIS T2_") + ax.plot(V_floris, P3_floris, "w-.", label="_FLORIS T3_") + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("turbine power, $P(V)$ (MW)") + ax.legend() + fig.savefig("multi_turb_Pconv.png", dpi=300, bbox_inches="tight") + + ## thrust plot + + fig, ax = plt.subplots() + + ax.plot([], [], "w--", label="T1") + ax.plot([], [], "w:", label="T2") + ax.plot([], [], "w-.", label="T3") + + for df_thrust, label_case in zip(df_curves_collection, label_collection): + pt0 = ax.plot([], [], "-", label=label_case) + ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf1, "--", c=pt0[-1].get_color(), label="__") + ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf2, ":", c=pt0[-1].get_color(), label="__") + ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf3, "-.", c=pt0[-1].get_color(), label="__") + # ax.plot(df_thrust.V, rho_fluid*df_thrust.Tf_tot/3, "-", c=pt0[-1].get_color(), label=label_case) + ax.plot([], [], "w-", label="FLORIS") + ax.plot(V_floris, Tf1_floris, "w--", label="_FLORIS T1_") + ax.plot(V_floris, Tf2_floris, "w:", label="_FLORIS T2_") + ax.plot(V_floris, Tf3_floris, "w-.", label="_FLORIS T3_") + ax.set_xlabel("hub-height farfield velocity, $V$ (m/s)") + ax.set_ylabel("turbine thrust, $T(V)$ (kN)") + ax.legend() + fig.savefig("multi_turb_Tconv.png", dpi=300, bbox_inches="tight") + + + +plt.show() diff --git a/demo/documented/studies/powercurve_validation/writeup/power_curve_validation.md b/demo/documented/studies/powercurve_validation/writeup/power_curve_validation.md new file mode 100644 index 00000000..b129b90e --- /dev/null +++ b/demo/documented/studies/powercurve_validation/writeup/power_curve_validation.md @@ -0,0 +1,50 @@ + +# WindSE Power Curve Validation + +`WindSE` has the ability to feather thrust and power coefficients in order to emulate a controlled turbine through multiple regions of the power curve. +In this study, we will validate the power curve behavior in isolation and with multiple inline turbines. +We consider the discretization dependence of the `WindSE` results, and the comparison of `WindSE` power estimates compared to other tools. + +## Tuning methodology + +Using the `IEA-3.4MW-130m` turbine, we run a calibration study at one region II velocity and one region III velocity to adjust the power curve model parameters to match the reference curves for the turbine. +A magic number scaling parameter is applied to the output power to achieve a match between the force/velocity product power and the expected output shaft power. + +## Single turbine: discretization dependence + +In this section we compare the `WindSE` results to the canonical IEA 3.4MW 130m turbine's power curve that is an input to `FLORIS`. + +![Power curve convergence plot](single_turb_Tconv.png) + +The first plot here shows the thrust curves compared to the reference turbine data, and demonstrates a strong match for a single turbine. +The thrust peaking behavior in the reference power curve is due to the lack of thrust peak-shaving control, and `WindSE` implicitly performs peak shaving by the smoothing method used to achieve differentiability. +Long run behavior in velocity exhibits some inexact scaling; this is likely due to inaccuracy-- or at least mismatch-- in the induction profiles assumed for `WindSE`'s actuator disk. + +![Power curve convergence plot](single_turb_Pconv.png) + +FLORIS includes cut-in behavior that is not modeled in region I or region I.5. +Comparison with the standard IEA-3.4MW-130m power curve used for FLORIS demonstrates similar power behavior in region II, with an slight apparent overprediction of $C_P$. +The FLORIS curve does not have thrust-peak shaving, while the `WindSE` power curve implicitly does (via the smoothing), so there is a power mismatch in region II.5. + +![$C_P$ curve convergence plot](single_turb_CPconv.png) + +Plotting the power coefficient indicates that $C_P$ is in fact initially overpredicted in region II, and moreover $C_P$ is not quite constant in region II. +As noted previously, region II.5 is not modeled in the FLORIS power curve, but long run behavior of $C_P$ appears to be accurate, in region III. + +## Multi-turbine: discretization dependence + +We now run a three-turbine case with waked turbines to validate that behavior matches expectation. +In this case, we have three directly-waked turbines at $7.5D$ spacing. + +![Farm thrust curve convergence plot](multi_turb_Tconv.png) + +As previously, FLORIS does not implement peak-shaving which is implicitly included in the WindSE results. +Outside of region II.5 (peak-shaving), there is a strong qualitative match between the thrust of the turbines in WindSE and the FLORIS simulation. +This indicates that WindSE is correctly estimating the thrust on a given turbine as well as passing the resulting thrust into the flowfield correctly. + +![Farm power curve convergence plot](multi_turb_Pconv.png) + +In the power curves, WindSE again makes a strong prediction that is closely matched with the FLORIS result, less thrust peak-shaving in region II.5. + + + diff --git a/windse/DomainManager.py b/windse/DomainManager.py index 534f5217..7edd6e33 100644 --- a/windse/DomainManager.py +++ b/windse/DomainManager.py @@ -1377,7 +1377,7 @@ def RecomputeBoundaryMarkers(self,inflow_angle): diagonal_ids = [outflow_id,outflow_id,inflow_id,inflow_id] ### Count the number of pi/2 sections in the new inflow_angle ### - turns = inflow_angle/(pi/2) + turns = ((inflow_angle % (2*np.pi) + (2*np.pi)) % (2*np.pi))/(pi/2) ### New order ### tol = 1e-3 diff --git a/windse/OptimizationManager.py b/windse/OptimizationManager.py index 8a5a8842..a82f7b1d 100644 --- a/windse/OptimizationManager.py +++ b/windse/OptimizationManager.py @@ -2,7 +2,7 @@ The OptimizationManager submodule contains all the required function for optimizing via dolfin-adjoint. To use dolfin-adjoin set:: - general: + general: dolfin_adjoint: True in the param.yaml file. @@ -21,7 +21,7 @@ main_file = os.path.basename(__main__.__file__) else: main_file = "ipython" - + ### This checks if we are just doing documentation ### if not main_file in ["sphinx-build", "__main__.py"]: from dolfin import * @@ -47,46 +47,46 @@ import matplotlib.pyplot as plt else: InequalityConstraint = object - + import openmdao.api as om class ObjComp(om.ExplicitComponent): """ OpenMDAO component to wrap the objective computation from dolfin. - + Specifically, we use the J and dJ (function and Jacobian) methods to compute the function value and derivative values as needed by the - OpenMDAO optimizers. + OpenMDAO optimizers. """ def initialize(self): self.options.declare('initial_DVs', types=np.ndarray) self.options.declare('J', types=object) self.options.declare('dJ', types=object) self.options.declare('callback', types=object) - + def setup(self): self.add_input('DVs', val=self.options['initial_DVs']) self.add_output('obj', val=0.) self.declare_partials('*', '*') - + def compute(self, inputs, outputs): m = list(inputs['DVs']) computed_output = self.options['J'](m) outputs['obj'] = computed_output self.options['callback'](m) - + def compute_partials(self, inputs, partials): m = list(inputs['DVs']) jac = self.options['dJ'](m) partials['obj', 'DVs'] = jac - - + + class ConsComp(om.ExplicitComponent): """ OpenMDAO component to wrap the constraint computation. - + A small wrapper used on the fenics methods for computing constraint and Jacobian values using the OpenMDAO syntax. """ @@ -95,26 +95,26 @@ def initialize(self): self.options.declare('J', types=object) self.options.declare('dJ', types=object) self.options.declare('con_name', types=str) - + def setup(self): self.con_name = self.options["con_name"] self.add_input('DVs', val=self.options['initial_DVs']) - + output = self.options['J'](self.options['initial_DVs']) self.add_output(self.con_name, val=output) self.declare_partials('*', '*') - + def compute(self, inputs, outputs): m = list(inputs['DVs']) computed_output = self.options['J'](m) outputs[self.con_name] = computed_output - + def compute_partials(self, inputs, partials): m = list(inputs['DVs']) jac = self.options['dJ'](m) partials[self.con_name, 'DVs'] = jac - + def gather(m): """ @@ -126,13 +126,13 @@ def gather(m): return m._ad_to_list(m) else: return m # Assume it is gathered already - + def om_wrapper(J, initial_DVs, dJ, H, bounds, **kwargs): """ Custom optimization wrapper to use OpenMDAO optimizers with dolfin-adjoint. - + Follows the API as defined by dolfin-adjoint. - + Parameters ---------- J : object @@ -146,30 +146,30 @@ def om_wrapper(J, initial_DVs, dJ, H, bounds, **kwargs): Function to compute the Hessian at a design point (not used). bounds : array Array of lower and upper bound values for the design variables. - + Returns ------- DVs : array The optimal design variable values. """ - + # build the model prob = om.Problem(model=om.Group()) - + if 'callback' in kwargs: callback = kwargs['callback'] else: callback = None prob.model.add_subsystem('obj_comp', ObjComp(initial_DVs=initial_DVs, J=J, dJ=dJ, callback=callback), promotes=['*']) - + constraint_types = [] if 'constraints' in kwargs: constraints = kwargs['constraints'] - + if not isinstance(constraints, list): constraints = [constraints] - + for idx, c in enumerate(constraints): if isinstance(c, InequalityConstraint): typestr = "ineq" @@ -183,13 +183,13 @@ def jac(x): return [gather(y) for y in out] constraint_types.append(typestr) - + con_name = f'con_{idx}' prob.model.add_subsystem(f'cons_comp_{idx}', ConsComp(initial_DVs=initial_DVs, J=c.function, dJ=jac, con_name=con_name), promotes=['*']) - + lower_bounds = bounds[:, 0] upper_bounds = bounds[:, 1] - + # set up the optimization if 'SLSQP' in kwargs['opt_routine']: prob.driver = om.ScipyOptimizeDriver() @@ -200,30 +200,30 @@ def jac(x): folder_output = kwargs["options"]["folder"] prob.driver.opt_settings["Summary file"] = os.path.join(folder_output, "SNOPT_summary.out") prob.driver.opt_settings["Print file"] = os.path.join(folder_output, "SNOPT_print.out") - + if kwargs["options"]["verify_snopt"]: prob.driver.opt_settings["Verify level"] = 0 else: prob.driver.opt_settings["Verify level"] = -1 - + # prob.driver.opt_settings["Major iterations limit"] = 3 prob.model.add_design_var('DVs', lower=lower_bounds, upper=upper_bounds) prob.model.add_objective('obj', ref=kwargs["options"]["obj_ref"], ref0=kwargs["options"]["obj_ref0"]) - + for idx, constraint_type in enumerate(constraint_types): con_name = f'con_{idx}' if constraint_type == "eq": prob.model.add_constraint(con_name, equals=0.) else: # Inequality means it's positive from scipy and dolfin - prob.model.add_constraint(con_name, lower=0.) + prob.model.add_constraint(con_name, lower=0.) prob.setup() - + prob.set_val('DVs', initial_DVs) - + # Optional debugging step to check the total derivatives if kwargs["options"]["check_totals"]: prob.run_model() @@ -232,15 +232,15 @@ def jac(x): else: # Run the optimization prob.run_driver() - + # Return the optimal design variables return(prob['DVs']) class Optimizer(object): """ A GenericProblem contains on the basic functions required by all problem objects. - - Args: + + Args: dom (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. """ def __init__(self, solver): @@ -271,15 +271,15 @@ def __init__(self, solver): else: self.layout_bounds = np.array([[0,0],[0,0]]) self.layout_bounds = self.layout_bounds*self.xscale - + self.iteration = 0 self.fprint("Setting Up Optimizer",special="header") - + self.fprint("Controls: {0}".format(self.control_types)) self.CreateControls() - - + + self.fprint("Define Bounds") self.CreateBounds() @@ -307,7 +307,7 @@ def SetupConstraints(self): min_dist_dict = self.constraint_types.pop("min_dist") if "layout" in self.control_types: x_inds = self.indexes[0] - y_inds = self.indexes[1] + y_inds = self.indexes[1] RD = np.max(self.farm.get_rotor_diameters()) constraints.append(MinimumDistanceConstraint(x_inds, y_inds, min_dist_dict["target"]*RD, min_dist_dict["scale"])) @@ -340,7 +340,7 @@ def DebugOutput(self,iteration=0, m=[]): if hasattr(self,"gradients"): for i, d in enumerate(self.gradients): self.tag_output("grad_"+self.names[i],float(d)) - + ### TODO: Output taylor convergence data if hasattr(self,"conv_rate"): pass @@ -362,7 +362,7 @@ def RecomputeReducedFunctional(self): self.Jcurrent = self.J def ReducedFunctionalCallback(self, j, m): - self.Jcurrent = j + self.Jcurrent = j def CreateControls(self): @@ -370,7 +370,7 @@ def CreateControls(self): self.controls = [] self.control_pointers = [] self.names = [] - self.indexes = [[],[],[],[],[],[],[],[]] + self.indexes = [[],[],[],[],[],[],[],[],[],[],[],[]] self.init_vals = [] j = 0 @@ -408,6 +408,36 @@ def CreateControls(self): self.control_pointers.append(("axial",i)) self.init_vals.append(Constant(float(self.farm.turbines[i].maxial))) + if "power_curve" in self.control_types: + for i in self.solver.opt_turb_id: + self.indexes[8].append(j) + j+=1 + self.names.append("Prated_"+repr(i)) + self.controls.append(Control(self.farm.turbines[i].mPrated)) + self.control_pointers.append(("Prated",i)) + self.init_vals.append(Constant(float(self.farm.turbines[i].mPrated))) + + self.indexes[9].append(j) + j+=1 + self.names.append("Trated_"+repr(i)) + self.controls.append(Control(self.farm.turbines[i].mTrated)) + self.control_pointers.append(("Trated",i)) + self.init_vals.append(Constant(float(self.farm.turbines[i].mTrated))) + + self.indexes[10].append(j) + j+=1 + self.names.append("CPprime0_"+repr(i)) + self.controls.append(Control(self.farm.turbines[i].mCPprime0)) + self.control_pointers.append(("CPprime0",i)) + self.init_vals.append(Constant(float(self.farm.turbines[i].mCPprime0))) + + self.indexes[11].append(j) + j+=1 + self.names.append("CTprime0_"+repr(i)) + self.controls.append(Control(self.farm.turbines[i].mCTprime0)) + self.control_pointers.append(("CTprime0",i)) + self.init_vals.append(Constant(float(self.farm.turbines[i].mCTprime0))) + if "lift" in self.control_types: for i in self.solver.opt_turb_id: for k in range(self.farm.turbines[0].num_actuator_nodes): @@ -569,7 +599,7 @@ def Gradient(self): capture_memory = False if capture_memory: mem_out, der = memory_usage(self.Jhat.derivative,max_usage=True,retval=True,max_iterations=1) - else: + else: mem_out = 2*mem0 der = self.Jhat.derivative() @@ -673,7 +703,7 @@ def SaveControls(self,m): else: m_old.append(float(getattr(self.farm.turbines[index],control_name))) # print(m_new) - # print(type(m_new)) + # print(type(m_new)) for i in range(len(m_new)): if "yaw" in self.names[i]: @@ -719,7 +749,7 @@ def SaveFunctions(self): self.params.Save(u,"velocity",subfolder="OptSeries/",val=self.iteration,file=self.velocity_file) self.params.Save(p,"pressure",subfolder="OptSeries/",val=self.iteration,file=self.pressure_file) # self.params.Save(tf,"tf",subfolder="OptSeries/",val=self.iteration,file=self.tf_file) - + @no_annotations def OptPrintFunction(self,m,test=None): if test is not None: @@ -736,8 +766,8 @@ def OptPrintFunction(self,m,test=None): self.problem.farm.plot_farm(filename="wind_farm_step_"+repr(self.iteration),objective_value=self.Jcurrent) # if "chord" in self.control_types: - # c_lower = np.array(self.bounds[0])[self.indexes[6]] - # c_upper = np.array(self.bounds[1])[self.indexes[6]] + # c_lower = np.array(self.bounds[0])[self.indexes[6]] + # c_upper = np.array(self.bounds[1])[self.indexes[6]] # self.problem.farm.plot_chord(filename="chord_step_"+repr(self.iteration),objective_value=self.Jcurrent,bounds=[c_lower,c_upper]) self.DebugOutput(self.iteration, m) @@ -774,7 +804,7 @@ def Optimize(self): options["verify_snopt"] = self.verify_snopt if hasattr(self, 'check_totals'): options["check_totals"] = self.check_totals - + if self.opt_type == "minimize": opt_function = minimize elif self.opt_type == "maximize": @@ -782,14 +812,14 @@ def Optimize(self): else: raise ValueError(f"Unknown optimization type: {self.opt_type}") - ### optimize + ### optimize if "SNOPT" in self.opt_routine or "OM_SLSQP" in self.opt_routine: m_opt=opt_function(self.Jhat, method="Custom", tol=1.0e-8, options = options, constraints = self.merged_constraint, bounds = self.bounds, callback = self.OptPrintFunction, algorithm=om_wrapper, opt_routine=self.opt_routine) else: m_opt=opt_function(self.Jhat, method=self.opt_routine, tol=1.0e-8, options = options, constraints = self.merged_constraint, bounds = self.bounds, callback = self.OptPrintFunction) self.m_opt = m_opt - + if self.num_controls == 1: self.m_opt = [self.m_opt] self.OptPrintFunction(self.m_opt) @@ -804,7 +834,7 @@ def Optimize(self): # if "axial" in self.control_types: # new_values["a"] = m_f[self.indexes[3]] # self.problem.farm.UpdateControls(**new_values) - + # self.fprint("Solving With New Values") # self.solver.Solve() @@ -813,7 +843,7 @@ def Optimize(self): return self.m_opt def TaylorTest(self): - + self.fprint("Beginning Taylor Test",special="header") h = [] @@ -975,5 +1005,5 @@ def jacobian(self, m): out = np.empty((0,len(m))) for con in self.constraint_list: out = np.vstack((out, con.jacobian(m))) - + return out diff --git a/windse/ProblemManager.py b/windse/ProblemManager.py index 4fa6f398..238a0c94 100644 --- a/windse/ProblemManager.py +++ b/windse/ProblemManager.py @@ -1,5 +1,5 @@ """ -The ProblemManager contains all of the +The ProblemManager contains all of the different classes of problems that windse can solve """ @@ -12,7 +12,7 @@ main_file = os.path.basename(__main__.__file__) else: main_file = "ipython" - + ### This checks if we are just doing documentation ### if not main_file in ["sphinx-build", "__main__.py"]: from dolfin import * @@ -32,8 +32,8 @@ class GenericProblem(object): """ A GenericProblem contains on the basic functions required by all problem objects. - - Args: + + Args: domain (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. windfarm (:meth:`windse.WindFarmManager.GenericWindFarmm`): a windse windfarm object. function_space (:meth:`windse.FunctionSpaceManager.GenericFunctionSpace`): a windse function space object. @@ -44,7 +44,7 @@ def __init__(self,domain,windfarm,function_space,boundary_data): self.params = windse_parameters self.dom = domain self.farm = windfarm - self.fs = function_space + self.fs = function_space self.bd = boundary_data self.tf_first_save = True self.fprint = self.params.fprint @@ -59,7 +59,7 @@ def __init__(self,domain,windfarm,function_space,boundary_data): if isinstance(self.record_time,str): self.record_time = 0.0 - self.extra_kwarg = {} + self.extra_kwarg = {} if self.params.dolfin_adjoint: self.extra_kwarg["annotate"] = False @@ -128,7 +128,7 @@ def DebugOutput(self): self.tag_output("max_cd", np.max(cd)) self.tag_output("avg_cd", np.mean(cd)) self.tag_output("num_actuator_nodes", np.mean(num_actuator_nodes)) - + def ComputeTurbineForceTerm(self,u,v,inflow_angle,**kwargs): tf_start = time.time() @@ -237,7 +237,7 @@ def ChangeWindAngle(self,inflow_angle): """ This function recomputes all necessary components for a new wind direction - Args: + Args: inflow_angle (float): The new wind angle in radians """ adj_start = time.time() @@ -281,17 +281,17 @@ def UpdateActuatorLineControls(self, c_lift = None, c_drag = None, chord = None, yaw = float(yaw) self.farm.yaw[turb_i] = yaw self.farm.myaw[turb_i] = Constant(yaw) - + self.CopyALMtoWindFarm() class StabilizedProblem(GenericProblem): """ - The StabilizedProblem setup everything required for solving Navier-Stokes with + The StabilizedProblem setup everything required for solving Navier-Stokes with a stabilization term - Args: + Args: domain (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. windfarm (:meth:`windse.WindFarmManager.GenericWindFarmm`): a windse windfarm object. function_space (:meth:`windse.FunctionSpaceManager.GenericFunctionSpace`): a windse function space object. @@ -299,7 +299,7 @@ class StabilizedProblem(GenericProblem): """ def __init__(self,domain,windfarm,function_space,boundary_conditions): super(StabilizedProblem, self).__init__(domain,windfarm,function_space,boundary_conditions) - + ### Create Functional ### self.ComputeFunctional(self.dom.inflow_angle) self.DebugOutput() @@ -336,7 +336,7 @@ def ComputeFunctional(self,inflow_angle): ### These constants will be moved into the params file ### f = Constant((0.0,)*self.dom.dim) f.rename("f","f") - + nu = self.viscosity vonKarman=0.41 eps=Constant(self.stability_eps) @@ -358,7 +358,7 @@ def ComputeFunctional(self,inflow_angle): self.F += - inner(div(v),self.p_k)*dx self.F += - inner(div(self.u_k),q)*dx self.F += - inner(f,v)*dx - self.F += - tf_term + self.F += - tf_term # Add body force to functional @@ -379,14 +379,14 @@ def ComputeFunctional(self,inflow_angle): ######################################################## # self.F_sans_tf = (1.0)*inner(grad(self.u_k), grad(v))*dx - inner(div(v),self.p_k)*dx - inner(div(self.u_k),q)*dx - inner(f,v)*dx - # self.F = inner(grad(self.u_k)*self.u_k, v)*dx + (nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx - inner(div(v),self.p_k)*dx - inner(div(self.u_k),q)*dx - inner(f,v)*dx - inner(self.tf*(self.u_k[0]**2+self.u_k[1]**2),v)*dx - # stab_sans_tf = - eps*inner(grad(q), grad(self.p_k))*dx + # self.F = inner(grad(self.u_k)*self.u_k, v)*dx + (nu+self.nu_T)*inner(grad(self.u_k), grad(v))*dx - inner(div(v),self.p_k)*dx - inner(div(self.u_k),q)*dx - inner(f,v)*dx - inner(self.tf*(self.u_k[0]**2+self.u_k[1]**2),v)*dx + # stab_sans_tf = - eps*inner(grad(q), grad(self.p_k))*dx # self.F_sans_tf += stab ### Add in the Stabilizing term ### if abs(float(eps)) >= 1e-14: self.fprint("Using Stabilization Term") - stab = - eps*inner(grad(q), grad(self.p_k))*dx - eps*inner(grad(q), dot(grad(self.u_k), self.u_k))*dx + stab = - eps*inner(grad(q), grad(self.p_k))*dx - eps*inner(grad(q), dot(grad(self.u_k), self.u_k))*dx self.F += stab @@ -410,9 +410,9 @@ def ComputeFunctional(self,inflow_angle): class TaylorHoodProblem(GenericProblem): """ - The TaylorHoodProblem sets up everything required for solving Navier-Stokes + The TaylorHoodProblem sets up everything required for solving Navier-Stokes - Args: + Args: domain (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. windfarm (:meth:`windse.WindFarmManager.GenericWindFarmm`): a windse windfarm object. function_space (:meth:`windse.FunctionSpaceManager.GenericFunctionSpace`): a windse function space object. @@ -469,12 +469,12 @@ def ComputeFunctional(self,inflow_angle): self.F += - inner(div(self.v),self.p_k)*dx self.F += - inner(div(self.u_k),self.q)*dx # self.F += - inner(f,v)*dx - self.F += - tf_term + self.F += - tf_term # ### Add in the Stabilizing term ### # if abs(float(eps)) >= 1e-14: # self.fprint("Using Stabilization Term") - # stab = - eps*inner(grad(q), grad(self.p_k))*dx - eps*inner(grad(q), dot(grad(self.u_k), self.u_k))*dx + # stab = - eps*inner(grad(q), grad(self.p_k))*dx - eps*inner(grad(q), dot(grad(self.u_k), self.u_k))*dx # self.F += stab # Add body force to functional @@ -506,7 +506,7 @@ class IterativeSteady(GenericProblem): The IterativeSteady sets up everything required for solving Navier-Stokes using the SIMPLE algorithm - Args: + Args: domain (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. windfarm (:meth:`windse.WindFarmManager.GenericWindFarmm`): a windse windfarm object. function_space (:meth:`windse.FunctionSpaceManager.GenericFunctionSpace`): a windse function space object. @@ -548,7 +548,7 @@ def ComputeFunctional(self,inflow_angle): self.dp = Function(self.fs.Q) self.p_k_old = Function(self.fs.Q) - U_CN = 0.5*(u + self.u_k) + U_CN = 0.5*(u + self.u_k) # Adams-Bashforth velocity # U_AB = 1.5*u_k - 0.5*u_k_old # Time level k+1/2 @@ -559,7 +559,7 @@ def ComputeFunctional(self,inflow_angle): # ================================================================ - # FIXME: This up_k function is only present to avoid errors + # FIXME: This up_k function is only present to avoid errors # during assignments in GenericSolver.__init__ # Create the combined function space @@ -581,7 +581,7 @@ def ComputeFunctional(self,inflow_angle): # Solve for u_hat, a velocity estimate which doesn't include pressure gradient effects F1 = inner(dot(self.u_k, nabla_grad(u)), v)*dx \ + (nu+self.nu_T)*inner(grad(u), grad(v))*dx \ - - tf_term + - tf_term # Add body force to functional F1 += inner(-self.mbody_force*self.bd.inflow_unit_vector,v)*dx @@ -632,7 +632,7 @@ class UnsteadyProblem(GenericProblem): The UnsteadyProblem sets up everything required for solving Navier-Stokes using a fractional-step method with an adaptive timestep size - Args: + Args: domain (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. windfarm (:meth:`windse.WindFarmManager.GenericWindFarmm`): a windse windfarm object. function_space (:meth:`windse.FunctionSpaceManager.GenericFunctionSpace`): a windse function space object. @@ -686,7 +686,7 @@ def ComputeFunctional(self,inflow_angle): self.u_k1.assign(self.bd.bc_velocity) self.u_k2.assign(self.bd.bc_velocity) - # Calculate Reynolds stress + # Calculate Reynolds stress self.uk_sum = Function(self.fs.V, name="uk_sum") self.uk_sum.assign(self.dt_c*self.u_k) self.vertKE = Function(self.fs.Q, name="vertKE") @@ -713,7 +713,7 @@ def ComputeFunctional(self,inflow_angle): # ================================================================ - # FIXME: This up_k function is only present to avoid errors + # FIXME: This up_k function is only present to avoid errors # during assignments in GenericSolver.__init__ # Create the combined function space diff --git a/windse/SolverManager.py b/windse/SolverManager.py index c76dff93..dfb0b1d2 100644 --- a/windse/SolverManager.py +++ b/windse/SolverManager.py @@ -12,7 +12,7 @@ main_file = os.path.basename(__main__.__file__) else: main_file = "ipython" - + ### This checks if we are just doing documentation ### if not main_file in ["sphinx-build", "__main__.py"]: from dolfin import * @@ -56,12 +56,12 @@ ### Improved dolfin parameters ### parameters["std_out_all_processes"] = False; - parameters['form_compiler']['cpp_optimize_flags'] = '-O3 -fno-math-errno -march=native' + parameters['form_compiler']['cpp_optimize_flags'] = '-O3 -fno-math-errno -march=native' parameters["form_compiler"]["optimize"] = True parameters["form_compiler"]["cpp_optimize"] = True parameters['form_compiler']['representation'] = default_representation parameters['form_compiler']['quadrature_degree'] = windse_parameters["function_space"]["quadrature_degree"] - + class GenericSolver(object): """ A GenericSolver contains on the basic functions required by all solver objects. @@ -87,7 +87,7 @@ def __init__(self,problem): for key, value in self.params["optimization"].items(): setattr(self,key,value) - self.extra_kwarg = {} + self.extra_kwarg = {} if self.params.dolfin_adjoint: self.extra_kwarg["annotate"] = False @@ -107,7 +107,7 @@ def __init__(self,problem): if isinstance(self.opt_turb_id,int): self.opt_turb_id = [self.opt_turb_id] elif self.opt_turb_id == "all": - self.opt_turb_id = range(self.problem.farm.numturbs) + self.opt_turb_id = range(self.problem.farm.numturbs) elif isinstance(self.opt_turb_id, str): self.opt_turb_id = [int(self.opt_turb_id)] @@ -181,7 +181,7 @@ def ChangeWindSpeed(self,inflow_speed): """ This function recomputes all necessary components for a new wind direction - Args: + Args: theta (float): The new wind angle in radians """ self.problem.ChangeWindSpeed(inflow_speed) @@ -190,7 +190,7 @@ def ChangeWindAngle(self,inflow_angle): """ This function recomputes all necessary components for a new wind direction - Args: + Args: theta (float): The new wind angle in radians """ self.problem.ChangeWindAngle(inflow_angle) @@ -224,13 +224,23 @@ def ChangeWindAngle(self,inflow_angle): @no_annotations def EvaulatePowerFunctional(self): kwargs = { - "iter_val": self.iter_val, + "iter_val": self.iter_val, "simTime": self.simTime } out = self.problem.farm.save_power(self.problem.u_k,self.problem.dom.inflow_angle, **kwargs) return out + @no_annotations + def EvaulateThrustFunctional(self): + kwargs = { + "iter_val": self.iter_val, + "simTime": self.simTime + } + out = self.problem.farm.save_thrust(self.problem.u_k,self.problem.dom.inflow_angle, **kwargs) + return out + + def EvaluateObjective(self,output_name="objective_data",opt_iter=-1): self.fprint("Evaluating Objective Data",special="header") start = time.time() @@ -239,7 +249,7 @@ def EvaluateObjective(self,output_name="objective_data",opt_iter=-1): if self.J_saved: first_call = False - annotate = self.params.dolfin_adjoint + annotate = self.params.dolfin_adjoint ### Iterate over objectives ### obj_list = [opt_iter, self.iter_val, self.simTime] @@ -254,7 +264,7 @@ def EvaluateObjective(self,output_name="objective_data",opt_iter=-1): kwargs.update(obj_kwargs) out = obj_funcs._annotated_objective(objective_func, *args, **kwargs) obj_list.append(out) - J = obj_list[3] #grab first objective + J = obj_list[3] #grab first objective # ### Flip the sign because the objective is minimized but these values are maximized # for i in range(1,len(obj_list)): @@ -295,7 +305,7 @@ class SteadySolver(GenericSolver): """ This solver is for solving the steady state problem - Args: + Args: problem (:meth:`windse.ProblemManager.GenericProblem`): a windse problem object. """ def __init__(self,problem): @@ -322,7 +332,7 @@ def Solve(self): # exit() #################################################################### - ### This is the better way to define a nonlinear problem but it + ### This is the better way to define a nonlinear problem but it ### doesn't play nice with dolfin_adjoint ### Define Jacobian ### # dU = TrialFunction(self.problem.fs.W) @@ -391,7 +401,7 @@ def Solve(self): "relative_tolerance": 1e-5, "error_on_nonconvergence": False, "krylov_solver": krylov_options} - + solver_parameters = {"nonlinear_solver": "newton", "newton_solver": newton_options} @@ -413,7 +423,7 @@ def Solve(self): "error_on_nonconvergence": False, "line_search": "bt", #[basic,bt,cp,l2,nleqerr] "krylov_solver": krylov_options - }} + }} else: raise ValueError("Unknown nonlinear solver type: {0}".format(self.nonlinear_solver)) @@ -484,12 +494,15 @@ def Solve(self): ### Evaluate the objectives ### if self.optimizing or self.save_objective: self.J += self.EvaluateObjective() - # self.J += self.objective_func(self,(self.iter_theta-self.problem.dom.inflow_angle)) + # self.J += self.objective_func(self,(self.iter_theta-self.problem.dom.inflow_angle)) self.J = self.control_updater(self.J, self.problem) if self.save_power: self.EvaulatePowerFunctional() + if self.save_thrust: + self.EvaulateThrustFunctional() + # print(self.outflow_markers) # self.J += -dot(self.problem.farm.rotor_disks,self.u_k)*dx @@ -509,14 +522,14 @@ def Solve(self): self.DebugOutput() -# +# # ================================================================ class IterativeSteadySolver(GenericSolver): """ This solver is for solving the iterative steady state problem - Args: + Args: problem (:meth:`windse.ProblemManager.GenericProblem`): a windse problem object. """ def __init__(self,problem): @@ -553,7 +566,7 @@ def set_solver_mode(solver_type): return sor_vel, sor_pr def get_relaxation_param(a_min, a_max, a_c, a_cw, x): - + # a_min = minimum sor value # a_max = maximum sor value # a_c = crossover location @@ -563,7 +576,7 @@ def get_relaxation_param(a_min, a_max, a_c, a_cw, x): alpha = 1.0/(1.0+np.exp(alpha_exp)) alpha = alpha*(a_max-a_min) + a_min alpha = float(alpha) - + return alpha ### Save Files before solve ### @@ -760,7 +773,7 @@ def get_relaxation_param(a_min, a_max, a_c, a_cw, x): ### Evaluate the objectives ### if self.optimizing or self.save_objective: self.J += self.EvaluateObjective() - # self.J += self.objective_func(self,(self.iter_theta-self.problem.dom.inflow_angle)) + # self.J += self.objective_func(self,(self.iter_theta-self.problem.dom.inflow_angle)) self.J = self.control_updater(self.J, self.problem) if self.save_power and "power": @@ -800,7 +813,7 @@ class UnsteadySolver(GenericSolver): This solver can only be used if an unsteady problem has been specified in the input file. - Args: + Args: problem (:meth:`windse.ProblemManager.GenericProblem`): a windse problem object. """ def __init__(self,problem): @@ -824,7 +837,7 @@ def Solve(self): #self.final_time = 10.0 # ================================================================ - + # Start a counter for the total simulation time self.fprint("dt: %.4f" % (self.problem.dt)) self.fprint("Final Time: %.1f" % (self.final_time)) @@ -1038,7 +1051,7 @@ def save_adj(adj): tic = time.time() # solve(self.problem.a2==self.problem.L2, self.problem.p_k, bcs=self.problem.bd.bcp) # solve(self.problem.a2==self.problem.L2, self.problem.p_k, bcs=self.problem.bd.bcp, solver_parameters={"linear_solver": "gmres","preconditioner": "petsc_amg"}) - + b2 = assemble(self.problem.L2, tensor=b2) [bc.apply(b2) for bc in self.problem.bd.bcp] if self.optimizing: @@ -1055,7 +1068,7 @@ def save_adj(adj): tic = time.time() # solve(self.problem.a3==self.problem.L3, self.problem.u_k) # solve(self.problem.a3==self.problem.L3, self.problem.u_k, solver_parameters={"linear_solver": "cg","preconditioner": "jacobi"}) - + b3 = assemble(self.problem.L3, tensor=b3) if self.optimizing: # solve(A3, self.problem.u_k.vector(), b3, 'gmres', 'default') @@ -1085,13 +1098,13 @@ def save_adj(adj): self.problem.uk_sum.assign(self.problem.uk_sum+self.problem.dt_c*self.problem.u_k) elif self.simTime >= self.record_time: self.problem.vertKE = (self.problem.u_k[0]-self.problem.uk_sum[0]/(self.record_time-self.u_avg_time))*(self.problem.u_k[2]-self.problem.uk_sum[2]/(self.record_time-self.u_avg_time))*(self.problem.uk_sum[0]/(self.record_time-self.u_avg_time)) - + if self.save_all_timesteps: self.SaveTimeSeries(self.simTime,simIter) elif save_next_timestep: # Read in new inlet values # bcu = self.updateInletVelocityFromFile(saveCount, bcu) - + # Clean up self.simTime to avoid accumulating round-off error saveCount += 1 self.simTime = self.save_interval*saveCount @@ -1124,10 +1137,10 @@ def save_adj(adj): # # t2 = time.time() # # print(t2-t1) - # Power [=] N*m*rads/s + # Power [=] N*m*rads/s # self.problem.alm_power = self.problem.rotor_torque*(2.0*np.pi*self.problem.rpm/60.0) # self.problem.alm_power_dolfin = self.problem.rotor_torque_dolfin*(2.0*np.pi*self.problem.rpm/60.0) - + # # self.problem.alm_power_sum += self.problem.alm_power*self.problem.dt # # self.problem.alm_power_average = self.problem.alm_power_sum/self.simTime @@ -1197,7 +1210,7 @@ def save_adj(adj): self.adj_time_list.append(self.simTime) self.J += float(self.problem.dt)*J_next - self.problem.dt_sum += self.problem.dt + self.problem.dt_sum += self.problem.dt J_new = float(self.J/self.problem.dt_sum) ### TODO, replace this with an actual stabilization criteria such as relative difference tolerance @@ -1441,14 +1454,14 @@ def rot_y(theta): Ry = np.array([[np.cos(theta), 0, np.sin(theta)], [0, 1, 0], [-np.sin(theta), 0, np.cos(theta)]]) - + return Ry def rot_z(theta): Rz = np.array([[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]]) - + return Rz #================================================================ @@ -1585,68 +1598,68 @@ def rot_z(theta): # Set if using local velocity around inidividual nodes using_local_velocity = False - + if using_local_velocity: # Generate the fluid velocity from the actual node locations in the flow u_fluid = np.zeros((3, num_actuator_nodes)) - + for k in range(num_actuator_nodes): u_fluid[:, k] = self.problem.u_k1(xblade_rotated[0, k], xblade_rotated[1, k], xblade_rotated[2, k]) - + else: # Generate the fluid velocity analytically using the hub height velocity # u_inf_vec = u_inf*np.ones(num_actuator_nodes) - + # u_fluid = np.vstack((u_inf_vec, # np.zeros(num_actuator_nodes), # np.zeros(num_actuator_nodes))) u_fluid = np.zeros((3, num_actuator_nodes)) - + for k in range(num_actuator_nodes): u_fluid[0, k] = 8.0*(xblade_rotated[2, k]/hub_height)**0.18 - + # Rotate the blade velocity in the global x, y, z, coordinate system blade_vel_rotated = np.dot(Rz, np.dot(Rx, -blade_vel)) - + # Form the total relative velocity vector (including velocity from rotating blade) u_rel = u_fluid + blade_vel_rotated - + # Create unit vectors in the direction of u_rel u_rel_mag = np.linalg.norm(u_rel, axis=0) u_rel_mag[u_rel_mag < 1e-6] = 1e-6 u_unit = u_rel/u_rel_mag - + # Calculate the lift and drag forces using the relative velocity magnitude lift = (0.5*cl*rho*c*w*u_rel_mag**2)/(eps**3 * np.pi**1.5) drag = (0.5*cd*rho*c*w*u_rel_mag**2)/(eps**3 * np.pi**1.5) - + # Calculate the force at every mesh point due to every node [numGridPts x NumActuators] nodal_lift = lift*np.exp(-dist2/eps**2) nodal_drag = drag*np.exp(-dist2/eps**2) - + # Calculate a vector in the direction of the blade - blade_unit = xblade_rotated[:, -1] - np.array([0.0, 0.0, hub_height]) - + blade_unit = xblade_rotated[:, -1] - np.array([0.0, 0.0, hub_height]) + for k in range(num_actuator_nodes): # The drag unit simply points opposite the relative velocity unit vector drag_unit = -u_unit[:, k] - + # The lift is normal to the plane generated by the blade and relative velocity lift_unit = np.cross(drag_unit, blade_unit) lift_unit_mag = np.linalg.norm(lift_unit) if lift_unit_mag < 1e-6: lift_unit_mag = 1e-6 lift_unit = lift_unit/lift_unit_mag - + vector_nodal_drag = np.outer(nodal_drag[:, k], drag_unit) vector_nodal_lift = np.outer(nodal_lift[:, k], lift_unit) drag_force += vector_nodal_drag lift_force += vector_nodal_lift - + # The total turbine force is the sum of lift and drag effects turbine_force = drag_force + lift_force @@ -1739,7 +1752,7 @@ def UpdateActuatorLineForceOld(self, simTime): for theta_0 in theta_vec: theta = theta_0 + theta_offset - + # Create rotation matrix for this turbine blade rotA = np.array([[1, 0, 0], [0, np.cos(theta), -np.sin(theta)], @@ -1786,7 +1799,7 @@ def UpdateActuatorLineForceOld(self, simTime): tf_vec[3*k+0] += -drag_sum tf_vec[3*k+1] += lift_sum*np.sin(theta) tf_vec[3*k+2] += -lift_sum*np.cos(theta) - + tf_vec[np.abs(tf_vec) < 1e-12] = 0.0 self.problem.tf.vector()[:] = tf_vec @@ -1821,7 +1834,7 @@ def UpdateTurbineForce(self, simTime, turbsPerPlatform): # Position of the kth turbune xpos = float(self.problem.farm.mx[k]) ypos = float(self.problem.farm.my[k]) - + if self.problem.dom.dim == 2: x0 = np.array([xpos, ypos]) else: @@ -2034,10 +2047,10 @@ class MultiAngleSolver(SteadySolver): This solver will solve the problem using the steady state solver for every angle in angles. - Args: + Args: problem (:meth:`windse.ProblemManager.GenericProblem`): a windse problem object. angles (list): A list of wind inflow directions. - """ + """ def __init__(self,problem): super(MultiAngleSolver, self).__init__(problem) @@ -2080,10 +2093,10 @@ class TimeSeriesSolver(SteadySolver): This solver will solve the problem using the steady state solver for every angle in angles. - Args: + Args: problem (:meth:`windse.ProblemManager.GenericProblem`): a windse problem object. angles (list): A list of wind inflow directions. - """ + """ def __init__(self,problem): super(TimeSeriesSolver, self).__init__(problem) @@ -2153,7 +2166,7 @@ def Solve(self): # W = self.problem.farm.thickness[0]*1.0 - # R = self.problem.farm.RD[0]/2.0 + # R = self.problem.farm.RD[0]/2.0 # T_norm = 2.0*gamma(7.0/6.0) # D_norm = pi*gamma(4.0/3.0) @@ -2173,7 +2186,7 @@ def Solve(self): # print(dim_scale) # dim_scale = 2*R*R*volNormalization_2D/volNormalization_3D # print(dim_scale) - # else: + # else: # dim_scale = 1.0 # ### Calculate Power ### @@ -2204,7 +2217,7 @@ def Solve(self): # # ### Set some Values ### # # yaw = self.problem.farm.myaw[i]+inflow_angle # # W = self.problem.farm.W[i] - # # R = self.problem.farm.RD[i]/2.0 + # # R = self.problem.farm.RD[i]/2.0 # # A = pi*R**2.0 # # a = self.problem.farm.ma[i] # # C_tprime = 4*a/(1-a) @@ -2268,7 +2281,7 @@ def Solve(self): # f.close() # return J - + # # def CalculatePowerFunctional_o(self,delta_yaw = 0.0): # # # def CalculatePowerFunctional(self,delta_yaw = 0.0): # # self.fprint("Computing Power Functional") @@ -2283,14 +2296,14 @@ def Solve(self): # # mz = self.problem.farm.mz[i] # # x0 = [mx,my,mz] # # W = self.problem.farm.thickness[i]*1.0 - # # R = self.problem.farm.RD[i]/2.0 + # # R = self.problem.farm.RD[i]/2.0 # # ma = self.problem.farm.ma[i] # # yaw = self.problem.farm.myaw[i]+delta_yaw # # u = self.u_next # # A = pi*R**2.0 # # C_tprime = 4*ma/(1-ma) # # C_pprime = 0.45/(1-ma)**3 - + # # ### Rotate and Shift the Turbine ### # # xs = self.problem.farm.YawTurbine(x,x0,yaw) # # u_d = u[0]*cos(yaw) + u[1]*sin(yaw) @@ -2360,7 +2373,7 @@ def Solve(self): # # ### Create the function that represents the force ### # # if self.problem.farm.force == "constant": - # # F = 0.5*self.problem.farm.RD[i]*C_tprime + # # F = 0.5*self.problem.farm.RD[i]*C_tprime # # elif self.problem.farm.force == "sine": # # F = 0.5*self.problem.farm.RD[i]*C_tprime*(r*sin(pi*r)+0.5)/(.81831) @@ -2386,12 +2399,12 @@ def Solve(self): # # ### LIST O PROBLEMS ### # # # 1. Unfortunately, this is not a form so cannot be used with dolfin_adjoint (this just returns a float) # # # 2. Where did WTGbase go? I know it need to be a scaler but does yaw angle not matter in 2D? - # # # 3. Should the C_pprime term be inside the numerator? + # # # 3. Should the C_pprime term be inside the numerator? # # # 4. Are you positive you want to divide by the total force (F*T*D) instead of just the space kernal (T*D) - + diff --git a/windse/default_parameters.yaml b/windse/default_parameters.yaml index 0b08e944..1d34785b 100644 --- a/windse/default_parameters.yaml +++ b/windse/default_parameters.yaml @@ -1,23 +1,23 @@ # These are all the available parameters. For more details visit https://windse.readthedocs.io/en/latest/params.html general: # parameters for ParameterManager - name: Null # The folder name for the current run + name: Null # The folder name for the current run preappend_datetime: False # if true, the time and date will be preappened to the output folder output: [solution] # These are the fields that can be saved during/after the simulation, choices: "mesh", "initial_guess", "height", "turbine_force", "solution", "debug" output_folder: output/ # This is the root folder for all the output files output_type: pvd # this is the filetype various fields are saved it, choices "pvd", "xdmf" (note "xdmf" is better for parallel, but "pvd" can be opened during simulation) - dolfin_adjoint: False # If true, this will import dolfin_adjoint which is required for calculating gradients. Must be true for optimizations + dolfin_adjoint: False # If true, this will import dolfin_adjoint which is required for calculating gradients. Must be true for optimizations debug_mode: False # If true, a file, tagged_output.yaml, in the root output folder. used for recursive tests domain: # parameters for DomainManager - type: Null # The type of domain + type: Null # The type of domain path: Null # The path to all domain files if using imported domain and standard naming mesh_path: Null # a specific path for the imported mesh file terrain_path: Null # a specific path for the terrain file if using complex/interpolated terrain bound_path: Null # a specific path to a MeshFunction file that store the boundary IDs filetype: xml.gz # The file type for imported domains. Choices: xml, xml.gz, h5 scaled: False # Attempt scale the domain to km instead of m. (Extremely experimental, do not use) - ground_reference: 0.0 # the z-coordinate of the ground. + ground_reference: 0.0 # the z-coordinate of the ground. streamwise_periodic: False # sets periodic boundary conditions along the x-direction (Not yet implemented) spanwise_periodic: False # sets periodic boundary conditions along the y-direction (Not yet implemented) x_range: Null # the extents, in meters, of the domain in the streamwise direction @@ -45,8 +45,8 @@ domain: # parameters for DomainManager my: Null # y slope wind_farm: # parameters for WindFarmManager - type: Null # type of wind farm (grid, random, imported, empty) - path: Null # location of imported wind farm + type: Null # type of wind farm (grid, random, imported, empty) + path: Null # location of imported wind farm display: False # If true, then matplotlib with show() the wind farm/chord profiles mid run ex_x: Null # extents of the farm in the x direction in meters ex_y: Null # extents of the farm in the y direction in meters @@ -71,10 +71,10 @@ turbines: # parameters for the turbine objects force: sine # distribution of force along the radial direction of an actuator disk. choices: constant, sine, ?chord? rpm: Null # rotations per minute for the alm method read_turb_data: Null # location of alm data - blade_segments: computed # number of nodes along the rotor radius + blade_segments: computed # number of nodes along the rotor radius use_local_velocity: True # use the velocity at the rotor to compute alm forces (otherwise use inflow) max_chord: 1000 # upper limit when optimizing chord - chord_factor: 1.0 # This multiplies all the chords by a constant factor, e.g., 2.0 makes a chord that's twice as thick everywhere + chord_factor: 1.0 # This multiplies all the chords by a constant factor, e.g., 2.0 makes a chord that's twice as thick everywhere gauss_factor: 2.0 # This is the factor that gets multiplied by the minimum mesh spacing to set the gaussian width, e.g., gaussian_width = 2.0*dx_min tip_loss: True # Determines whether or not a tip-loss model is used in the calculation of the ALM force (False means no tip loss is modeled) hub_rad: 0.0 # The radius of the hub. If non-zero, actuator nodes will still be placed in the range [0, rotor_radius], but the lift/drag properties in the range [0, hub_rad] will be modified to reflect the blade root @@ -84,7 +84,7 @@ turbines: # parameters for the turbine objects motion_file: Null # Location to the platform motion data motion_type: Null # Type of motion to apply can be single string or list of: 'surge', 'sway', 'heave', 'roll', 'pitch', and/or 'yaw' use_gauss_vel_probe: False # Prob velocity at ALM nodes using an gaussian sphere rather than eval() - use_ap_linear_interp: False # Uses linear interpolation when building the airfoil polars. + use_ap_linear_interp: False # Uses linear interpolation when building the airfoil polars. refine: # parameters for RefinementManager warp_type: Null # warping will shift the nodes along the z direction concentrating them near the ground. choices: "smooth", "split" @@ -93,17 +93,17 @@ refine: # parameters for RefinementManager warp_height: Null # for split warps, where the spit happens farm_num: 0 # number of farm level refinements farm_type: square # type of farm level refinements. choices: "full", "box", "cylinder", "stream" - farm_factor: 1.0 # scaling factor for the size of the refinement + farm_factor: 1.0 # scaling factor for the size of the refinement turbine_num: 0 # number of turbine level refinements turbine_type: simple # type of turbine refinement. choices: "simple", "tear", "wake" - turbine_factor: 1.0 # scaling factor for the size of the refinement + turbine_factor: 1.0 # scaling factor for the size of the refinement refine_custom: Null # allows for a list of custom refine commands refine_power_calc: False # bare minimum refinement around turbines to increase power calculation accuracy function_space: # parameters for the FunctionSpaceManager type: Null # type of function space quadrature_degree: 6 # used when calculating integrals - turbine_space: Quadrature # used with numpy turbine_method, sets the space the turbine are calculate on + turbine_space: Quadrature # used with numpy turbine_method, sets the space the turbine are calculate on turbine_degree: 6 # used with numpy turbine_method, sets degree boundary_conditions: # parameters for the BoundaryManager @@ -123,7 +123,7 @@ boundary_conditions: # parameters for the BoundaryManager top: Null # positive z inflow: Null # inflow, used in cylinder/circle outflow: Null # outflow, used in cylinder/circle - boundary_types: # used for changing the boundary types + boundary_types: # used for changing the boundary types inflow: Null # list of inflow bc names no_slip: Null # list of no_slip bc names free_slip: Null # list of free_slip bc names @@ -149,9 +149,10 @@ solver: # parameters for the SolverManager endpoint: False # include the final wind angle in the sweep velocity_path: Null # location of inflow velocities for multivelocity save_power: True # save the power data to the data folder + save_thrust: True # save the thrust data to the data folder nonlinear_solver: snes # type of nonlinear solver. choices: snes, newton newton_relaxation: 1.0 # relaxation parameter (0,1] for newton - cfl_target: 0.5 # target cfl number for unsteady solve + cfl_target: 0.5 # target cfl number for unsteady solve cl_iterator: 0 # debugging tool, do not use save_all_timesteps: False # Save fields at every time step @@ -166,13 +167,13 @@ optimization: # parameters for OptimizationManager scale: 1 # sets the scale for the constraint, use to match objective magnitude. use -1 to flip constraint to target-c(m)>=0 ==> c(m) <= target kwargs: Null # if constraint is based on an objective function with kwargs, set them here save_objective: False # save the objective do a file in data/ - opt_turb_id : all # which turbines to optimize, int or list or "all" - record_time: 0.0 # when to start recording for unsteady objectives float or "computed" + opt_turb_id : all # which turbines to optimize, int or list or "all" + record_time: 0.0 # when to start recording for unsteady objectives float or "computed" alm_DELs: DEL_start_time: 0 u_avg_time: 0 # when to start averaging velocity for use in objective functions opt_routine: SLSQP # optimization method, choices: SLSQP, L-BFGS-B, OM_SLSQP, SNOPT (requires custom install) - obj_ref: 1. # Sets the value of the objective function that will be treated as 1 by the SNOPT driver + obj_ref: 1. # Sets the value of the objective function that will be treated as 1 by the SNOPT driver obj_ref0: 0. # Sets the value of the objective function that will be treated as 0 by the SNOPT driver taylor_test: False # Run the taylor test optimize: False # optimize the problem diff --git a/windse/input_schema.yaml b/windse/input_schema.yaml index a0579247..49c17cd2 100644 --- a/windse/input_schema.yaml +++ b/windse/input_schema.yaml @@ -444,13 +444,13 @@ properties: units: "dimensionless" use_local_tf_dx: default: False - type: + type: - "boolean" description: "A flag to assemble the tf using a local dx measure for each turbine." units: "dimensionless" local_dx_scaling: default: 2.5 - type: + type: - "number" description: "Controls how big the local domain integration is in rotor diameters." units: "dimensionless" @@ -470,6 +470,7 @@ properties: - "2D_disk" - "numpy_disk" - "disk" + - "power_disk" - "hybrid_disk" - "line" - "expr_disk" @@ -970,6 +971,10 @@ properties: default: True type: "boolean" description: "Save the power data to the data folder." + save_thrust: + default: False + type: "boolean" + description: "Save the thrust data to the data folder." nonlinear_solver: type: "string" description: "Type of nonlinear solver, options are ``snes`` or ``newton``." diff --git a/windse/turbine_types/ActuatorDisk.py b/windse/turbine_types/ActuatorDisk.py index fbdf81f1..522a6053 100644 --- a/windse/turbine_types/ActuatorDisk.py +++ b/windse/turbine_types/ActuatorDisk.py @@ -12,12 +12,12 @@ class ActuatorDisk(GenericTurbine): def __init__(self, i,x,y,dom,imported_params=None): - # Define the acceptable column of an wind_farm.csv imput file + # Define the acceptable column of an wind_farm.csv imput file self.yaml_inputs = ["HH", "RD", "yaw", "thickness", "axial", "force"] # Init turbine super(ActuatorDisk, self).__init__(i,x,y,dom,imported_params) - + def load_parameters(self): self.HH = self.params["turbines"]["HH"] @@ -26,7 +26,7 @@ def load_parameters(self): self.thickness = self.params["turbines"]["thickness"] self.axial = self.params["turbines"]["axial"] self.force = self.params["turbines"]["force"] - + def create_controls(self): self.mx = Constant(self.x, name="x_{:d}".format(self.index)) self.my = Constant(self.y, name="y_{:d}".format(self.index)) @@ -34,7 +34,7 @@ def create_controls(self): self.maxial = Constant(self.axial, name="axial_{:d}".format(self.index)) # The control list is very important for extracting information from the optimzation step. - self.controls_list = ["x","y","yaw","axial"] + self.controls_list = ["x","y","yaw","axial"] def build_actuator_disk(self): @@ -52,11 +52,11 @@ def build_actuator_disk(self): T_norm = 2.0*gamma(7.0/6.0) if self.dom.dim == 3: WTGbase = as_vector((cos(yaw),sin(yaw),0.0)) - A = np.pi*R**2.0 + A = np.pi*R**2.0 D_norm = np.pi*gamma(4.0/3.0) else: WTGbase = as_vector((cos(yaw),sin(yaw))) - A = 2*R + A = 2*R D_norm = 2.0*gamma(7.0/6.0) ### Rotate and Shift the Turbine ### @@ -98,11 +98,11 @@ def build_actuator_disk(self): def turbine_force(self,u,inflow_angle,fs): """ - This function creates a turbine force by applying - a spacial kernel to each turbine. This kernel is + This function creates a turbine force by applying + a spacial kernel to each turbine. This kernel is created from the turbines location, yaw, thickness, diameter, and force density. Currently, force density is limit to a scaled - version of + version of .. math:: @@ -167,7 +167,7 @@ def prepare_saved_functions(self, func_list): - #### These two chord based functions are still very experimental + #### These two chord based functions are still very experimental def radial_chord_force(r,chord): cx = np.linspace(-0.1,1.25,len(chord)) @@ -182,7 +182,7 @@ def radial_chord_force(r,chord): force = 0 int_force = 0 for i in range(len(chord)): - Li = 1.0 + Li = 1.0 for j in range(len(chord)): if i!=j: Li *= (r - cx[j])/(cx[i]-cx[j]) diff --git a/windse/turbine_types/ActuatorDiskSimplePowerCurve.py b/windse/turbine_types/ActuatorDiskSimplePowerCurve.py new file mode 100644 index 00000000..e56e77e9 --- /dev/null +++ b/windse/turbine_types/ActuatorDiskSimplePowerCurve.py @@ -0,0 +1,219 @@ +# import windse function +from typing import Generic +from . import GenericTurbine + +# import dolfin functions +from . import Constant, SpatialCoordinate, as_vector, cos, sin, exp, ln, sqrt, dot, assemble, dx + +from . import ActuatorDisk + +import numpy as np +from scipy.special import gamma + +class ActuatorDiskSimplePowerCurve(GenericTurbine): + + def __init__(self, i, x, y, dom, imported_params=None): + # Define the acceptable column of an wind_farm.csv imput file + self.yaml_inputs = ["HH", "RD", "yaw", "thickness", "CTprime0", "CPprime0", "Prated", "force"] + + # Init turbine + super(ActuatorDiskSimplePowerCurve, self).__init__(i,x,y,dom,imported_params) + + def load_parameters(self): + self.HH = self.params["turbines"]["HH"] + self.RD = self.params["turbines"]["RD"] + self.yaw = self.params["turbines"]["yaw"] + self.thickness = self.params["turbines"]["thickness"] + self.axial = -1 # needed for generic turbine + self.CPprime0 = self.params["turbines"]["CPprime0"] + self.CTprime0 = self.params["turbines"]["CTprime0"] + self.Prated = self.params["turbines"]["Prated"] + self.force = self.params["turbines"]["force"] + + # calibration factors + self.calibration_factor_CTprime0 = 1.0 # 0.872882971129596 # based on single-turbine thrust-power-curve calibration + self.calibration_factor_CPprime0 = 1.0 # based on single-turbine thrust-power-curve calibration + self.calibration_factor_Prated = 1.0 # based on single-turbine thrust-power-curve calibration + self.calibration_factor_CPprime0 = 1.0 # 0.8777933671305554 # based on single-turbine thrust-power-curve calibration + + def create_controls(self): + self.mx = Constant(self.x, name="x_{:d}".format(self.index)) + self.my = Constant(self.y, name="y_{:d}".format(self.index)) + self.myaw = Constant(self.yaw, name="yaw_{:d}".format(self.index)) + self.mCPprime0 = Constant(self.CPprime0, name="CPprime0_{:d}".format(self.index)) + self.mCTprime0 = Constant(self.CTprime0, name="CTprime0_{:d}".format(self.index)) + self.mPrated = Constant(self.Prated, name="Prated_{:d}".format(self.index)) + + # The control list is very important for extracting information from the optimzation step. + self.controls_list = ["x","y","yaw","CPprime0","CTprime0","Prated"] + + def build_actuator_disk(self): + + ### Alias useful values ### + x0 = [self.mx,self.my,self.mz] + yaw = self.myaw+self.inflow_angle + W = self.thickness*1.0 + R = self.RD/2.0 + CTprime0 = self.mCTprime0 + # CPprime0 = self.mCPprime0 + x=SpatialCoordinate(self.dom.mesh) + + ### Set up some dim dependent values ### + S_norm = (2.0+np.pi)/(2.0*np.pi) + T_norm = 2.0*gamma(7.0/6.0) + if self.dom.dim == 3: + WTGbase = as_vector((cos(yaw),sin(yaw),0.0)) + A = np.pi*R**2.0 + D_norm = np.pi*gamma(4.0/3.0) + else: + WTGbase = as_vector((cos(yaw),sin(yaw))) + A = 2*R + D_norm = 2.0*gamma(7.0/6.0) + + ### Rotate and Shift the Turbine ### + xs = self.yaw_turbine(x,x0,yaw) + + ### Create the function that represents the Thickness of the turbine ### + T = exp(-pow((xs[0]/W),6.0)) + + ### Create the function that represents the Disk of the turbine + r = sqrt(xs[1]**2.0+xs[2]**2.0)/R + D = exp(-pow(r,6.0)) + + ### Create the function that represents the force ### + if self.force == "constant": + force = 1.0 + elif self.force == "sine": + force = (r*sin(np.pi*r)+0.5)/S_norm + elif self.force == "chord": + raise NotImplementedError("not implemented! -cfrontin") + # self._setup_chord_force() + # chord = self.mchord[i] + # force = self.radial_chord_force(r,chord) + F = -0.5*A*(self.calibration_factor_CTprime0*CTprime0)*force + + ### Calculate normalization constant ### + volNormalization = T_norm*D_norm*W*R**(self.dom.dim-1) + # volNormalization = assemble(T*D*dx) + + # self.fprint(f"Turbine {self.index}:") + # self.fprint(f"Assembled Volume: {assemble(T*D*dx)}",offset=1) + # self.fprint(f"Computed Volume: {volNormalization}",offset=1) + # self.fprint(f"Assembled Force: {assemble(F*T*D/volNormalization*dx)}",offset=1) + + # compute disk averaged velocity in yawed case and don't project + actuator_disks=F*T*D*WTGbase/volNormalization + # actuator_disks=project(actuator_disks,self.fs.V,solver_type='cg',preconditioner_type="hypre_amg") + + + return actuator_disks + + def turbine_force(self,u,inflow_angle,fs): + """ + This function creates a turbine force by applying + a spacial kernel to each turbine. This kernel is + created from the turbines location, yaw, thickness, diameter, + and force density. Currently, force density is limit to a scaled + version of + + .. math:: + + r\\sin(r), + + where :math:`r` is the distance from the center of the turbine. + + Args: + u: the velocity field + inflow_angle: direction the wind is blowing (0 is west to east and positive rotates clockwise) + fs: the function space manager object + + Returns: + tf (dolfin.Function): the turbine force. + + Todo: + * Setup a way to get the force density from file + """ + + ### compute the space kernal and radial force + self.fs = fs + self.inflow_angle.assign(inflow_angle) + self.actuator_disk = self.build_actuator_disk() + + ### Expand the dot product + yaw = self.myaw+self.inflow_angle + self.tf1 = self.actuator_disk * cos(yaw)**2 + self.tf2 = self.actuator_disk * sin(yaw)**2 + self.tf3 = self.actuator_disk * 2.0 * cos(yaw) * sin(yaw) + + # compute power curve adjustment as a function of velocity + CPprime0 = self.mCPprime0 + CTprime0 = self.mCTprime0 + # `-> created to match a single-turbine test case + Prated = self.mPrated + A = np.pi/4.0*self.RD**2.0 + Vrated3 = Prated/(0.5*CPprime0*A) # cubed rated velo + + # precompute key values + beta_smooth = 128.0 + vel_magnitude = sqrt(u[0]**2 + u[1]**2 + u[2]**2) + + f2 = Vrated3/(vel_magnitude**3) + + # smooth once: one in power space and one in coefficient space + + # # first attempt: boltzmann operator + # # first blend: smoothmin Region II.5 and Region III power + # blend1 = 1.0*( + # f2*exp(-beta_smooth*f2) + f3*exp(-beta_smooth*f3) + # )/( + # exp(-beta_smooth*f2) + exp(-beta_smooth*f3) + # ) + # + # # second blend: smoothmin Region II and first blend coefficient + # CTp_factor = ( + # f1/f0*exp(-beta_smooth*f1/f0) + blend1/f0*exp(-beta_smooth*blend1/f0) + # )/( + # exp(-beta_smooth*f1/f0) + exp(-beta_smooth*blend1/f0) + # ) + + # # second attempt: mellowmax operator + # # first blend: justRegion III power + # blend1 = f2 + + # second blend: smoothmin Region II and first blend coefficient + CTp_factor = -1.0/beta_smooth*ln( + exp(-beta_smooth) + exp(-beta_smooth*f2) + ) + + ### Compose full turbine force + self.tf = CTp_factor*(self.tf1*u[0]**2+self.tf2*u[1]**2+self.tf3*u[0]*u[1]) + + return self.tf + + def force_gradient(self): + pass + + def power(self, u, inflow_angle): + # adjust for turbine inefficiency + magic_number = 1.1291826687642552 # comes from calibration studies + return magic_number*self.calibration_factor_CPprime0*self.mCPprime0/(self.calibration_factor_CTprime0*self.mCTprime0)*dot(-self.tf,u)/1.0e6 # report in megawatts + + def thrust(self, u, inflow_angle): + # adjust for turbine inefficiency + nHat = u/sqrt(u[0]**2+u[1]**2+u[2]**2) + return dot(-self.tf,nHat)/1.0e3 # report in kilonewtons + + def power_gradient(self): + pass + + def prepare_saved_functions(self, func_list): + if len(func_list) == 0: + func_list = [ + [self.actuator_disk,"actuator_disk"], + [self.tf,"turbine_force"] + ] + else: + func_list[0][0] += self.actuator_disk + func_list[1][0] += self.tf + + return func_list diff --git a/windse/turbine_types/ActuatorDiskTSPowerCurve.py b/windse/turbine_types/ActuatorDiskTSPowerCurve.py new file mode 100644 index 00000000..1f74dbb8 --- /dev/null +++ b/windse/turbine_types/ActuatorDiskTSPowerCurve.py @@ -0,0 +1,213 @@ +# import windse function +from typing import Generic +from . import GenericTurbine + +# import dolfin functions +from . import Constant, SpatialCoordinate, as_vector, cos, sin, exp, ln, sqrt, dot, assemble, dx + +from . import ActuatorDisk + +import numpy as np +from scipy.special import gamma + +class ActuatorDiskTSPowerCurve(GenericTurbine): + + def __init__(self, i, x, y, dom, imported_params=None): + # Define the acceptable column of an wind_farm.csv imput file + self.yaml_inputs = ["HH", "RD", "yaw", "thickness", "CTprime0", "CPprime0", "Vrated", "force"] + + # Init turbine + super(ActuatorDiskTSPowerCurve, self).__init__(i,x,y,dom,imported_params) + + def load_parameters(self): + self.HH = self.params["turbines"]["HH"] + self.RD = self.params["turbines"]["RD"] + self.yaw = self.params["turbines"]["yaw"] + self.thickness = self.params["turbines"]["thickness"] + self.axial = -1 # needed for generic turbine + self.CPprime0 = self.params["turbines"]["CPprime0"] + self.CTprime0 = self.params["turbines"]["CTprime0"] + self.Prated = self.params["turbines"]["Prated"] + self.Trated = self.params["turbines"]["Trated"] + self.force = self.params["turbines"]["force"] + + def create_controls(self): + self.mx = Constant(self.x, name="x_{:d}".format(self.index)) + self.my = Constant(self.y, name="y_{:d}".format(self.index)) + self.myaw = Constant(self.yaw, name="yaw_{:d}".format(self.index)) + self.mCPprime0 = Constant(self.CPprime0, name="CPprime0_{:d}".format(self.index)) + self.mCTprime0 = Constant(self.CTprime0, name="CTprime0_{:d}".format(self.index)) + self.mPrated = Constant(self.Prated, name="Prated_{:d}".format(self.index)) + self.mTrated = Constant(self.Trated, name="Trated_{:d}".format(self.index)) + + # The control list is very important for extracting information from the optimzation step. + self.controls_list = ["x","y","yaw","CPprime0","CTprime0","Prated","Trated"] + + def build_actuator_disk(self): + + ### Alias useful values ### + x0 = [self.mx,self.my,self.mz] + yaw = self.myaw+self.inflow_angle + W = self.thickness*1.0 + R = self.RD/2.0 + CTprime0 = self.mCTprime0 + # CPprime0 = self.mCPprime0 + x=SpatialCoordinate(self.dom.mesh) + + ### Set up some dim dependent values ### + S_norm = (2.0+np.pi)/(2.0*np.pi) + T_norm = 2.0*gamma(7.0/6.0) + if self.dom.dim == 3: + WTGbase = as_vector((cos(yaw),sin(yaw),0.0)) + A = np.pi*R**2.0 + D_norm = np.pi*gamma(4.0/3.0) + else: + WTGbase = as_vector((cos(yaw),sin(yaw))) + A = 2*R + D_norm = 2.0*gamma(7.0/6.0) + + ### Rotate and Shift the Turbine ### + xs = self.yaw_turbine(x,x0,yaw) + + ### Create the function that represents the Thickness of the turbine ### + T = exp(-pow((xs[0]/W),6.0)) + + ### Create the function that represents the Disk of the turbine + r = sqrt(xs[1]**2.0+xs[2]**2.0)/R + D = exp(-pow(r,6.0)) + + ### Create the function that represents the force ### + if self.force == "constant": + force = 1.0 + elif self.force == "sine": + force = (r*sin(np.pi*r)+0.5)/S_norm + elif self.force == "chord": + raise NotImplementedError("not implemented! -cfrontin") + # self._setup_chord_force() + # chord = self.mchord[i] + # force = self.radial_chord_force(r,chord) + F = -0.5*A*CTprime0*force + + ### Calculate normalization constant ### + volNormalization = T_norm*D_norm*W*R**(self.dom.dim-1) + # volNormalization = assemble(T*D*dx) + + # self.fprint(f"Turbine {self.index}:") + # self.fprint(f"Assembled Volume: {assemble(T*D*dx)}",offset=1) + # self.fprint(f"Computed Volume: {volNormalization}",offset=1) + # self.fprint(f"Assembled Force: {assemble(F*T*D/volNormalization*dx)}",offset=1) + + # compute disk averaged velocity in yawed case and don't project + actuator_disks=F*T*D*WTGbase/volNormalization + # actuator_disks=project(actuator_disks,self.fs.V,solver_type='cg',preconditioner_type="hypre_amg") + + + return actuator_disks + + def turbine_force(self,u,inflow_angle,fs): + """ + This function creates a turbine force by applying + a spacial kernel to each turbine. This kernel is + created from the turbines location, yaw, thickness, diameter, + and force density. Currently, force density is limit to a scaled + version of + + .. math:: + + r\\sin(r), + + where :math:`r` is the distance from the center of the turbine. + + Args: + u: the velocity field + inflow_angle: direction the wind is blowing (0 is west to east and positive rotates clockwise) + fs: the function space manager object + + Returns: + tf (dolfin.Function): the turbine force. + + Todo: + * Setup a way to get the force density from file + """ + + ### compute the space kernal and radial force + self.fs = fs + self.inflow_angle.assign(inflow_angle) + self.actuator_disk = self.build_actuator_disk() + + ### Expand the dot product + yaw = self.myaw+self.inflow_angle + self.tf1 = self.actuator_disk * cos(yaw)**2 + self.tf2 = self.actuator_disk * sin(yaw)**2 + self.tf3 = self.actuator_disk * 2.0 * cos(yaw) * sin(yaw) + + # compute power curve adjustment as a function of velocity + CPprime0 = self.mCPprime0 + CTprime0 = self.mCTprime0 + Prated = self.mPrated + Trated = self.mTrated + Vrated = Prated*CTprime0/(Trated*CPprime0) + A = np.pi/4.0*self.RD**2.0 + Vps = sqrt(2*Prated/(A*Vrated*CPprime0)) + + # precompute key values + beta_smooth = 16.0 + vel_magnitude = sqrt(u[0]**2 + u[1]**2 + u[2]**2) + + f0 = (0.5*A*vel_magnitude**3)/1e6 + f1 = (0.5*A*vel_magnitude**3)*CTprime0/1e6 + f2 = (0.5*A*vel_magnitude**3)*CTprime0*Vps**2/(vel_magnitude**2)/1e6 + f3 = (0.5*A*vel_magnitude**3)*CTprime0*Vps**2*Vrated/(vel_magnitude**3)/1e6 + + # smooth twice: one in power space and one in coefficient space + + # # first attempt: boltzmann operator + # # first blend: smoothmin Region II.5 and Region III power + # blend1 = 1.0*( + # f2*exp(-beta_smooth*f2) + f3*exp(-beta_smooth*f3) + # )/( + # exp(-beta_smooth*f2) + exp(-beta_smooth*f3) + # ) + # + # # second blend: smoothmin Region II and first blend coefficient + # CTp_factor = ( + # f1/f0*exp(-beta_smooth*f1/f0) + blend1/f0*exp(-beta_smooth*blend1/f0) + # )/( + # exp(-beta_smooth*f1/f0) + exp(-beta_smooth*blend1/f0) + # ) + + # second attempt: mellowmax operator + # first blend: smoothmin Region II.5 and Region III power + blend1 = -1.0/beta_smooth*ln(exp(-beta_smooth*f2) + exp(-beta_smooth*f3)) + + # second blend: smoothmin Region II and first blend coefficient + CTp_factor = -1.0/beta_smooth*ln( + exp(-beta_smooth*f1/f0) + exp(-beta_smooth*blend1/f0) + )/CTprime0 + + ### Compose full turbine force + self.tf = CTp_factor*(self.tf1*u[0]**2+self.tf2*u[1]**2+self.tf3*u[0]*u[1]) + + return self.tf + + def force_gradient(self): + pass + + def power(self, u, inflow_angle): + # adjust for turbine inefficiency + return self.mCPprime0/self.mCTprime0*dot(-self.tf,u)/1.0e6 + + def power_gradient(self): + pass + + def prepare_saved_functions(self, func_list): + if len(func_list) == 0: + func_list = [ + [self.actuator_disk,"actuator_disk"], + [self.tf,"turbine_force"] + ] + else: + func_list[0][0] += self.actuator_disk + func_list[1][0] += self.tf + + return func_list diff --git a/windse/turbine_types/__init__.py b/windse/turbine_types/__init__.py index 81314df6..ad618f5e 100644 --- a/windse/turbine_types/__init__.py +++ b/windse/turbine_types/__init__.py @@ -11,6 +11,8 @@ from .ActuatorDisk2D import ActuatorDisk2D from .ActuatorDiskExpr import ActuatorDiskExpr from .ActuatorDiskNumpy import ActuatorDiskNumpy +from .ActuatorDiskSimplePowerCurve import ActuatorDiskSimplePowerCurve +from .ActuatorDiskTSPowerCurve import ActuatorDiskTSPowerCurve from .ActuatorHybridDisk import ActuatorHybridDisk from .ActuatorLine import ActuatorLine from .ActuatorLineDolfin import ActuatorLineDolfin @@ -28,6 +30,10 @@ "numpy_disks": ActuatorDiskNumpy, "hybrid_disk": ActuatorHybridDisk, "hybrid_disks": ActuatorHybridDisk, + "power_disk": ActuatorDiskSimplePowerCurve, + "power_disks": ActuatorDiskSimplePowerCurve, + "power_ts_disk": ActuatorDiskTSPowerCurve, + "power_ts_disks": ActuatorDiskTSPowerCurve, "line": ActuatorLine, "lines": ActuatorLine, "dolfin_line": ActuatorLineDolfin, diff --git a/windse/wind_farm_types/GenericWindFarm.py b/windse/wind_farm_types/GenericWindFarm.py index 26d5e6d7..20a81999 100644 --- a/windse/wind_farm_types/GenericWindFarm.py +++ b/windse/wind_farm_types/GenericWindFarm.py @@ -4,14 +4,14 @@ import time, os from . import MeshFunction, CompiledSubDomain, Measure, cells, project, inner, FiniteElement, FunctionSpace, MixedElement, assemble, dx, parameters, Form, File import matplotlib.pyplot as plt -from pyadjoint.tape import stop_annotating +from pyadjoint.tape import stop_annotating from pyadjoint import AdjFloat class GenericWindFarm(object): """ A GenericWindFarm contains on the basic functions and attributes required by all wind farm objects. - - Args: + + Args: dom (:meth:`windse.DomainManager.GenericDomain`): a windse domain object. """ def __init__(self): @@ -34,6 +34,7 @@ def __init__(self): # Set any global flags self.func_first_save = True self.power_first_save = True + self.thrust_first_save = True # Init blank turbine list self.turbines = [] @@ -45,7 +46,7 @@ def __init__(self): self.control_types = self.params["optimization"]["control_types"] self.optimizing = True - ### extra_kwargs will reduce overhead on operations that we don't want dolfin_adjoint to track + ### extra_kwargs will reduce overhead on operations that we don't want dolfin_adjoint to track self.extra_kwarg = {} if self.params.dolfin_adjoint: self.extra_kwarg["annotate"] = False @@ -70,23 +71,29 @@ def Finalize(self, dom): def setup(self): """ This function builds the wind farm as well as sets up the turbines - """ - self.load_parameters() - self.compute_parameters() + """ + self.load_parameters() + self.compute_parameters() self.fprint("Number of Turbines: {:d}".format(self.numturbs)) self.fprint("Type of Turbines: {}".format(self.turbine_type)) self.initial_turbine_locations = self.initialize_turbine_locations() + self.fprint("Setting Up Turbines: ",special="header") + self.setup_turbines() + self.fprint("Turbines Set up",special="footer") + + self.debug_output() + def load_parameters(self): """ This function will parse the parameters from the yaml file - """ + """ raise NotImplementedError(type(self)) def compute_parameters(self): """ This function will compute any additional parameters - """ + """ raise NotImplementedError(type(self)) def initialize_turbine_locations(): @@ -142,7 +149,7 @@ def get_yaw_angles(self): def calculate_farm_bounding_box(self): """ - This functions takes into consideration the turbine locations, diameters, + This functions takes into consideration the turbine locations, diameters, and hub heights to create lists that describe the extent of the windfarm. These lists are append to the parameters object. """ @@ -150,7 +157,7 @@ def calculate_farm_bounding_box(self): RD = self.get_rotor_diameters() ground = self.get_ground_heights() - ### Locate the extreme turbines ### + ### Locate the extreme turbines ### x_min = np.argmin(x) x_max = np.argmax(x) y_min = np.argmin(y) @@ -162,12 +169,12 @@ def calculate_farm_bounding_box(self): self.ex_x = [x[x_min]-RD[x_min]/2.0,x[x_max]+RD[x_max]/2.0] self.ex_y = [y[y_min]-RD[y_min]/2.0,y[y_max]+RD[y_max]/2.0] self.ex_z = [min(ground),z[z_max]+RD[z_max]/2.0] - + return [self.ex_x,self.ex_y,self.ex_z] def update_heights(self): """ - updates the hub height for all turbines in farm + updates the hub height for all turbines in farm """ for turb in self.turbines: turb.calculate_heights() @@ -205,7 +212,7 @@ def compute_turbine_force(self,u,v,inflow_angle,fs,**kwargs): # create a measure for the turbines self.local_dx = Measure('dx', subdomain_data=self.turbine_subdomains) - # check if the number of turbines might trigger a python recursion error. + # check if the number of turbines might trigger a python recursion error. if self.numturbs > self.too_many_turbines: final_tf_list, final_local_dx = self.compress_turbine_function(u,fs) else: @@ -239,7 +246,7 @@ def compress_turbine_function(self, u, fs): final_local_dx = [] counter = 0 group = 1 - + # Create a cell function that will be used to store the subdomains temp_subdomain = MeshFunction("size_t", self.dom.mesh, self.dom.mesh.topology().dim()) temp_subdomain.set_all(0) @@ -265,7 +272,7 @@ def compress_turbine_function(self, u, fs): if self.dom.dim == 3: subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)+(x[2]-HH)*(x[2]-HH)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, HH=turb.z, x0=turb.x, y0=turb.y) else: - subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) + subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) # add the subdomain to the marker subdomain_marker.mark(temp_subdomain,group) @@ -279,7 +286,7 @@ def compress_turbine_function(self, u, fs): self.fprint(f"Projecting tf3 for group {group}") proj_tf3 = project(tf3,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") combined_tf_list.append(proj_tf1*u[0]**2+proj_tf2*u[1]**2+proj_tf3*u[0]*u[1]) - + # reset loop values tf1 = 0 tf2 = 0 @@ -319,7 +326,7 @@ def compress_turbine_function(self, u, fs): subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)+(x[2]-HH)*(x[2]-HH)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, HH=turb.z, x0=turb.x, y0=turb.y) else: subdomain_marker = CompiledSubDomain("(x[0]-x0)*(x[0]-x0)+(x[1]-y0)*(x[1]-y0)<=r*r",r=self.local_dx_scaling*turb.RD/2.0, x0=turb.x, y0=turb.y) - + # add the subdomain to the marker subdomain_marker.mark(temp_subdomain,group) @@ -328,7 +335,7 @@ def compress_turbine_function(self, u, fs): self.fprint(f"Projecting tf for group {group}") proj_tf = project(tf,fs.V,solver_type='gmres',preconditioner_type="hypre_amg") combined_tf_list.append(proj_tf) - + # reset loop values tf = 0 counter = 0 @@ -346,12 +353,12 @@ def compress_turbine_function(self, u, fs): # create a measure for the final turbine group final_local_dx = Measure('dx', subdomain_data=temp_subdomain) - # create the final turbine force + # create the final turbine force return combined_tf_list, final_local_dx def update_controls(self): """ - iterates over the controls list assigns self. to the control, self.m + iterates over the controls list assigns self. to the control, self.m """ for turb in self.turbines: turb.update_controls() @@ -363,13 +370,38 @@ def update_turbine_force(self, u, inflow_angle, fs, **kwargs): for turb in self.turbines: turb.update_turbine_force(u, inflow_angle, fs, **kwargs) + def compute_thrust(self, u, inflow_angle): + """ + Computes the per-turbine net thrust for the full farm + """ + + ### Build the integrand + val = 0 + for turb in self.turbines: + temp = turb.thrust(u, inflow_angle) + if not isinstance(temp,(float,AdjFloat)): + if self.use_local_tf_dx: + val += temp*self.local_dx(turb.index+1) + else: + val += temp*dx + else: + val += temp + + ### Assemble if needed + if not isinstance(val,(float,AdjFloat)): + J = assemble(val) + else: + J = val + + return J + def compute_power(self, u, inflow_angle): """ Computes the power for the full farm """ ### Build the integrand - val = 0 + val = 0.0 for turb in self.turbines: temp = turb.power(u, inflow_angle) if not isinstance(temp,(float,AdjFloat)): @@ -388,6 +420,41 @@ def compute_power(self, u, inflow_angle): return J + def save_thrust(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): + """ + saves the thrust output for each turbine + """ + + J_list=np.zeros(self.numturbs+3) + J_list[0]=iter_val + J_list[1]=simTime + with stop_annotating(): + for i,turb in enumerate(self.turbines): + val = turb.thrust(u,inflow_angle) + + ### Assemble if needed + if not isinstance(val,(float,AdjFloat)): + if self.use_local_tf_dx: + J = assemble(val*self.local_dx(turb.index+1)) + else: + J = assemble(val*dx) + else: + J = val + + J_list[i+2] = J + + J_list[-1]=sum(J_list[2:]) + + folder_string = self.params.folder+"data/" + if not os.path.exists(folder_string): os.makedirs(folder_string) + + if self.thrust_first_save: + header = str("Iter_val "+"Time "+"Turbine_%d "*self.numturbs % tuple(range(self.numturbs))+"Sum") + self.params.save_csv("thrust_data",header=header,data=[J_list],subfolder=self.params.folder+"data/",mode='w') + self.thrust_first_save = False + else: + self.params.save_csv("thrust_data",data=[J_list],subfolder=self.params.folder+"data/",mode='a') + def save_power(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): """ saves the power output for each turbine @@ -408,7 +475,7 @@ def save_power(self, u, inflow_angle, iter_val = 0.0, simTime = 0.0): J = assemble(val*dx) else: J = val - + J_list[i+2] = J J_list[-1]=sum(J_list[2:]) @@ -576,10 +643,10 @@ def plot_chord(self,show=False,filename="chord_profiles",objective_value=None, b elif isinstance(objective_value,(list,np.ndarray)): plt.title("Objective Value: {: 5.6f}".format(sum(objective_value))) else: - plt.title("Objective Value: {: 5.6f}".format(float(objective_value))) - plt.xlabel("Blade Span") + plt.title("Objective Value: {: 5.6f}".format(float(objective_value))) + plt.xlabel("Blade Span") plt.ylabel("Chord") - plt.legend() + plt.legend() plt.savefig(file_string, transparent=True) @@ -591,7 +658,7 @@ def plot_chord(self,show=False,filename="chord_profiles",objective_value=None, b def save_functions(self,val=0): """ This function call the prepare_saved_functions from each turbine, combines the functions and saves them. - It then check to see if it can save the function out right and if not it projects. + It then check to see if it can save the function out right and if not it projects. "val" can be the time, or angle, its just the iterator for saving mulitple steps Note: this function is way over engineered! """ @@ -621,7 +688,7 @@ def save_functions(self,val=0): # Note: this is overkill if not hasattr(func,"_cpp_object"): - # choose an appropriate function space + # choose an appropriate function space if func.geometric_dimension() == 1: FS = self.fs.Q else: @@ -639,7 +706,7 @@ def save_functions(self,val=0): self.func_files.append(out_file) else: self.params.Save(func, func_name,subfolder="functions/",val=val,file=self.func_files[i]) - + # save the farm level turbine subdomain function if self.numturbs > 0: if self.func_first_save: @@ -717,7 +784,7 @@ def finalize_farm(self): ######################################################################################################## ############## It would be smart to eventually move these to a separate refinement module ############## ######################################################################################################## - + def SimpleRefine(self,radius,expand_factor=1): if self.numturbs == 0: return @@ -730,7 +797,7 @@ def SimpleRefine(self,radius,expand_factor=1): ### Create the cell markers ### cell_f = MeshFunction('bool', self.dom.mesh, self.dom.mesh.geometry().dim(),False) - + ### Get Dimension ### n = self.numturbs d = self.dom.dim @@ -882,7 +949,7 @@ def TearRefine(self,radius,theta=0.0,expand_factor=1): ### Create the cell markers ### cell_f = MeshFunction('bool', self.dom.mesh, self.dom.mesh.geometry().dim(),False) - + ### Get Dimension ### n = self.numturbs d = self.dom.dim @@ -981,7 +1048,7 @@ def SphereRefine(self,radius,expand_factor=1): ### Create the cell markers ### cell_f = MeshFunction('bool', self.dom.mesh, self.dom.mesh.geometry().dim(),False) - + ### Get Dimension ### n = self.numturbs d = self.dom.dim @@ -992,7 +1059,7 @@ def SphereRefine(self,radius,expand_factor=1): turb_y = turb_locs[:,1] if self.dom.dim == 3: turb_z = turb_locs[:,2] - + self.fprint("Marking Near Turbine") mark_start = time.time() @@ -1030,4 +1097,3 @@ def SphereRefine(self,radius,expand_factor=1): refine_stop = time.time() self.fprint("Mesh Refinement Finished: {:1.2f} s".format(refine_stop-refine_start),special="footer") - \ No newline at end of file