Support multidimensional turbine definitions in all wake models (#812)

* Allow cc solver to accept multidim conditions.

* Emgauss runs with mutlidim conditions.

* Remove multidim_cp_ct option.

* Remove references to multidim_cp_ct wake model.

* Add multidim capabilities to full_flow solvers.

* Turbopark solver compatibility.

* add parameters for other models to multidim input yaml for exposition.

* Add regression tests for full flow solvers

* CC full flow solver bug fix

---------

Co-authored-by: Rafael M Mudafort <rafmudaf@gmail.com>

examples/30_multi_dimensional_cp_ct.py

@@ -19,8 +19,7 @@
 height. For every combination of Tp and Hs defined, a Cp/Ct/Wind speed table of values is also
 defined. It is required for this .csv file to have the last 3 columns be ws, Cp, and Ct. In order
 for this table to be used, the flag 'multi_dimensional_cp_ct' must be present and set to true in
-the turbine definition. Also of note is the 'velocity_model' must be set to 'multidim_cp_ct' in
-the main input file. With both of these values provided, the solver will downselect to use the
+the turbine definition. With this flag enabled, the solver will downselect to use the
 interpolant defined at the closest conditions. The user must supply these conditions in the
 main input file under the 'flow_field' section, e.g.:
 

examples/inputs/gch.yaml

@@ -139,7 +139,7 @@ flow_field:
   # The conditions that are specified for use with the multi-dimensional Cp/Ct capbility.
   # These conditions are external to FLORIS and specified by the user. They are used internally
   # through a nearest-neighbor selection process to choose the correct Cp/Ct interpolants
-  # to use. These conditions are only used with the ``multidim_cp_ct`` velocity deficit model.
+  # to use.
   multidim_conditions:
     Tp: 2.5
     Hs: 3.01

examples/inputs/gch_multi_dim_cp_ct.yaml

@@ -47,7 +47,7 @@ wake:
     combination_model: sosfs
     deflection_model: gauss
     turbulence_model: crespo_hernandez
-    velocity_model: multidim_cp_ct
+    velocity_model: gauss
 
   enable_secondary_steering: true
   enable_yaw_added_recovery: true
@@ -66,6 +66,12 @@ wake:
       ad: 0.0
       bd: 0.0
       kd: 0.05
+    empirical_gauss:
+      horizontal_deflection_gain_D: 3.0
+      vertical_deflection_gain_D: -1
+      deflection_rate: 30
+      mixing_gain_deflection: 0.0
+      yaw_added_mixing_gain: 0.0
 
   wake_velocity_parameters:
     cc:
@@ -77,17 +83,31 @@ wake:
       b_f: -0.68
       c_f: 2.41
       alpha_mod: 1.0
-    multidim_cp_ct:
+    gauss:
       alpha: 0.58
       beta: 0.077
       ka: 0.38
       kb: 0.004
     jensen:
       we: 0.05
+    turbopark:
+      A: 0.04
+      sigma_max_rel: 4.0
+    empirical_gauss:
+      wake_expansion_rates:
+      - 0.023
+      - 0.008
+      breakpoints_D:
+      - 10
+      sigma_0_D: 0.28
+      smoothing_length_D: 2.0
+      mixing_gain_velocity: 2.0
 
   wake_turbulence_parameters:
     crespo_hernandez:
       initial: 0.01
       constant: 0.9
       ai: 0.83
       downstream: -0.25
+    wake_induced_mixing:
+      atmospheric_ti_gain: 0.0

floris/simulation/solver.py

@@ -333,6 +333,9 @@ def full_flow_sequential_solver(
             turbine_type_map=turbine_grid_farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             ix_filter=[i],
+            average_method=turbine_grid.average_method,
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the thrust_coefficient function,
         # get the first index here (0:1)
@@ -349,6 +352,9 @@ def full_flow_sequential_solver(
             turbine_type_map=turbine_grid_farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             ix_filter=[i],
+            average_method=turbine_grid.average_method,
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the axial induction function,
         # get the first index here (0:1)
@@ -502,7 +508,8 @@ def cc_solver(
             turbine_type_map=farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         turb_Cts = turb_Cts[:, :, None, None]
         turb_aIs = axial_induction(
@@ -518,7 +525,8 @@ def cc_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         turb_aIs = turb_aIs[:, :, None, None]
 
@@ -538,7 +546,8 @@ def cc_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
 
         axial_induction_i = axial_induction_i[:, :, None, None]
@@ -670,7 +679,7 @@ def cc_solver(
 def full_flow_cc_solver(
     farm: Farm,
     flow_field: FlowField,
-    flow_field_grid: FlowFieldGrid,
+    flow_field_grid: FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid,
     model_manager: WakeModelManager,
 ) -> None:
     # Get the flow quantities and turbine performance
@@ -740,7 +749,7 @@ def full_flow_cc_solver(
         turb_avg_vels = average_velocity(turbine_grid_flow_field.u_sorted)
         turb_Cts = thrust_coefficient(
             velocities=turb_avg_vels,
-            air_density=flow_field_grid.air_density,
+            air_density=turbine_grid_flow_field.air_density,
             yaw_angles=turbine_grid_farm.yaw_angles_sorted,
             tilt_angles=turbine_grid_farm.tilt_angles_sorted,
             power_setpoints=turbine_grid_farm.power_setpoints_sorted,
@@ -750,7 +759,8 @@ def full_flow_cc_solver(
             turbine_type_map=turbine_grid_farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             average_method=turbine_grid.average_method,
-            cubature_weights=turbine_grid.cubature_weights
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         turb_Cts = turb_Cts[:, :, None, None]
 
@@ -767,7 +777,8 @@ def full_flow_cc_solver(
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=turbine_grid.average_method,
-            cubature_weights=turbine_grid.cubature_weights
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         axial_induction_i = axial_induction_i[:, :, None, None]
 
@@ -905,7 +916,8 @@ def turbopark_solver(
             turbine_type_map=farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
 
         ct_i = thrust_coefficient(
@@ -921,7 +933,8 @@ def turbopark_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the thrust coefficient function,
         # get the first index here (0:1)
@@ -939,7 +952,8 @@ def turbopark_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the axial induction function,
         # get the first index here (0:1)
@@ -984,7 +998,8 @@ def turbopark_solver(
                     turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
                     ix_filter=[ii],
                     average_method=grid.average_method,
-                    cubature_weights=grid.cubature_weights
+                    cubature_weights=grid.cubature_weights,
+                    multidim_condition=flow_field.multidim_conditions,
                 )
                 ct_ii = ct_ii[:, 0:1, None, None]
                 rotor_diameter_ii = farm.rotor_diameters_sorted[:, ii:ii+1, None, None]
@@ -1158,7 +1173,8 @@ def empirical_gauss_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the thrust coefficient function,
         # get the first index here (0:1)
@@ -1176,7 +1192,8 @@ def empirical_gauss_solver(
             turbine_power_thrust_tables=farm.turbine_power_thrust_tables,
             ix_filter=[i],
             average_method=grid.average_method,
-            cubature_weights=grid.cubature_weights
+            cubature_weights=grid.cubature_weights,
+            multidim_condition=flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the axial induction function,
         # get the first index here (0:1)
@@ -1359,6 +1376,9 @@ def full_flow_empirical_gauss_solver(
             turbine_type_map=turbine_grid_farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             ix_filter=[i],
+            average_method=turbine_grid.average_method,
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the thrust coefficient function,
         # get the first index here (0:1)
@@ -1375,6 +1395,9 @@ def full_flow_empirical_gauss_solver(
             turbine_type_map=turbine_grid_farm.turbine_type_map_sorted,
             turbine_power_thrust_tables=turbine_grid_farm.turbine_power_thrust_tables,
             ix_filter=[i],
+            average_method=turbine_grid.average_method,
+            cubature_weights=turbine_grid.cubature_weights,
+            multidim_condition=turbine_grid_flow_field.multidim_conditions,
         )
         # Since we are filtering for the i'th turbine in the axial induction function,
         # get the first index here (0:1)

floris/simulation/wake.py

@@ -53,7 +53,6 @@
         "jensen": JensenVelocityDeficit,
         "turbopark": TurbOParkVelocityDeficit,
         "empirical_gauss": EmpiricalGaussVelocityDeficit,
-        "multidim_cp_ct": GaussVelocityDeficit
     },
 }
 

tests/reg_tests/cumulative_curl_regression_test.py

@@ -137,6 +137,48 @@
     ]
 )
 
+full_flow_baseline = np.array(
+    [
+        [
+            [
+                [7.88772361, 8.00000000, 8.10178821],
+                [7.88772361, 8.00000000, 8.10178821],
+                [7.88772361, 8.00000000, 8.10178821],
+                [7.88772361, 8.00000000, 8.10178821],
+                [7.88772361, 8.00000000, 8.10178821],
+            ],
+            [
+                [7.88772361, 8.        , 8.10178821],
+                [7.85396979, 7.96487892, 8.06803439],
+                [4.19559099, 4.28925565, 4.40965558],
+                [7.85396979, 7.96487892, 8.06803439],
+                [7.88772361, 8.        , 8.10178821],
+            ],
+            [
+                [7.88769642, 7.99997223, 8.10176102],
+                [7.58415314, 7.69072103, 7.79821773],
+                [4.16725762, 4.26342392, 4.38132221],
+                [7.58415314, 7.69072103, 7.79821773],
+                [7.88769642, 7.99997223, 8.10176102],
+            ],
+            [
+                [7.88513176, 7.99737618, 8.09919636],
+                [7.21888868, 7.32333558, 7.43301511],
+                [4.30201226, 4.40270245, 4.51689213],
+                [7.21888868, 7.32333558, 7.43301511],
+                [7.88513176, 7.99737618, 8.09919636],
+            ],
+            [
+                [7.86539121, 7.97748824, 8.0794561 ],
+                [7.0723371 , 7.1790733 , 7.28645574],
+                [5.8436738 , 5.95178931, 6.05791862],
+                [7.0723371 , 7.1790733 , 7.28645574],
+                [7.86539121, 7.97748824, 8.0794561 ],
+            ]
+        ]
+    ]
+)
+
 
 # Note: compare the yawed vs non-yawed results. The upstream turbine
 # power should be lower in the yawed case. The following turbine
@@ -624,3 +666,32 @@ def test_regression_small_grid_rotation(sample_inputs_fixture):
     assert np.allclose(farm_powers[8,0:5], farm_powers[8,15:20])
     assert np.allclose(farm_powers[8,20], farm_powers[8,0])
     assert np.allclose(farm_powers[8,21], farm_powers[8,21:25])
+
+
+def test_full_flow_solver(sample_inputs_fixture):
+    """
+    Full flow solver test with the flow field planar grid.
+    This requires one wind condition, and the grid is deliberately coarse to allow for
+    visually comparing results, as needed.
+    The u-component of velocity is compared, and the array has the shape
+    (n_findex, n_turbines, n grid points in x, n grid points in y, 3 grid points in z).
+    """
+
+    sample_inputs_fixture.floris["wake"]["model_strings"]["velocity_model"] = VELOCITY_MODEL
+    sample_inputs_fixture.floris["wake"]["model_strings"]["deflection_model"] = DEFLECTION_MODEL
+    sample_inputs_fixture.floris["solver"] = {
+        "type": "flow_field_planar_grid",
+        "normal_vector": "z",
+        "planar_coordinate": sample_inputs_fixture.floris["farm"]["turbine_type"][0]["hub_height"],
+        "flow_field_grid_points": [5, 5],
+        "flow_field_bounds": [None, None],
+    }
+    sample_inputs_fixture.floris["flow_field"]["wind_directions"] = [270.0]
+    sample_inputs_fixture.floris["flow_field"]["wind_speeds"] = [8.0]
+
+    floris = Floris.from_dict(sample_inputs_fixture.floris)
+    floris.solve_for_viz()
+
+    velocities = floris.flow_field.u_sorted
+
+    assert_results_arrays(velocities, full_flow_baseline)