Enable cc model plotting (#335)

* adding full flow solver for CC model

* enabling switching between full flow solvers for plotting

* formatting

src/floris/simulation/__init__.py

@@ -39,7 +39,7 @@
 from .grid import Grid, TurbineGrid, FlowFieldGrid, FlowFieldPlanarGrid
 from .flow_field import FlowField
 from .wake import WakeModelManager
-from .solver import sequential_solver, full_flow_sequential_solver, cc_solver
+from .solver import sequential_solver, full_flow_sequential_solver, cc_solver, full_flow_cc_solver
 from .floris import Floris
 
 

src/floris/simulation/floris.py

@@ -31,7 +31,8 @@
     FlowFieldPlanarGrid,
     sequential_solver,
     cc_solver,
-    full_flow_sequential_solver
+    full_flow_sequential_solver,
+    full_flow_cc_solver
 )
 from attrs import define, field
 
@@ -156,11 +157,14 @@ def solve_for_viz(self):
         # This function call should be for a single wind direction and wind speed
         # since the memory consumption is very large.
 
-        # self.steady_state_atmospheric_condition()
-
         self.flow_field.initialize_velocity_field(self.grid)
 
-        full_flow_sequential_solver(self.farm, self.flow_field, self.turbine, self.grid, self.wake)
+        vel_model = self.wake.model_strings["velocity_model"]
+
+        if vel_model=="cc":
+            full_flow_cc_solver(self.farm, self.flow_field, self.turbine, self.grid, self.wake)
+        else:
+            full_flow_sequential_solver(self.farm, self.flow_field, self.turbine, self.grid, self.wake)
 
 
     ## I/O

src/floris/simulation/solver.py

@@ -496,3 +496,131 @@ def cc_solver(farm: Farm, flow_field: FlowField, turbine: Turbine, grid: Turbine
 
     flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3,4))
     flow_field.turbulence_intensity_field = flow_field.turbulence_intensity_field[:,:,:,None,None]
+
+
+def full_flow_cc_solver(farm: Farm, flow_field: FlowField, turbine: Turbine, flow_field_grid: FlowFieldGrid, model_manager: WakeModelManager) -> None:
+
+    # Get the flow quantities and turbine performance
+    turbine_grid_farm = copy.deepcopy(farm)
+    turbine_grid_flow_field = copy.deepcopy(flow_field)
+    turbine_grid = TurbineGrid(
+        turbine_coordinates=turbine_grid_farm.coordinates,
+        reference_turbine_diameter=turbine.rotor_diameter,
+        wind_directions=turbine_grid_flow_field.wind_directions,
+        wind_speeds=turbine_grid_flow_field.wind_speeds,
+        grid_resolution=3,
+    )
+    turbine_grid_flow_field.initialize_velocity_field(turbine_grid)
+    cc_solver(turbine_grid_farm, turbine_grid_flow_field, turbine, turbine_grid, model_manager)
+
+    ### Referring to the quantities from above, calculate the wake in the full grid
+
+    # Use full flow_field here to use the full grid in the wake models
+    deflection_model_args = model_manager.deflection_model.prepare_function(flow_field_grid, flow_field, turbine)
+    deficit_model_args = model_manager.velocity_model.prepare_function(flow_field_grid, flow_field, turbine)
+
+    v_wake = np.zeros_like(flow_field.v_initial)
+    w_wake = np.zeros_like(flow_field.w_initial)
+    turb_u_wake = np.zeros_like(flow_field.u_initial)
+
+    shape = (farm.n_turbines,) + np.shape(flow_field.u_initial)
+    Ctmp = np.zeros((shape))
+
+    # Calculate the velocity deficit sequentially from upstream to downstream turbines
+    for i in range(flow_field_grid.n_turbines):
+
+        # Get the current turbine quantities
+        x_i = np.mean(turbine_grid.x[:, :, i:i+1], axis=(3, 4))
+        x_i = x_i[:, :, :, None, None]
+        y_i = np.mean(turbine_grid.y[:, :, i:i+1], axis=(3, 4))        
+        y_i = y_i[:, :, :, None, None]
+        z_i = np.mean(turbine_grid.z[:, :, i:i+1], axis=(3, 4))
+        z_i = z_i[:, :, :, None, None]
+
+        u_i = turbine_grid_flow_field.u[:, :, i:i+1]
+        v_i = turbine_grid_flow_field.v[:, :, i:i+1]
+
+        turb_avg_vels = average_velocity(turbine_grid_flow_field.u)
+        turb_Cts = Ct(
+            velocities=turb_avg_vels,
+            yaw_angle=turbine_grid_farm.yaw_angles,
+            fCt=turbine.fCt_interp,
+        )
+        turb_Cts = turb_Cts[:, :, :, None, None]
+
+        axial_induction_i = axial_induction(
+            velocities=turbine_grid_flow_field.u,
+            yaw_angle=turbine_grid_farm.yaw_angles,
+            fCt=turbine.fCt_interp,
+            ix_filter=[i],
+        )
+        axial_induction_i = axial_induction_i[:, :, :, None, None]
+
+        turbulence_intensity_i = turbine_grid_flow_field.turbulence_intensity_field[:, :, i:i+1]
+        yaw_angle_i = turbine_grid_farm.yaw_angles[:, :, i:i+1, None, None]
+
+        effective_yaw_i = np.zeros_like(yaw_angle_i)
+        effective_yaw_i += yaw_angle_i
+
+        if model_manager.enable_secondary_steering:
+            added_yaw = wake_added_yaw(
+                u_i,
+                v_i,
+                turbine_grid_flow_field.u_initial,
+                turbine_grid.y[:, :, i:i+1] - y_i,
+                turbine_grid.z[:, :, i:i+1],
+                turbine.rotor_diameter,
+                turbine.hub_height,
+                turb_Cts[:, :, i:i+1],
+                turbine.TSR,
+                axial_induction_i,
+                scale=2.0
+            )
+            effective_yaw_i += added_yaw
+
+        # Model calculations
+        # NOTE: exponential
+        deflection_field = model_manager.deflection_model.function(
+            x_i,
+            y_i,
+            effective_yaw_i,
+            turbulence_intensity_i,
+            turb_Cts[:, :, i:i+1],
+            **deflection_model_args
+        )
+
+        if model_manager.enable_transverse_velocities:
+            v_wake, w_wake = calculate_transverse_velocity(
+                u_i,
+                flow_field.u_initial,
+                flow_field_grid.x - x_i,
+                flow_field_grid.y - y_i,
+                flow_field_grid.z,
+                turbine.rotor_diameter,
+                turbine.hub_height,
+                yaw_angle_i,
+                turb_Cts[:, :, i:i+1],
+                turbine.TSR,
+                axial_induction_i,
+                scale=2.0
+            )
+
+        # NOTE: exponential
+        turb_u_wake, Ctmp = model_manager.velocity_model.function(
+            i,
+            x_i,
+            y_i,
+            z_i,
+            u_i,
+            deflection_field,
+            yaw_angle_i,
+            turbine_grid_flow_field.turbulence_intensity_field,
+            turb_Cts,
+            turb_u_wake,
+            Ctmp,
+            **deficit_model_args
+        )
+
+        flow_field.v += v_wake
+        flow_field.w += w_wake
+    flow_field.u = flow_field.u_initial - turb_u_wake