Skip to content

Commit

Permalink
Save cubature weights on Grids
Browse files Browse the repository at this point in the history
Rather than the full cubature coefficient data
  • Loading branch information
rafmudaf committed May 11, 2023
1 parent cc406de commit 92ab0df
Show file tree
Hide file tree
Showing 5 changed files with 34 additions and 35 deletions.
32 changes: 16 additions & 16 deletions floris/simulation/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ class Grid(ABC):
x_sorted: NDArrayFloat = field(init=False)
y_sorted: NDArrayFloat = field(init=False)
z_sorted: NDArrayFloat = field(init=False)
cubature_coefficients: dict = field(init=False, default=None)
cubature_weights: NDArrayFloat = field(init=False, default=None)

def __attrs_post_init__(self) -> None:
self.turbine_coordinates_array = np.array([c.elements for c in self.turbine_coordinates])
Expand Down Expand Up @@ -304,15 +304,17 @@ def set_grid(self) -> None:
)

# Coefficients
self.cubature_coefficients = CubatureGrid.get_cubature_coefficients(self.grid_resolution)
cubature_coefficients = CubatureGrid.get_cubature_coefficients(self.grid_resolution)

# Generate grid points
yv = np.kron(self.cubature_coefficients["r"], self.cubature_coefficients["q"])
zv = np.kron(self.cubature_coefficients["r"], self.cubature_coefficients["t"])
yv = np.kron(cubature_coefficients["r"], cubature_coefficients["q"])
zv = np.kron(cubature_coefficients["r"], cubature_coefficients["t"])

# Calculate weighing terms for the grid points
self.cubature_coefficients["CUBcoeff"] = \
np.kron(self.cubature_coefficients["A"],np.ones((1, N))) * self.cubature_coefficients["B"] / np.pi
# Calculate weighting terms for the grid points
self.cubature_weights = (
np.kron(cubature_coefficients["A"], np.ones((1, self.grid_resolution)))
* cubature_coefficients["B"] / np.pi
)

# Here, the coordinates are already rotated to the correct orientation for each
# wind direction
Expand Down Expand Up @@ -424,15 +426,13 @@ def get_cubature_coefficients(cls, N: int):
q = [ 0.1564344650402308690101053, 0.4539904997395467915604084, 0.7071067811865475244008444, 0.8910065241883678623597096, 0.9876883405951377261900402, 0.9876883405951377261900402, 0.8910065241883678623597096, 0.7071067811865475244008444, 0.4539904997395467915604084, 0.1564344650402308690101053] # noqa: E501
A = [ 0.0592317212640472718785660, 0.1196571676248416170103229, 0.1422222222222222222222222, 0.1196571676248416170103229, 0.0592317212640472718785660, 0.0592317212640472718785660, 0.1196571676248416170103229, 0.1422222222222222222222222, 0.1196571676248416170103229, 0.0592317212640472718785660] # noqa: E501

cubature_coefficients = {
"r": np.array(r, dtype=float),
"t": np.array(t, dtype=float),
"q": np.array(q, dtype=float),
"A": np.array(A, dtype=float),
"B": np.pi/N,
}

return cubature_coefficients
return {
"r": np.array(r, dtype=float),
"t": np.array(t, dtype=float),
"q": np.array(q, dtype=float),
"A": np.array(A, dtype=float),
"B": np.pi/N,
}

@define
class FlowFieldGrid(Grid):
Expand Down
14 changes: 7 additions & 7 deletions floris/simulation/solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ def sequential_solver(
average_velocities = average_velocity(
flow_field.u_sorted,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)

ct_i = Ct(
Expand Down Expand Up @@ -491,7 +491,7 @@ def cc_solver(
turb_avg_vels = average_velocity(
turb_inflow_field,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)
turb_Cts = Ct(
turb_avg_vels,
Expand Down Expand Up @@ -523,7 +523,7 @@ def cc_solver(
different_average = average_velocity(
flow_field.u_sorted,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)
axial_induction_i = axial_induction(
average_velocities=different_average,
Expand Down Expand Up @@ -739,7 +739,7 @@ def full_flow_cc_solver(
turb_avg_vels = average_velocity(
turbine_grid_flow_field.u_sorted,
method=turbine_grid.average_method,
cubature_coefficients=turbine_grid.cubature_coefficients
cubature_weights=turbine_grid.cubature_weights
)
turb_Cts = Ct(
velocities=turb_avg_vels,
Expand Down Expand Up @@ -889,7 +889,7 @@ def turbopark_solver(
average_velocities = average_velocity(
flow_field.u_sorted,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)

Cts = Ct(
Expand Down Expand Up @@ -972,7 +972,7 @@ def turbopark_solver(
average_velocities = average_velocity(
flow_field.u_sorted,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)

yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None]
Expand Down Expand Up @@ -1211,7 +1211,7 @@ def empirical_gauss_solver(
average_velocities = average_velocity(
flow_field.u_sorted,
method=grid.average_method,
cubature_coefficients=grid.cubature_coefficients
cubature_weights=grid.cubature_weights
)

ct_i = Ct(
Expand Down
18 changes: 9 additions & 9 deletions floris/simulation/turbine.py
Original file line number Diff line number Diff line change
Expand Up @@ -430,22 +430,22 @@ def simple_mean(array, axis=0):
def cubic_mean(array, axis=0):
return np.cbrt(np.mean(array ** 3.0, axis=axis))

def simple_cubature(array, cubature_coefficients, axis=0):
weights = cubature_coefficients["CUBcoeff"].flatten()
def simple_cubature(array, cubature_weights, axis=0):
weights = cubature_weights.flatten()
weights = weights * len(weights) / np.sum(weights)
product = (array * weights[None, None, None, :, None])
return simple_mean(product, axis)

def cubic_cubature(array, cubature_coefficients, axis=0):
weights = cubature_coefficients["CUBcoeff"].flatten()
def cubic_cubature(array, cubature_weights, axis=0):
weights = cubature_weights.flatten()
weights = weights * len(weights) / np.sum(weights)
return np.cbrt(np.mean((array**3.0 * weights[None, None, None, :, None]), axis=axis))

def average_velocity(
velocities: NDArrayFloat,
method: str,
ix_filter: NDArrayFilter | Iterable[int] | None = None,
cubature_coefficients: dict | None = None
cubature_weights: dict | None = None
) -> NDArrayFloat:
"""This property calculates and returns the cube root of the
mean cubed velocity in the turbine's rotor swept area (m/s).
Expand All @@ -471,16 +471,16 @@ def average_velocity(

axis = tuple([3 + i for i in range(velocities.ndim - 3)])
if method == "simple-mean":
return simple_mean(velocities, axis=axis)
return simple_mean(velocities, axis)

elif method == "cubic-mean":
return cubic_mean(velocities, axis=axis)
return cubic_mean(velocities, axis)

elif method == "simple-cubature":
return simple_cubature(velocities, cubature_coefficients=cubature_coefficients, axis=axis)
return simple_cubature(velocities, cubature_weights, axis)

elif method == "cubic-cubature":
return cubic_cubature(velocities, cubature_coefficients=cubature_coefficients, axis=axis)
return cubic_cubature(velocities, cubature_weights, axis)

else:
raise ValueError("Incorrect method given.")
Expand Down
2 changes: 1 addition & 1 deletion floris/tools/floris_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -655,7 +655,7 @@ def turbine_average_velocities(self) -> NDArrayFloat:
return average_velocity(
velocities=self.floris.flow_field.u,
method=self.floris.grid.average_method,
cubature_coefficients=self.floris.grid.cubature_coefficients
cubature_weights=self.floris.grid.cubature_weights
)

@property
Expand Down
3 changes: 1 addition & 2 deletions tests/reg_tests/gauss_regression_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -274,8 +274,7 @@ def test_regression_tandem(sample_inputs_fixture):

farm_avg_velocities = average_velocity(
floris.flow_field.u,
method="cubic-mean",
cubature_coefficients=floris.grid.cubature_coefficients
method="cubic-mean"
)
farm_eff_velocities = rotor_effective_velocity(
floris.flow_field.air_density,
Expand Down

0 comments on commit 92ab0df

Please sign in to comment.