Skip to content

Instantly share code, notes, and snippets.

@ufechner7
Last active August 29, 2015 14:18
Show Gist options
  • Save ufechner7/538269c3c20e84665f60 to your computer and use it in GitHub Desktop.
Save ufechner7/538269c3c20e84665f60 to your computer and use it in GitHub Desktop.
Demo code to show the speed problem with record arrays in numba 0.18.
# -*- coding: utf-8 -*-
"""
Demo code, showing the problem with using record arrays with Numba 0.18.
I replaced just three of the elements of vec3 with a record array with three elements. The execution time
increases from 22us to 45 us.
Model of a kite-power system in implicit form: residual = f(y, yd, sw)
This model implements a 3D mass-spring system with reel-out. It uses five tether segments (the number can be
configured in the file Settings.py). The kite is modelled as additional mass at the end of the tether.
The spring constant and the damping decrease with the segment length. The aerodynamic kite forces are
calculated, depending on reel-out speed, depower and steering settings.
Scientific background: http://arxiv.org/abs/1406.6218
"""
DEBUG = False
import numpy as np
# pylint: disable=E0611
from numba import double, jit, autojit
import linalg_3d_v2 as la
from linalg_3d_v2 import sum2, sum3, sub2, sub3, norm, mul2, mul3, div3, copy2, dot, \
normalize1, normalize2, cross3, neg_sum, calc_alpha
import math
from math import radians, exp
from Settings import C_SPRING, DAMPING, SEGMENTS, L_0, V_REEL_OUT, MASS, KITE_MASS, KCU_MASS, ELEVATION, \
V_WIND, REL_SIDE_AREA, STEERING_COEFFICIENT, PERIOD_TIME, WINCH_MODEL, \
ALPHA_ZERO, BRIDLE_DRAG, AREA, C2_COR, WINCH_D_TETHER
from scipy.interpolate import InterpolatedUnivariateSpline
from WinchModel import calcAcceleration
from KCU_Sim import calcAlphaDepower
from Timer import Timer
from Approximate_CLCD import approximate_cl, approximate_cd
G_EARTH = 9.81 # gravitational acceleration
RHO_0 = 1.225 # kg / m³
C_0 = -0.0032 # steering offset
C_D_TETHER = 0.958 # tether drag coefficient
L_BRIDLE = 33.4 # sum of the lengths of the bridle lines [m]
K_ds = 1.5 # influence of the depower angle on the steering sensitivity
MAX_ALPHA_DEPOWER = 31.0 # was: 44
ALPHA = 1/7
# pylint: disable=E1101
# pylint: disable=C0326
""" Calculate the aerodynamic kite properties. """
ALPHA_CL = [-180.0, -160.0, -90.0, -20.0, -10.0, -5.0, 0.0, 20.0, 40.0, 90.0, 160.0, 180.0]
CL_LIST = [ 0.0, 0.5, 0.0, 0.08, 0.125, 0.15, 0.2, 1.0, 1.0, 0.0, -0.5, 0.0]
ALPHA_CD = [-180.0, -170.0, -140.0, -90.0, -20.0, 0.0, 20.0, 90.0, 140.0, 170.0, 180.0]
CD_LIST = [ 0.5, 0.5, 0.5, 1.0, 0.2, 0.1, 0.2, 1.0, 0.5, 0.5, 0.5]
if False:
calc_cl = InterpolatedUnivariateSpline(ALPHA_CL, CL_LIST)
calc_cd = InterpolatedUnivariateSpline(ALPHA_CD, CD_LIST)
else:
calc_cl = approximate_cl
calc_cd = approximate_cd
@autojit(nopython = True)
def calcRho(height):
return RHO_0 * exp(-height / 8550.0)
def form(number):
return "{:.2f}".format(number)
@jit(nopython=True)
def calcWindFactor(height):
""" Calculate the wind speed at a given height and reference height. """
return (height / 6.0)**ALPHA
# constants for accessing the array of scalars
ParamCD = 0
ParamCL = 1
Length = 2 # length of one tether segment at zero force
C_spring = 3 # spring constant of one tether segment
Damping = 4 # damping constant of one tether segment
Depower = 5
Steering = 6
Alpha_depower = 7
Alpha_zero = 8
Last_alpha = 9 # unused
Psi = 10 # heading angle in radian
Beta = 11 # elevation angle in radian; initial value about 70 degrees
Cor_steering = 12
D_tether = 13
V_app_norm = 14
Area = 15
Last_v_app_norm_tether = 16
# constants for accessing the array of vec3
Segment = 3
Rel_vel = 4
Av_vel = 5
Unit_vector = 7
Force = 8
Last_force = 9
Lift_force = 10
Drag_force = 11
Spring_force = 12
Total_forces = 13
Last_tether_drag = 14
V_apparent = 15
V_app_perp = 16
Kite_y = 17
Steering_force = 18
Acc = 19
Temp = 20
Vec_z = 21
class FastKPS(object):
""" Class with the inner properties of the residual calculations. """
def __init__(self):
x_dt = np.dtype([('v_wind', np.float64),
('v_wind_gnd', np.float64),
('v_wind_tether', np.float64)])
self.buf = np.zeros((3, 3)) # initialize the record array with zeros
self.rec3 = np.recarray(3, dtype=x_dt, buf=self.buf)
self.scalars = np.zeros(17)
self.scalars[ParamCD] = 1.0
self.scalars[ParamCL] = 0.2
self.scalars[Last_alpha] = double (0.1)
self.scalars[Beta] = double(1.220) # elevation angle in radian; initial value about 70 degrees
self.scalars[D_tether] = double(WINCH_D_TETHER)
self.vec3 = np.zeros((22, 3))
self.rec3.v_wind[0] = V_WIND # (westwind, downwind direction to the east)
self.rec3.v_wind_gnd[0] = V_WIND # (westwind, downwind direction to the east)
self.rec3.v_wind_tether[0] = V_WIND # wind at half of the height of the kite
self.initial_masses = (np.ones(SEGMENTS + 1)) # array of the initial masses of the particles
self.initial_masses *= double(MASS)
self.masses = (np.ones(SEGMENTS + 1))
@jit(nopython = True)
def calcDrag(vec3, rec3, v_segment, unit_vector, rho, last_tether_drag, v_app_perp, area):
""" Calculate the tether drag of a segment. """
sub3(rec3.v_wind_tether, v_segment, vec3[V_apparent])
v_app_norm = norm(vec3[V_apparent])
mul3(dot(vec3[V_apparent], unit_vector), unit_vector, v_app_perp)
sub3(vec3[V_apparent], v_app_perp, v_app_perp)
mul3(- 0.5 * C_D_TETHER * rho * norm(v_app_perp) * area, v_app_perp, last_tether_drag)
return v_app_norm
@jit(nopython = True)
def calcAeroForces(scalars, vec3, rec3, pos_kite, v_kite, rho, rel_steering, v_apparent):
"""
pos_kite: position of the kite
rho: air density [kg/m^3]
paramCD: drag coefficient (function of power settings)
paramCL: lift coefficient (function of power settings)
rel_steering: value between -1.0 and +1.0
"""
sub3(rec3.v_wind, v_kite, v_apparent)
scalars[V_app_norm] = norm(vec3[V_apparent])
normalize2(vec3[V_apparent], vec3[Drag_force])
cross3(pos_kite, vec3[Drag_force], vec3[Kite_y])
normalize1(vec3[Kite_y])
K = 0.5 * rho * scalars[V_app_norm]**2 * AREA
cross3(vec3[Drag_force], vec3[Kite_y], vec3[Temp])
normalize1(vec3[Temp])
mul3(K * scalars[ParamCL], vec3[Temp], vec3[Lift_force])
# some additional drag is created while steering
mul2( K * scalars[ParamCD] * BRIDLE_DRAG * (1.0 + 0.6 * abs(rel_steering)), vec3[Drag_force])
scalars[Cor_steering] = C2_COR / scalars[V_app_norm] * math.sin(scalars[Psi]) * math.cos(scalars[Beta])
mul3(- K * REL_SIDE_AREA * STEERING_COEFFICIENT * (rel_steering + scalars[Cor_steering]), vec3[Kite_y], \
vec3[Steering_force])
neg_sum(vec3[Lift_force], vec3[Drag_force], vec3[Steering_force], vec3[Last_force])
@jit(nopython = True)
def calcRes(scalars, vec3, rec3, pos1, pos2, vel1, vel2, mass, veld, result, i):
""" Calculate the vector res1, that depends on the velocity and the acceleration.
The drag force of each segment is distributed equaly on both particles. """
sub3(pos1, pos2, vec3[Segment])
height = double((pos1[2] + pos2[2]) * 0.5)
rho = double(calcRho(height)) # calculate the air density
sub3(vel1, vel2, vec3[Rel_vel]) # calculate the relative velocity
sum3(vel1, vel2, vec3[Av_vel]) # dito (continuation)
mul2(0.5, vec3[Av_vel]) # calculate the relative velocity
norm1 = norm(vec3[Segment]) # length of the tether segment
div3(norm1, vec3[Segment], vec3[Unit_vector]) # unit vector in the direction of the tether
# look at: http://en.wikipedia.org/wiki/Vector_projection
# calculate the relative velocity in the direction of the spring (=segment)
spring_vel = dot(vec3[Unit_vector], vec3[Rel_vel])
k2 = 0.05 * scalars[C_spring] # compression stiffness tether segments
if norm1 - scalars[Length] > 0.0:
mul3(scalars[C_spring] * (norm1 - scalars[Length]) + scalars[Damping] * spring_vel, \
vec3[Unit_vector], vec3[Spring_force])
else:
mul3(k2 * (norm1 - scalars[Length]) + scalars[Damping] * spring_vel, \
vec3[Unit_vector], vec3[Spring_force])
scalars[Area] = norm1 * scalars[D_tether]
scalars[Last_v_app_norm_tether] = calcDrag(vec3, rec3, vec3[Av_vel], vec3[Unit_vector], \
rho, vec3[Last_tether_drag], vec3[V_app_perp], scalars[Area])
# TODO: create a copy of the drag force !!!
if i == SEGMENTS:
scalars[Area] = L_BRIDLE * scalars[D_tether]
scalars[Last_v_app_norm_tether] = calcDrag(vec3, rec3, vec3[Av_vel], vec3[Unit_vector], \
rho, vec3[Last_tether_drag], vec3[V_app_perp], scalars[Area])
mul3(0.5, vec3[Last_tether_drag], vec3[Temp])
sum2(vec3[Spring_force], vec3[Temp])
sum3(vec3[Temp], vec3[Last_tether_drag], vec3[Force])
else:
mul3(0.5, vec3[Last_tether_drag], vec3[Temp])
sum3(vec3[Temp], vec3[Spring_force], vec3[Force])
sum3(vec3[Force], vec3[Last_force], vec3[Total_forces])
mul3(0.5, vec3[Last_tether_drag], vec3[Temp])
sub2(vec3[Spring_force], vec3[Temp])
copy2(vec3[Temp], vec3[Last_force])
div3(mass, vec3[Total_forces], vec3[Acc]) # create the vector of the spring acceleration
# the derivative of the velocity must be equal to the total acceleration
copy2(vec3[Acc], vec3[Temp])
vec3[Temp][2] -= G_EARTH
sub3(veld, vec3[Temp], result)
@jit(nopython = True)
def loop(scalars, vec3, rec3, initial_masses, masses, pos, vel, posd, veld, res0, res1):
""" Calculate the vector res0 using a vector expression, and calculate res1 using a loop
that iterates over all tether segments. """
mul3(scalars[Length] / L_0, initial_masses, masses)
masses[SEGMENTS] += (KITE_MASS + KCU_MASS)
copy2(pos[0], res0[0]) # res0[0] = pos[0]
copy2(vel[0], res1[0]) # res1[0] = vel[0]
# res0[1:SEGMENTS+1] = vel[1:SEGMENTS+1] - posd[1:SEGMENTS+1]
for i in xrange(1, SEGMENTS+1):
sub3(vel[i], posd[i], res0[i])
for i in xrange(SEGMENTS, 0, -1): # count down from particle = SEGMENTS to 1
calcRes(scalars, vec3, rec3, pos[i], pos[i-1], vel[i], vel[i-1], masses[i], veld[i], res1[i], i)
@jit(nopython = True)
def setCL_CD(scalars, alpha):
""" Calculate the lift and drag coefficient as a function of the relative depower setting. """
angle = alpha * 180.0 / math.pi + ALPHA_ZERO
if angle > 180.0:
angle -= 360.0
if angle < -180.0:
angle += 360.0
scalars[ParamCL] = calc_cl(angle)
scalars[ParamCD] = calc_cd(angle)
@jit(nopython = True)
def calcLoD(scalars, vec3, vec_c, v_app):
"""
Calculate the lift over drag ratio as a function of the direction vector of the last tether
segment, the current depower setting and the apparent wind speed.
"""
normalize2(vec_c, vec3[Vec_z]) # vec_z = la.normalize(vec_c)
alpha = calc_alpha(v_app, vec3[Vec_z]) - scalars[Alpha_depower]
setCL_CD(scalars, alpha)
class KPS(FastKPS):
""" Class, that defines the physics of a kite-power system. """
def __init__(self):
FastKPS.__init__(self)
self.res = np.zeros((SEGMENTS + 1) * 6).reshape((2, -1, 3))
if WINCH_MODEL:
self.res = np.append(self.res, 0.0) # res_length
self.res = np.append(self.res, 0.0) # res_v_reel_out
self.t_0 = 0.0 # relative start time of the current time interval
self.v_reel_out = 0.0
self.last_v_reel_out = 0.0
self.sync_speed = 0.0
self.rec3.v_wind[0] = V_WIND
self.l_tether = L_0 * SEGMENTS
self.pos_kite, self.v_kite = np.zeros(3), np.zeros(3)
self.rho = RHO_0
def clear(self):
""" Clear all member variables. """
self.__init__()
def residual(self, time, state_y, der_yd):
"""
N-point tether model:
Inputs:
State vector state_y = pos0, pos1, ..., posn-1, vel0, vel1, ..., veln-1
Derivative der_yd = vel0, vel1, ..., veln-1, acc0, acc1, ..., accn-1
Output:
Residual res = res0, res1 = pos0, ..., vel0, ...
"""
if WINCH_MODEL:
length = state_y[-2]
v_reel_out = state_y[-1]
lengthd = der_yd[-2]
v_reel_outd = der_yd[-1]
part = state_y[0:-2].reshape((2, -1, 3)) # reshape the state vector to a vector of particles
# reshape the state derivative to a vector of particles derivatives
partd = der_yd[0:-2].reshape((2, -1, 3))
res0, res1 = self.res[0:-2].reshape((2, -1, 3))[0], self.res[0:-2].reshape((2, -1, 3))[1]
else:
part = state_y.reshape((2, -1, 3)) # reshape the state vector to a vector of particles
# reshape the state derivative to a vector of particles derivatives
partd = der_yd.reshape((2, -1, 3))
res0, res1 = self.res[0], self.res[1]
pos, vel = part[0], part[1]
self.pos_kite, self.v_kite = pos[SEGMENTS], vel[SEGMENTS]
posd, veld = partd[0], partd[1]
if WINCH_MODEL:
sync_speed = self.sync_speed
self.scalars[Length] = length / SEGMENTS
else:
delta_t = time - self.t_0
delta_v = self.v_reel_out - self.last_v_reel_out
self.scalars[Length] = (self.l_tether + self.last_v_reel_out * delta_t + 0.5 * delta_v * delta_t**2) \
/ SEGMENTS
self.scalars[C_spring] = C_SPRING * L_0 / self.scalars[Length]
self.scalars[Damping] = DAMPING * L_0 / self.scalars[Length]
calcLoD(self.scalars, self.vec3, pos[SEGMENTS - 1] - self.pos_kite , self.rec3.v_wind - self.v_kite)
calcAeroForces(self.scalars, self.vec3, self.rec3, self.pos_kite, self.v_kite, self.rho, self.scalars[Steering], \
self.vec3[V_apparent]) # force at the kite
loop(self.scalars, self.vec3, self.rec3, self.initial_masses, self.masses, pos, vel, posd, veld, res0, res1)
if WINCH_MODEL:
self.res[-2] = lengthd - v_reel_out
self.res[-1] = v_reel_outd - calcAcceleration(sync_speed, v_reel_out, la.norm(self.vec3[Last_force]))
return self.res
else:
return self.res.flatten()
def setV_ReelOut(self, v_reel_out, t_0, period_time = PERIOD_TIME):
""" Setter for the reel-out speed. Must be called every 50 ms (before each simulation).
It also updates the tether length, therefore it must be called even if v_reelout has
not changed. """
if v_reel_out is None:
v_reel_out = 0.0
self.l_tether += 0.5 * (v_reel_out + self.last_v_reel_out) * period_time
self.last_v_reel_out = self.v_reel_out
self.v_reel_out = v_reel_out
self.t_0 = t_0
def setSyncSpeed(self, sync_speed, t_0):
""" Setter for the reel-out speed. Must be called every 50 ms (before each simulation).
It also updates the tether length, therefore it must be called even if v_reelout has
not changed. """
# self.last_sync_speed = self.sync_speed
self.sync_speed = sync_speed
self.t_0 = t_0
if t_0 <= 5.0:
self.scalars[Alpha_zero] = t_0/5.0 * ALPHA_ZERO
def setDepowerSteering(self, depower, steering):
""" Setter depower and the steering model inputs. Valid range for steering: -1.0 .. 1.0.
Valid range for depower: 0 .. 1.0 """
self.scalars[Steering] = steering
self.scalars[Depower] = depower
self.scalars[Alpha_depower] = calcAlphaDepower(depower) * (MAX_ALPHA_DEPOWER / 31.0)
# print "depower, alpha_depower", form(depower), form(degrees(self.scalars[Alpha_depower]))
# print "v_app_norm, CL, rho: ", form(self.scalars[V_app_norm]),form(self.scalars[ParamCL]), form(self.rho)
self.scalars[Steering] = (steering - C_0) / (1.0 + K_ds * (self.scalars[Alpha_depower] \
/ radians(MAX_ALPHA_DEPOWER)))
# print "LoD: ", self.scalars[ParamCL]/ self.scalars[ParamCD]
def setBetaPsi(self, beta, psi):
self.scalars[Beta] = beta
self.scalars[Psi] = psi
def setL_Tether(self, l_tether):
""" Setter for the tether reel-out lenght (at zero force). """
self.l_tether = l_tether
def getL_Tether(self):
""" Getter for the tether reel-out lenght (at zero force). """
return self.l_tether
def getForce(self):
""" Return the absolute value of the force at the winch as calculated during the last
simulation. """
return la.norm(self.vec3[Last_force]) # last_force is the force at the whinch!
def getSpringForces(self, pos):
""" Return an array of the scalar spring forces of all tether segements.
Input: The vector pos of the positions of the point masses that belong to the tether. """
forces = np.zeros(SEGMENTS)
for i in range(SEGMENTS):
forces[i] = self.scalars[C_spring] * (la.norm(pos[i+1] - pos[i]) - self.scalars[Length])
return forces
def getLiftDrag(self):
return la.norm(self.lift_force), la.norm(self.vec3[Drag_force])
def getV_Wind(self):
""" Return the vector of the wind velocity at the height of the kite. """
return self.rec3.v_wind
def setV_WindGround(self, height, v_wind_gnd=V_WIND, wind_dir=0.0, v_wind = None, rel_turbulence = 0.0):
""" Set the vector of the wind-velocity at the height of the kite. As parameter the height,
the ground wind speed, the wind direction, the turbulence vector and the relative turbulence are needed.
Must be called every 50 ms.
"""
if height < 6.0:
height = 6.0
v_wind0 = v_wind_gnd * minilib.calcWindFactor(height) * math.cos(wind_dir)
v_wind1 = v_wind_gnd * minilib.calcWindFactor(height) * math.sin(wind_dir)
if v_wind is not None:
self.rec3.v_wind[0] = v_wind0 * (1-rel_turbulence) + v_wind[0] * rel_turbulence
self.rec3.v_wind[1] = v_wind1 * (1-rel_turbulence) + v_wind[1] * rel_turbulence
self.rec3.v_wind[2] = v_wind[2] * rel_turbulence
else:
self.rec3.v_wind[0] = v_wind0
self.rec3.v_wind[1] = v_wind1
self.rec3.v_wind[2] = 0.0
# print 'v_wind[0]', self.vec3[V_wind][0], calcWindHeight(v_wind_gnd, height)
self.rec3.v_wind_gnd[0] = v_wind_gnd * math.cos(wind_dir)
self.rec3.v_wind_gnd[1] = v_wind_gnd * math.sin(wind_dir)
# self.vec3[V_wind_tether][0] = calcWindHeight(v_wind_gnd, height / 2.0)
self.rec3.v_wind_tether[0] = v_wind_gnd * minilib.calcWindFactor(height / 2.0) * math.cos(wind_dir)
self.rec3.v_wind_tether[1] = v_wind_gnd * minilib.calcWindFactor(height / 2.0) * math.sin(wind_dir)
self.rho = calcRho(height)
def initTether(self):
print "KPS4.initTether()..................", WINCH_D_TETHER, "m"
self.scalars[D_tether] = WINCH_D_TETHER
def getLoD(self):
lift, drag = self.getLiftDrag()
return lift / drag
def init(self):
""" Calculate the initial conditions y0, yd0 and sw0. Tether with the given elevation angle,
particle zero fixed at origin. """
DELTA = 1e-6
setCL_CD(self.scalars, 10.0/180.0 * math.pi)
print 'paramCL, paramCD: ', self.scalars[ParamCL], self.scalars[ParamCD]
pos, vel, acc = [], [], []
state_y = DELTA
vel_incr = 0
sin_el, cos_el = math.sin(ELEVATION / 180.0 * math.pi), math.cos(ELEVATION / 180.0 * math.pi)
for i in range (SEGMENTS + 1):
radius = - i * L_0
if i == 0:
pos.append(np.array([-cos_el * radius, DELTA, -sin_el * radius]))
vel.append(np.array([DELTA, DELTA, DELTA]))
else:
pos.append(np.array([-cos_el * radius, state_y, -sin_el * radius]))
if i < SEGMENTS:
vel.append(np.array([DELTA, DELTA, -sin_el * vel_incr*i]))
else:
vel.append(np.array([DELTA, DELTA, -sin_el * vel_incr*(i-1.0)]))
acc.append(np.array([DELTA, DELTA, -9.81]))
state_y0, yd0 = np.array([]), np.array([])
for i in range (SEGMENTS + 1):
state_y0 = np.append(state_y0, pos[i]) # Initial state vector
yd0 = np.append(yd0, vel[i]) # Initial state vector derivative
for i in range (SEGMENTS + 1):
state_y0 = np.append(state_y0, vel[i]) # Initial state vector
yd0 = np.append(yd0, acc[i]) # Initial state vector derivative
self.setL_Tether(L_0 * SEGMENTS)
if WINCH_MODEL:
self.setSyncSpeed(V_REEL_OUT, 0.0)
# append the initial reel-out length and it's derivative
state_y0 = np.append(state_y0, L_0 * SEGMENTS)
yd0 = np.append(yd0, V_REEL_OUT)
# append reel-out velocity and it's derivative
state_y0 = np.append(state_y0, V_REEL_OUT)
yd0 = np.append(yd0, DELTA)
else:
self.setV_ReelOut(V_REEL_OUT, 0.0)
self.setV_ReelOut(V_REEL_OUT, 0.0)
return state_y0, yd0
if __name__ == '__main__':
PRINT = True
MY_KPS = KPS()
Y0, YD0 = MY_KPS.init() # basic consistency test; the result should be zero
if PRINT:
print 'Y0:' , Y0
print 'Yd0:', YD0
print Y0.shape, YD0.shape
LOOPS = 100
RES = MY_KPS.residual(0.0, Y0, YD0)
with Timer() as t:
for j in range(LOOPS):
RES = MY_KPS.residual(0.0, Y0, YD0)
if PRINT:
if WINCH_MODEL:
print 'res:', (RES[0:-2].reshape((2, -1, 3)))
print 'res_winch', RES[-2], RES[-1]
else:
print 'res:', (RES.reshape((2, -1, 3)))
print "\nExecution time (µs): ", t.secs / LOOPS * 1e6
# baseline excecution time: 610 us
# was 440 us (without array access)
# now 490 us (with array access)
# now 480 us (removed obj in one function)
# do not use obj in calcRes, use scalar and vec instead: 350 us
# do not use obj in loop, use scalar, vec and intial_masses instead: 310 us
# converted calcRes to nopython mode. Now 122 us.
# replaced obj with scalar as parameter in two functions. Now 79 us
# modified function "loop" to nopython = True. Now 66 us
# replaced calc_cl and calc_cd with approximate_cl and approximate_cd. Now 23 us
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment