-
-
Save sklam/11049399 to your computer and use it in GitHub Desktop.
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 -*- | |
| """ | |
| Test and wrapper script for simulator core functions. | |
| To recompile the C++ library type ! ./make.sh from the ipython console and restart | |
| the ipython kernel. | |
| """ | |
| import os | |
| import sys | |
| path = os.path.dirname(os.path.realpath(__file__)) | |
| sys.path.append(path + '/..') | |
| from scipy.interpolate import InterpolatedUnivariateSpline | |
| import numpy.linalg as la | |
| import numpy as np | |
| import math | |
| from math import pi, acos | |
| from Timer import Timer | |
| from cffi import FFI | |
| from numba import int_, double, void, jit | |
| from Settings import SEGMENTS, PARTICLES, SPRINGS_F, D_TETHER, ELEVATION, KITE_PARTICLES, L_0, \ | |
| REL_SIDE_AREA, ALPHA_ZTIP | |
| ffi = FFI() | |
| # Simulator core functions as defined in sim_core_.cpp. | |
| ffi.cdef("double calc_cl(double alpha);") | |
| ffi.cdef("double calc_cd(double alpha);") | |
| ffi.cdef("void init_cl(double *alpha_cl, double *cl_list, int size);") | |
| ffi.cdef("void init_cd(double *alpha_cl, double *cl_list, int size);") | |
| ffi.cdef("void calcDrag(double *v_wt, double *vp, double *uv, double l_segment, double rho, \ | |
| double d_tether, double *drag);") | |
| ffi.cdef("void calcParticleForces(double *pos1, double *pos2, double *vel1, double *vel2, \ | |
| double *v_wind_tether, double* spring, double* forces, \ | |
| double stiffnes_factor, int segments, double d_tether, double rho, int i);") | |
| ffi.cdef("void calcAeroForces4(double *pos, double *vel, double *v_wind, double rho, double alpha_depower, \ | |
| double rel_steering, double alpha_zero, int segments, double KS, \ | |
| double ALPHA_ZTIP, double AREA, double REL_SIDE_AREA, double* forces, double* xyz, double* \ | |
| last_alphas, double *lift, double *drag);") | |
| ffi.cdef("void innerLoop(double *pos, double *vel, double *v_wind, double *v_wind_tether, double* SPRINGS, double* forces, \ | |
| double stiffnes_factor, int segments, double d_tether);") | |
| # linalg helpers | |
| ffi.cdef("double dot3(double a0, double a1, double a2, double b0, double b1, double b2);") | |
| ffi.cdef("double norm3(double in1, double in2, double in3);") | |
| ffi.cdef("void normalize3(double in1, double in2, double in3, double *out);") | |
| ffi.cdef("double calcWindHeight(double v_wind_gnd, double height);") | |
| ffi.cdef("double calcWindFactor(double height);") | |
| # load the C++ library | |
| lib = ffi.dlopen(path + "/sim_core_.so") | |
| # create a pointer to an output buffer array | |
| buf0 = np.zeros(3) | |
| buf = ffi.cast("double *", buf0.ctypes.data) | |
| res1 = np.zeros(3) | |
| K_D = 1.0 | |
| """ 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*K_D, 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*K_D, 1.0, 0.5, 0.5, 0.5] | |
| calc_cl_ = InterpolatedUnivariateSpline(ALPHA_CL, CL_LIST) | |
| calc_cd_ = InterpolatedUnivariateSpline(ALPHA_CD, CD_LIST) | |
| # random values for alpha in the range of -5.0 to 25.0 degrees | |
| RANDOM = np.random.random_sample((10000,)) * 30.0 - 5.0 | |
| SPRINGS_ = ffi.cast("double *", SPRINGS_F.ctypes.data) | |
| AREA = 16.5 # projected kite area [m^2] | |
| KD = 40.0 / 180.0 * math.pi # max depower | |
| KS = 20.0 / 180.0 * math.pi # max steering | |
| K = 1 - REL_SIDE_AREA # correction factor for the drag | |
| N = 1024 | |
| alpha_cl = np.linspace(-180.0, 180.0, N, endpoint=True) | |
| cl_list = np.zeros(N) | |
| alpha_cd = np.linspace(-180.0, 180.0, N, endpoint=True) | |
| cd_list = np.zeros(N) | |
| # pylint: disable=E1101 | |
| d_array = double[:] | |
| dd_array = double[:, :] | |
| def init(): | |
| """ Calculate the initial conditions y0, yd0 and sw0. Tether with the given elevation angle, | |
| particel zero fixed at origin. """ | |
| DELTA = 1e-6 | |
| # self.setCL_CD(10.0 / 180.0 * math.pi) | |
| 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])) | |
| pos_array, vel_array = np.array(pos, order='F'), np.array(vel, order='F') | |
| # kite particles (pos and vel) | |
| for i in range (KITE_PARTICLES): | |
| pos_array = np.vstack((pos_array, (PARTICLES[i+2]))) # Initial state vector | |
| # kite particles (vel and acc) | |
| for i in range (KITE_PARTICLES): | |
| vel_array = np.vstack((vel_array, (np.array([DELTA, DELTA, DELTA])))) # Initial state vector | |
| # return state_y0 | |
| return( pos_array, vel_array) | |
| def init_cl(): | |
| for i in range(N): | |
| cl_list[i] = calc_cl_(alpha_cl[i]) | |
| lib.init_cl(ffi.cast("double *", alpha_cl.ctypes.data), ffi.cast("double *", cl_list.ctypes.data), N) | |
| def init_cd(): | |
| for i in range(N): | |
| cd_list[i] = calc_cd_(alpha_cd[i]) | |
| lib.init_cd(ffi.cast("double *", alpha_cd.ctypes.data), ffi.cast("double *", cd_list.ctypes.data), N) | |
| def dot_(vec1, vec2): | |
| """ Calculate the dot product of two 3d vectors. | |
| Not very fast. Only for testing the C function. """ | |
| return lib.dot3(vec1[0], vec1[1], vec1[2], vec2[0], vec2[1], vec2[2]) | |
| dot = jit(double(double[:], double[:]))(dot_) | |
| def norm_(vec3): | |
| return lib.norm3(vec3[0], vec3[1], vec3[2]) | |
| norm = jit(double(double[:]))(norm_) | |
| def normalize_(vec3): | |
| res = np.zeros(3) | |
| lib.normalize3(vec3[0], vec3[1], vec3[2], buf) | |
| res[:] = buf0[:] | |
| return res | |
| normalize = jit(double[:](double[:]))(normalize_) | |
| def cross_(vec1, vec2): | |
| """ Calculate the dot product of two 3d vectors. """ | |
| result = np.zeros(3) | |
| a1, a2, a3 = double(vec1[0]), double(vec1[1]), double(vec1[2]) | |
| b1, b2, b3 = double(vec2[0]), double(vec2[1]), double(vec2[2]) | |
| result[0] = a2 * b3 - a3 * b2 | |
| result[1] = a3 * b1 - a1 * b3 | |
| result[2] = a1 * b2 - a2 * b1 | |
| return result | |
| cross = jit(double[:](double[:], double[:]))(cross_) | |
| def calcDrag_(v_wind_tether, v_particle, unit_vector, l_segment, rho): | |
| v_apparent = v_wind_tether - v_particle | |
| v_app_norm = norm(v_apparent) | |
| area = l_segment * D_TETHER * (1 - 0.9 * abs(dot(unit_vector, v_apparent) / v_app_norm)) | |
| return -0.5 * rho * v_app_norm * area * v_apparent | |
| #extern "C" void calcDrag(double *v_wt, double *vp, double *uv, double l_segment, double rho, \ | |
| # double d_tether, double *drag) | |
| def calcDrag2(v_wind_tether, v_particle, unit_vector, l_segment, rho, res): | |
| d_tether = D_TETHER | |
| v_wt = ffi.cast("double *", v_wind_tether.ctypes.data) | |
| vp = ffi.cast("double *", v_particle.ctypes.data) | |
| uv = ffi.cast("double *", unit_vector.ctypes.data) | |
| lib.calcDrag(v_wt, vp, uv, l_segment, rho, d_tether, buf) | |
| res[:] = buf0[:] | |
| return res | |
| def calcParticleForces_(pos1, pos2, vel1, vel2, v_wind_tether, spring, forces, stiffnes_factor, segments, \ | |
| d_tether, i): | |
| """ Calculate the drag force of the tether segment, defined by the parameters po1, po2, vel1 and vel2 | |
| and distribute it equally on the two particles, that are attached to the segment. | |
| The result is stored in the array self.forces. """ | |
| p_1 = int_(spring[0]) # Index of point nr. 1 | |
| p_2 = int_(spring[1]) # Index of point nr. 2 | |
| l_0 = spring[2] # Unstressed length | |
| k = spring[3] * stiffnes_factor # Spring constant | |
| c = spring[4] # Damping coefficient | |
| #print 'k, c: ', k, c | |
| segment = pos1 - pos2 | |
| rel_vel = vel1 - vel2 | |
| av_vel = 0.5 * (vel1 + vel2) | |
| norm1 = norm(segment) | |
| unit_vector = segment / norm1 | |
| k1 = 0.25 * k # compression stiffness kite segments | |
| k2 = 0.1 * k # compression stiffness tether segments | |
| # look at: http://en.wikipedia.org/wiki/Vector_projection | |
| # calculate the relative velocity in the direction of the spring (=segment) | |
| spring_vel = dot(unit_vector, rel_vel) | |
| #print 'norm1, l_0: ', norm1, l_0 | |
| #print "spring_vel: ", spring_vel | |
| #print "unit_vector: ", unit_vector | |
| if (norm1 - l_0) > 0.0: | |
| spring_force = (k * (norm1 - l_0) + (c * spring_vel)) * unit_vector | |
| elif i >= segments: # kite spring | |
| spring_force = (k1 * (norm1 - l_0) + (c * spring_vel)) * unit_vector | |
| else: | |
| spring_force = (k2 * (norm1 - l_0) + (c * spring_vel)) * unit_vector | |
| #print "spring_force: ", spring_force | |
| # Aerodynamic damping for particles of the tether and kite | |
| v_apparent = v_wind_tether - av_vel | |
| v_app_norm = norm(v_apparent) | |
| area = norm1 * d_tether * (1 - 0.9 * abs(dot(unit_vector, v_apparent) / v_app_norm)) | |
| half_drag_force = -0.25 * 1.25 * v_app_norm * area * v_apparent | |
| forces[p_1] += half_drag_force + spring_force | |
| forces[p_2] += half_drag_force - spring_force | |
| calcParticleForces2 = jit(void(double[:], double[:], double[:], double[:], double[:], double[:], \ | |
| double[:,:], double, int_, double, int_)) (calcParticleForces_) | |
| def calcParticleForces(pos1, pos2, vel1, vel2, v_wind_tether, spring, forces, stiffnes_factor, segments, \ | |
| d_tether, i): | |
| pos1_ = ffi.cast("double *", pos1.ctypes.data) | |
| pos2_ = ffi.cast("double *", pos2.ctypes.data) | |
| vel1_ = ffi.cast("double *", vel1.ctypes.data) | |
| vel2_ = ffi.cast("double *", vel2.ctypes.data) | |
| v_wind_tether_ = ffi.cast("double *", v_wind_tether.ctypes.data) | |
| spring_ = ffi.cast("double *", spring.ctypes.data) | |
| forces_ = ffi.cast("double *", forces.ctypes.data) | |
| lib.calcParticleForces(pos1_, pos2_, vel1_, vel2_, \ | |
| v_wind_tether_, spring_, forces_, \ | |
| stiffnes_factor, segments, d_tether, i) | |
| def innerLoop2_(pos, vel, v_wind, v_wind_tether, forces, stiffnes_factor, segments, d_tether): | |
| for i in xrange(SPRINGS_F.shape[0]): | |
| p_1 = int_(SPRINGS_F[i, 0]) # First point nr. | |
| p_2 = int_(SPRINGS_F[i, 1]) # Second point nr. | |
| height = 0.5 * (pos[p_1][2] + pos[p_2][2]) | |
| v_wind_tether[0] = lib.calcWindHeight(v_wind[0], height) | |
| v_wind_tether[1] = lib.calcWindHeight(v_wind[1], height) | |
| v_wind_tether[2] = lib.calcWindHeight(v_wind[2], height) | |
| # print "height, v_wind_tether: ", height, v_wind_tether | |
| calcParticleForces_(pos[p_1], pos[p_2], vel[p_1], vel[p_2], v_wind_tether, SPRINGS_F[i], forces, \ | |
| stiffnes_factor, segments, d_tether, i) | |
| innerLoop2 = jit(void(double[:, :], double[:, :], double[:], double[:], double[:, :], double, int_, double))(innerLoop2_) | |
| @jit | |
| class FastPSS(object): | |
| """ Fast particle system solver. """ | |
| @void() | |
| def __init__(self): | |
| self.pos = np.zeros((SEGMENTS + KITE_PARTICLES + 1,3), order = 'F') | |
| self.vel = np.zeros((SEGMENTS + KITE_PARTICLES + 1,3), order = 'F') | |
| self.v_wind_tether = np.zeros(3) | |
| self.v_wind = np.zeros(3) | |
| self.last_alphas = d_array(np.array((0.1, ALPHA_ZTIP, ALPHA_ZTIP))) | |
| self.forces = dd_array(np.zeros((SEGMENTS + PARTICLES.shape[0] - 1, 3), order='F')) | |
| self.xyz = dd_array(np.zeros((3,3), order='F')) | |
| self.pos_ = ffi.cast("double *", self.pos.ctypes.data) | |
| self.vel_ = ffi.cast("double *", self.vel.ctypes.data) | |
| self.v_wind_tether_ = ffi.cast("double *", self.v_wind_tether.ctypes.data) | |
| self.v_wind_ = ffi.cast("double *", self.v_wind.ctypes.data) | |
| self.forces_ = ffi.cast("double *", self.forces.ctypes.data) | |
| self.xyz_ = ffi.cast("double *", self.xyz.ctypes.data) | |
| self.last_alphas_ = ffi.cast("double *", self.last_alphas.ctypes.data) | |
| self.lift_ = ffi.new("double *") | |
| self.lift_[0] = 0.0 | |
| self.drag_ = ffi.new("double *") | |
| self.drag_[0] = 0.0 | |
| @void(double[:, :], double[:, :], double[:], double[:], double, int_, double) | |
| def innerLoop(self, pos, vel, v_wind, v_wind_tether, stiffnes_factor, segments, d_tether): | |
| self.pos[:] = pos | |
| self.vel[:] = vel | |
| self.v_wind[:] = v_wind | |
| self.v_wind_tether[:] = v_wind_tether | |
| lib.innerLoop(self.pos_, self.vel_, self.v_wind_, self.v_wind_tether_, SPRINGS_, self.forces_, stiffnes_factor, segments, d_tether) | |
| @void( double[:,:], double[:,:], double[:], double, double, double, double, int_) | |
| def calcAeroForces4_(self, pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS): | |
| """ | |
| pos0, pos3, pos4: position of the kite particles P0, P3, and P4 | |
| v2, v3, v4: velocity of the kite particles P2, P3, and P4 | |
| rho: air density [kg/m^3] | |
| rel_depower: value between 0.0 and 1.0 | |
| rel_steering: value between -1.0 and +1.0 | |
| """ | |
| self.pos[:] = pos | |
| self.vel[:] = vel | |
| pos2, pos3, pos4 = pos[SEGMENTS+2], pos[SEGMENTS+3], pos[SEGMENTS+4] | |
| v2, v3, v4 = vel[SEGMENTS+2], vel[SEGMENTS+3], vel[SEGMENTS+4] | |
| va_2, va_3, va_4 = v_wind - v2, v_wind - v3, v_wind - v4 | |
| pos_centre = 0.5 * (pos3 + pos4) | |
| delta = pos2 - pos_centre | |
| z = -normalize(delta) | |
| y = normalize(pos3 - pos4) | |
| x = cross(y, z) | |
| va_xz2 = va_2 - dot(va_2, y) * y | |
| va_xy3 = va_3 - dot(va_3, z) * z | |
| va_xy4 = va_4 - dot(va_4, z) * z | |
| alpha_2 = (pi - acos(dot(normalize(va_xz2), x)) - alpha_depower ) * 360.0 / pi + alpha_zero | |
| alpha_3 = (pi - acos(dot(normalize(va_xy3), x)) - rel_steering * KS) * 360.0 / pi + ALPHA_ZTIP | |
| alpha_4 = (pi - acos(dot(normalize(va_xy4), x)) + rel_steering * KS) * 360.0 / pi + ALPHA_ZTIP | |
| CL2, CD2 = double(lib.calc_cl(alpha_2)), double(lib.calc_cd(alpha_2)) | |
| CL3, CD3 = double(lib.calc_cl(alpha_3)), double(lib.calc_cd(alpha_3)) | |
| CL4, CD4 = double(lib.calc_cl(alpha_4)), double(lib.calc_cd(alpha_4)) | |
| L2 = -0.5 * rho * (norm(va_xz2))**2 * AREA * CL2 * normalize(cross(va_2, y)) | |
| L3 = -0.5 * rho * (norm(va_xy3))**2 * AREA * REL_SIDE_AREA * CL3 * normalize(cross(va_3, z)) | |
| L4 = -0.5 * rho * (norm(va_xy4))**2 * AREA * REL_SIDE_AREA * CL4 * normalize(cross(z, va_4)) | |
| D2 = -0.5 * K * rho * norm(va_2) * AREA * CD2 * va_2 | |
| D3 = -0.5 * K * rho * norm(va_3) * AREA * REL_SIDE_AREA * CD3 * va_3 | |
| D4 = -0.5 * K * rho * norm(va_4) * AREA * REL_SIDE_AREA * CD4 * va_4 | |
| self.forces[SEGMENTS + 2] += (L2 + D2) | |
| self.forces[SEGMENTS + 3] += (L3 + D3) | |
| self.forces[SEGMENTS + 4] += (L4 + D4) | |
| @void( double[:,:], double[:,:], double[:], double, double, double, double, int_, \ | |
| double, double, double, double) | |
| def calcAeroForces4(self, pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS, \ | |
| KS, ALPHA_ZTIP, AREA, REL_SIDE_AREA): | |
| """ | |
| pos0, pos3, pos4: position of the kite particles P0, P3, and P4 | |
| v2, v3, v4: velocity of the kite particles P2, P3, and P4 | |
| rho: air density [kg/m^3] | |
| rel_depower: value between 0.0 and 1.0 | |
| rel_steering: value between -1.0 and +1.0 | |
| """ | |
| self.pos[:] = pos | |
| self.vel[:] = vel | |
| self.v_wind[:] = v_wind | |
| self.xyz.fill(0.0) | |
| lib.calcAeroForces4(self.pos_, self.vel_, self.v_wind_, rho, alpha_depower, \ | |
| rel_steering, alpha_zero, SEGMENTS, KS, \ | |
| ALPHA_ZTIP, AREA, REL_SIDE_AREA, self.forces_, self.xyz_, self.last_alphas_, self.lift_, self.drag_) | |
| PSS = FastPSS() | |
| def test_calcAeroForces4(): | |
| pos, vel = init() | |
| rho = 1.25 | |
| v_wind = np.array((8.0, 0.2, 0.0)) | |
| alpha_depower = 0.1 | |
| rel_steering = -0.1 | |
| alpha_zero = 5.0 | |
| PSS.forces.fill(0.0) | |
| PSS.calcAeroForces4_(pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS) | |
| # print PSS.forces | |
| forces1=np.array(PSS.forces) | |
| PSS.forces.fill(0.0) | |
| PSS.calcAeroForces4(pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS,\ | |
| KS, ALPHA_ZTIP, AREA, REL_SIDE_AREA) | |
| # print PSS.forces | |
| assert la.norm(PSS.forces - forces1) <= 1e-10 | |
| print "Fished test_calcAeroForces4()" | |
| def innerLoop_(pos, vel, v_wind_tether, forces, stiffnes_factor, segments, d_tether): | |
| pos_ = ffi.cast("double *", pos.ctypes.data) | |
| vel_ = ffi.cast("double *", vel.ctypes.data) | |
| v_wind_tether_ = ffi.cast("double *", v_wind_tether.ctypes.data) | |
| forces_ = ffi.cast("double *", forces.ctypes.data) | |
| lib.innerLoop(pos_, vel_, v_wind_tether_, SPRINGS_, forces_, stiffnes_factor, segments, d_tether) | |
| innerLoop = jit(void(double[:, :], double[:, :], double[:], double[:, :], double, int_, double))(innerLoop_) | |
| def test_calcParticleForces(): | |
| pos1 = np.array((1.0, 2.0, 3.0)) | |
| pos2 = np.array((2.0, 3.0, 4.0)) | |
| vel1 = np.array((3.0, 4.0, 5.0)) | |
| vel2 = np.array((4.0, 5.0, 6.0)) | |
| for i in range(SPRINGS_F.shape[0]): | |
| spring = SPRINGS_F[i] | |
| forces = dd_array(np.zeros((SEGMENTS + PARTICLES.shape[0] - 1, 3), order='F')) | |
| stiffnes_factor = 0.5 | |
| v_wind_tether = np.array((8.0, 0.1, 0.0)) | |
| calcParticleForces_(pos1, pos2, vel1, vel2, v_wind_tether, spring, forces, stiffnes_factor, SEGMENTS, \ | |
| D_TETHER, 1) | |
| forces1=np.array(forces) | |
| #print forces1 | |
| forces.fill(0.0) | |
| calcParticleForces(pos1, pos2, vel1, vel2, v_wind_tether, spring, forces, stiffnes_factor, SEGMENTS, \ | |
| D_TETHER, 1) | |
| assert la.norm(forces - forces1) <= 1e-10 | |
| print "Fished test_calcParticleForces()" | |
| def test_innerLoop(): | |
| pos, vel = init() | |
| # print 'pos.shape: ', pos.shape | |
| # print SEGMENTS + KITE_PARTICLES + 1 | |
| # print "pos.shape", pos.shape | |
| v_wind = np.array((7.0, 0.1, 0.0)) | |
| v_wind_tether = np.array((8.0, 0.1, 0.0)) | |
| stiffnes_factor = 0.5 | |
| forces = dd_array(np.zeros((SEGMENTS + PARTICLES.shape[0] - 1, 3), order='F')) | |
| innerLoop2_(pos, vel, v_wind, v_wind_tether, forces, stiffnes_factor, SEGMENTS, D_TETHER) | |
| forces1=np.array(forces) | |
| #print forces1 | |
| PSS.forces.fill(0.0) | |
| PSS.innerLoop(np.ascontiguousarray(pos), vel, v_wind, v_wind_tether, stiffnes_factor, SEGMENTS, D_TETHER) | |
| #print forces | |
| assert la.norm(forces - forces1) <= 1e-10 | |
| print "Fished test_innerLoop()" | |
| init_cl() | |
| init_cd() | |
| # test_calcParticleForces() | |
| test_innerLoop() | |
| test_calcAeroForces4() | |
| if __name__ == '__main__': | |
| vec1 = np.array((1.0, 2.0, 3.0)) | |
| vec2 = np.array((2.0, 3.0, 4.0)) | |
| vec3 = np.array((3.0, 4.0, 5.0)) | |
| vec4 = np.array((4.0, 5.0, 6.0)) | |
| spring = SPRINGS_F[0] | |
| stiffnes_factor = 0.5 | |
| v_wind_tether = vec1 | |
| v_wind = np.array((8.0, 0.2, 0.0)) | |
| alpha_depower = 0.1 | |
| rel_steering = -0.1 | |
| alpha_zero = 5.0 | |
| v_particle = np.array((-1.0, 2.5, 3.0)) | |
| l_segment = 100.0 | |
| rho = 1.25 | |
| v_wind_tether = vec1 | |
| v_particle = np.array((-1.0, 2.5, 3.0)) | |
| unit_vector = normalize(np.array((0.5, 0.1, 9.0))) | |
| l_segment = 100.0 | |
| rho = 1.25 | |
| drag1 = calcDrag_(v_wind_tether, v_particle, unit_vector, l_segment, rho) | |
| drag2 = calcDrag2(v_wind_tether, v_particle, unit_vector, l_segment, rho, res1) | |
| # print drag1 | |
| # print drag2 | |
| pos, vel = init() | |
| # print "pos.shape", pos.shape | |
| v_wind_tether = np.array((8.0, 0.1, 0.0)) | |
| forces = dd_array(np.zeros((SEGMENTS + PARTICLES.shape[0] - 1, 3), order='F')) | |
| stiffnes_factor = 0.5 | |
| print "calcWindHeight(8.0, 500.0): ", lib.calcWindHeight(8.0, 500.0) | |
| if True: | |
| with Timer() as t0: | |
| for i in xrange(10000): | |
| PSS.forces.fill(0.0) | |
| print "time for empty loop [µs] ", t0.secs / 10000 * 1e6 | |
| with Timer() as t1: | |
| for i in xrange(10000): | |
| PSS.forces.fill(0.0) | |
| PSS.calcAeroForces4_(pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS) | |
| print "calcAeroForces4: [µs] ", (t1.secs - t0.secs) / 10000 * 1e6 | |
| with Timer() as t1: | |
| for i in xrange(10000): | |
| PSS.forces.fill(0.0) | |
| PSS.calcAeroForces4(pos, vel, v_wind, rho, alpha_depower, rel_steering, alpha_zero, SEGMENTS, \ | |
| KS, ALPHA_ZTIP, AREA, REL_SIDE_AREA) | |
| print "calcAeroForces4 cffi: [µs] ", (t1.secs - t0.secs) / 10000 * 1e6 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment