Last active
August 29, 2015 14:18
-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -*- 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