Commit 6a6aad1e authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

parallelize aep4smart_start

parent c1c27749
Pipeline #32796 passed with stages
in 20 minutes and 15 seconds
......@@ -297,23 +297,17 @@ class ConvexBoundaryComp(BoundaryBaseComp):
unit_normals = self.unit_normals
# initialize array to hold distances from each point to each face
face_distance = np.zeros([nPoints, nVertices])
from numpy import newaxis as na
# loop through points and find distances to each face
for i in range(0, nPoints):
# define the vector from the point of interest to the first point of the face
PA = (vertices[:, na] - points[na])
# determine if point is inside or outside of each face, and distances from each face
for j in range(0, nVertices):
# define the vector from the point of interest to the first point of the face
pa = np.array([vertices[j, 0] - points[i, 0], vertices[j, 1] - points[i, 1]])
# find perpendicular distances from point to current surface (vector projection)
d_vec = np.vdot(pa, unit_normals[j]) * unit_normals[j]
# calculate the sign of perpendicular distances from point to current face (+ is inside, - is outside)
face_distance[i, j] = np.vdot(d_vec, unit_normals[j])
# print (face_distance)
return face_distance
# find perpendicular distances from point to current surface (vector projection)
dist = np.sum(PA * unit_normals[:, na], 2)
# calculate the sign of perpendicular distances from point to current face (+ is inside, - is outside)
d_vec = dist[:, :, na] * unit_normals[:, na]
face_distance = np.sum(d_vec * unit_normals[:, na], 2)
return face_distance.T
def distances(self, x, y):
return self.calculate_distance_to_boundary(np.array([x, y]).T)
......@@ -384,7 +378,8 @@ class PolygonBoundaryComp(BoundaryBaseComp):
dEdgeDist_dy = dx / length
edge_vect_j = np.array([x2 - x1, y2 - y1])
edge_vect_len_j = np.linalg.norm(edge_vect_j, axis=0)
return (xy_boundary, x1, y1, x2, y2, dEdgeDist_dx, dEdgeDist_dy, edge_unit_vec, dx, xy1_vec, xy2_vec, edge_vect_j, edge_vect_len_j)
return (xy_boundary, x1, y1, x2, y2, dEdgeDist_dx, dEdgeDist_dy,
edge_unit_vec, dx, xy1_vec, xy2_vec, edge_vect_j, edge_vect_len_j)
def _calc_distance_and_gradients(self, x, y, boundary_properties=None):
"""
......@@ -640,7 +635,8 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
turns_left_ij = np.broadcast_to(np.cross(np.roll(edge_vect_j, 1, axis=1), edge_vect_j, axis=0) > 0, shape_ij)
turns_right_ij = np.logical_not(turns_left_ij)
De_ij = np.abs((x2[na, :] - x1[na, :]) * (y1[na, :] - y[:, na]) - (x1[na, :] - x[:, na]) * (y2[na, :] - y1[na, :])) / np.sqrt((x2[na, :] - x1[na, :]) ** 2 + (y2[na, :] - y1[na, :]) ** 2)
De_ij = np.abs((x2[na, :] - x1[na, :]) * (y1[na, :] - y[:, na]) - (x1[na, :] - x[:, na]) *
(y2[na, :] - y1[na, :])) / np.sqrt((x2[na, :] - x1[na, :]) ** 2 + (y2[na, :] - y1[na, :]) ** 2)
Dp_ij = np.sqrt((x[:, na] - x1[na, :]) ** 2 + (y[:, na] - y1[na, :]) ** 2)
D_ij[overlapping_ij] = De_ij[overlapping_ij]
......@@ -658,17 +654,20 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
D_ij[before_ij & turns_left_ij] *= -1
D_ij[after_ij & np.roll(turns_left_ij, -1, axis=1)] *= -1
D_ij[before_ij & turns_right_ij & np.roll(outside_edge_ij, 1, axis=1) & outside_edge_ij] *= -1
D_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) & np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
D_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) &
np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
dDdx_ij[before_ij & turns_left_ij] *= -1
dDdx_ij[after_ij & np.roll(turns_left_ij, -1, axis=1)] *= -1
dDdx_ij[before_ij & turns_right_ij & np.roll(outside_edge_ij, 1, axis=1) & outside_edge_ij] *= -1
dDdx_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) & np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
dDdx_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) &
np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
dDdy_ij[before_ij & turns_left_ij] *= -1
dDdy_ij[after_ij & np.roll(turns_left_ij, -1, axis=1)] *= -1
dDdy_ij[before_ij & turns_right_ij & np.roll(outside_edge_ij, 1, axis=1) & outside_edge_ij] *= -1
dDdy_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) & np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
dDdy_ij[after_ij & np.roll(turns_right_ij, -1, axis=1) &
np.roll(outside_edge_ij, -1, axis=1) & outside_edge_ij] *= -1
return D_ij, dDdx_ij, dDdy_ij
......@@ -805,7 +804,11 @@ def main():
plt.figure()
ax = plt.axes(projection='3d')
ax.contour3D(x.reshape(N_points, N_points), y.reshape(N_points, N_points), distances.reshape(N_points, N_points), 50)
ax.contour3D(
x.reshape(
N_points, N_points), y.reshape(
N_points, N_points), distances.reshape(
N_points, N_points), 50)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
......
......@@ -11,7 +11,7 @@ import topfarm
import numpy as np
from scipy.interpolate.interpolate import RegularGridInterpolator
import warnings
from py_wake.flow_map import HorizontalGrid
from py_wake.flow_map import HorizontalGrid, Points, XYGrid
import xarray as xr
from py_wake.utils.gradients import autograd
......@@ -58,6 +58,7 @@ class PyWakeAEP(AEPCalculator):
class PyWakeAEPCostModelComponent(AEPCostModelComponent):
def __init__(self, windFarmModel, n_wt, wd=None, ws=None, max_eval=None, grad_method=autograd, n_cpu=1, **kwargs):
self.windFarmModel = windFarmModel
self.n_cpu = n_cpu
def aep(**kwargs):
return self.windFarmModel.aep(x=kwargs[topfarm.x_key],
......@@ -94,13 +95,9 @@ class PyWakeAEPCostModelComponent(AEPCostModelComponent):
def get_aep4smart_start(self, ws=[6, 8, 10, 12, 14], wd=np.arange(360)):
def aep4smart_start(X, Y, wt_x, wt_y, type=0):
sim_res = self.windFarmModel(wt_x, wt_y, type=type, wd=wd, ws=ws)
x = np.sort(np.unique(X))
y = np.sort(np.unique(Y))
aep_map = sim_res.flow_map(HorizontalGrid(x, y)).aep_xy()
if isinstance(aep_map, xr.DataArray):
aep_map = aep_map.values[:, :, 0]
return RegularGridInterpolator((y, x), aep_map)(np.array([Y, X]).T)
sim_res = self.windFarmModel(wt_x, wt_y, type=type, wd=wd, ws=ws, n_cpu=self.n_cpu)
H = np.full(X.shape, self.windFarmModel.windTurbines.hub_height())
return sim_res.aep_map(Points(X, Y, H), n_cpu=self.n_cpu).values
return aep4smart_start
......
......@@ -72,6 +72,25 @@ def testMultiPolygon():
plot_comp.show()
def test_calculate_distance_to_boundary():
import matplotlib.pyplot as plt
boundary = np.array([(0, 0), (10, 0), (20, 10), (20, 20), (0, 20)])
points = np.array([(2, 10), (10, 21), (14, 6)])
boundary_constr = XYBoundaryConstraint(boundary, 'convex_hull').get_comp(10)
import numpy.testing as npt
if 0:
plt.plot(*boundary.T, )
plt.plot(*points.T, '.')
plt.axis('equal')
plt.grid()
plt.show()
npt.assert_array_almost_equal(boundary_constr.calculate_distance_to_boundary(points),
[[18., 10., 2., 10., 12.72792206],
[10., -1., 10., 21., 14.8492424],
[6., 14., 14., 6., 1.41421356]])
def testDistanceRelaxation():
boundary = [([(0, 0), (5, 0), (5, 2), (3, 2), (3, 1), (2, 1), (2, 2), (0, 2), (0, 0)], 1),
([(3.5, 0.5), (4.5, 0.5), (4.5, 1.5), (3.5, 1.5)], 1),
......
......@@ -4,6 +4,7 @@ from numpy import newaxis as na
import multiprocessing
import threading
import time
from tqdm import tqdm
def smart_start(XX, YY, ZZ, N_WT, min_space, radius=None, random_pct=0, plot=False, seed=None):
......@@ -56,12 +57,15 @@ def smart_start(XX, YY, ZZ, N_WT, min_space, radius=None, random_pct=0, plot=Fal
np.random.seed(seed)
seed = np.random.randint(0, 2**31)
np.random.seed(seed)
for i in range(N_WT):
for i in tqdm(range(N_WT), desc='Smartstart'):
if arr.shape[1] == 0:
raise Exception('No feasible positions for wt %d' % i)
if ZZ_is_func:
z = ZZ(arr[0], arr[1], xs, ys)
if random_pct < 100:
z = ZZ(arr[0], arr[1], xs, ys)
else:
z = np.zeros_like(arr[0])
else:
z = arr[2]
......@@ -137,7 +141,8 @@ def smooth_max_gradient(X, alpha, axis=0):
Matrix of smooth maximum derivatives.
'''
return np.exp(alpha * X) / np.expand_dims(np.exp(alpha * X).sum(axis=axis), axis) * (1 + alpha * (X - np.expand_dims(smooth_max(X, alpha, axis=axis), axis)))
return np.exp(alpha * X) / np.expand_dims(np.exp(alpha * X).sum(axis=axis), axis) * \
(1 + alpha * (X - np.expand_dims(smooth_max(X, alpha, axis=axis), axis)))
def main():
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment