Implementation of a fast simulation core for kite power systems using numba and 2d arrays

FastKPS_01.py

# 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
V_wind           =  0 # (westwind, downwind direction to the east)
V_wind_gnd       =  1
V_wind_tether    =  2 # wind at half of the height of the kite
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):
        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(pro._winch.d_tether)

        self.vec3 = np.zeros((22, 3))
        self.vec3[V_wind, 0]        = V_WIND # (westwind, downwind direction to the east)
        self.vec3[V_wind_gnd, 0]    = V_WIND # (westwind, downwind direction to the east)
        self.vec3[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, v_segment, unit_vector, rho, last_tether_drag, v_app_perp, area):
    """ Calculate the tether drag of a segment. """
    sub3(vec3[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, 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(vec3[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, 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, 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, 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, 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, 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)