Skip to content

Commit

Permalink
I set the fastmath argument to False. This fixes issue #14.
Browse files Browse the repository at this point in the history
  • Loading branch information
@camUrban
camUrban committed Apr 6, 2022
1 parent 5947bb8 commit f17fb1b
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 46 deletions.
23 changes: 3 additions & 20 deletions .idea/watcherTasks.xml
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

42 changes: 21 additions & 21 deletions pterasoftware/aerodynamics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@
# dramatically affect the stability of the result. I'm using this value, as cited for
# use in flapping-wing vehicles in "Role of Filament Strain in the Free-Vortex
# Modeling of Rotor Wakes" (Ananthan and Leishman, 2004). It is unitless.
squire = 10 ** -4
squire = 10**-4

# Set the value of Lamb's constant that will be used by the induced velocity
# functions. Lamb's constant relates to the size of the vortex cores and the rate at
Expand Down Expand Up @@ -335,7 +335,7 @@ def update_position(
)


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def collapsed_velocities_from_horseshoe_vortices(
points,
back_right_vortex_vertices,
Expand Down Expand Up @@ -419,7 +419,7 @@ def collapsed_velocities_from_horseshoe_vortices(
return induced_velocities


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def expanded_velocities_from_horseshoe_vortices(
points,
back_right_vortex_vertices,
Expand Down Expand Up @@ -503,7 +503,7 @@ def expanded_velocities_from_horseshoe_vortices(
return induced_velocities


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def collapsed_velocities_from_ring_vortices(
points,
back_right_vortex_vertices,
Expand Down Expand Up @@ -589,7 +589,7 @@ def collapsed_velocities_from_ring_vortices(
return induced_velocities


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def expanded_velocities_from_ring_vortices(
points,
back_right_vortex_vertices,
Expand Down Expand Up @@ -675,7 +675,7 @@ def expanded_velocities_from_ring_vortices(
return induced_velocities


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def collapsed_velocities_from_line_vortices(
points,
origins,
Expand Down Expand Up @@ -759,10 +759,10 @@ def collapsed_velocities_from_line_vortices(
r_0_z = termination[2] - origin[2]

# Find the r_0 vector's length.
r_0 = math.sqrt(r_0_x ** 2 + r_0_y ** 2 + r_0_z ** 2)
r_0 = math.sqrt(r_0_x**2 + r_0_y**2 + r_0_z**2)

c_1 = strength / (4 * math.pi)
c_2 = r_0 ** 2 * r_c ** 2
c_2 = r_0**2 * r_c**2

for point_id in range(num_points):
point = points[point_id]
Expand All @@ -783,24 +783,24 @@ def collapsed_velocities_from_line_vortices(
r_3_z = r_1_x * r_2_y - r_1_y * r_2_x

# Find the r_1, r_2, and r_3 vectors' lengths.
r_1 = math.sqrt(r_1_x ** 2 + r_1_y ** 2 + r_1_z ** 2)
r_2 = math.sqrt(r_2_x ** 2 + r_2_y ** 2 + r_2_z ** 2)
r_3 = math.sqrt(r_3_x ** 2 + r_3_y ** 2 + r_3_z ** 2)
r_1 = math.sqrt(r_1_x**2 + r_1_y**2 + r_1_z**2)
r_2 = math.sqrt(r_2_x**2 + r_2_y**2 + r_2_z**2)
r_3 = math.sqrt(r_3_x**2 + r_3_y**2 + r_3_z**2)

c_3 = r_1_x * r_2_x + r_1_y * r_2_y + r_1_z * r_2_z

# If part of the vortex is so close to the point that they are touching (
# within machine epsilon), there is a removable discontinuity. In this
# case, continue to the next point because there is no velocity induced
# by the current vortex at this point.
if r_1 < eps or r_2 < eps or r_3 ** 2 < eps:
if r_1 < eps or r_2 < eps or r_3**2 < eps:
continue
else:
c_4 = (
c_1
* (r_1 + r_2)
* (r_1 * r_2 - c_3)
/ (r_1 * r_2 * (r_3 ** 2 + c_2))
/ (r_1 * r_2 * (r_3**2 + c_2))
)
velocities[point_id, 0] += c_4 * r_3_x
velocities[point_id, 1] += c_4 * r_3_y
Expand All @@ -809,7 +809,7 @@ def collapsed_velocities_from_line_vortices(
return velocities


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def expanded_velocities_from_line_vortices(
points,
origins,
Expand Down Expand Up @@ -893,10 +893,10 @@ def expanded_velocities_from_line_vortices(
r_0_z = termination[2] - origin[2]

# Find the r_0 vector's length.
r_0 = math.sqrt(r_0_x ** 2 + r_0_y ** 2 + r_0_z ** 2)
r_0 = math.sqrt(r_0_x**2 + r_0_y**2 + r_0_z**2)

c_1 = strength / (4 * math.pi)
c_2 = r_0 ** 2 * r_c ** 2
c_2 = r_0**2 * r_c**2

for point_id in range(num_points):
point = points[point_id]
Expand All @@ -917,24 +917,24 @@ def expanded_velocities_from_line_vortices(
r_3_z = r_1_x * r_2_y - r_1_y * r_2_x

# Find the r_1, r_2, and r_3 vectors' lengths.
r_1 = math.sqrt(r_1_x ** 2 + r_1_y ** 2 + r_1_z ** 2)
r_2 = math.sqrt(r_2_x ** 2 + r_2_y ** 2 + r_2_z ** 2)
r_3 = math.sqrt(r_3_x ** 2 + r_3_y ** 2 + r_3_z ** 2)
r_1 = math.sqrt(r_1_x**2 + r_1_y**2 + r_1_z**2)
r_2 = math.sqrt(r_2_x**2 + r_2_y**2 + r_2_z**2)
r_3 = math.sqrt(r_3_x**2 + r_3_y**2 + r_3_z**2)

c_3 = r_1_x * r_2_x + r_1_y * r_2_y + r_1_z * r_2_z

# If part of the vortex is so close to the point that they are touching (
# within machine epsilon), there is a removable discontinuity. In this
# case, set the velocity components to their true values, which are 0.0
# meters per second.
if r_1 < eps or r_2 < eps or r_3 ** 2 < eps:
if r_1 < eps or r_2 < eps or r_3**2 < eps:
continue
else:
c_4 = (
c_1
* (r_1 + r_2)
* (r_1 * r_2 - c_3)
/ (r_1 * r_2 * (r_3 ** 2 + c_2))
/ (r_1 * r_2 * (r_3**2 + c_2))
)
velocities[point_id, vortex_id, 0] = c_4 * r_3_x
velocities[point_id, vortex_id, 1] = c_4 * r_3_y
Expand Down
10 changes: 5 additions & 5 deletions pterasoftware/functions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,19 +119,19 @@ def angle_axis_rotation_matrix(angle, axis, axis_already_normalized=False):
rotation_matrix = np.array(
[
[
cos_theta + u_x ** 2 * (1 - cos_theta),
cos_theta + u_x**2 * (1 - cos_theta),
u_x * u_y * (1 - cos_theta) - u_z * sin_theta,
u_x * u_z * (1 - cos_theta) + u_y * sin_theta,
],
[
u_y * u_x * (1 - cos_theta) + u_z * sin_theta,
cos_theta + u_y ** 2 * (1 - cos_theta),
cos_theta + u_y**2 * (1 - cos_theta),
u_y * u_z * (1 - cos_theta) - u_x * sin_theta,
],
[
u_z * u_x * (1 - cos_theta) - u_y * sin_theta,
u_z * u_y * (1 - cos_theta) + u_x * sin_theta,
cos_theta + u_z ** 2 * (1 - cos_theta),
cos_theta + u_z**2 * (1 - cos_theta),
],
]
)
Expand All @@ -140,7 +140,7 @@ def angle_axis_rotation_matrix(angle, axis, axis_already_normalized=False):
return rotation_matrix


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def numba_centroid_of_quadrilateral(
front_left_vertex, front_right_vertex, back_left_vertex, back_right_vertex
):
Expand Down Expand Up @@ -616,7 +616,7 @@ def calculate_steady_freestream_wing_influences(steady_solver):
)


@njit(cache=True, fastmath=True)
@njit(cache=True, fastmath=False)
def numba_1d_explicit_cross(vectors_1, vectors_2):
"""This function takes in two arrays, each of which contain N vectors of 3
components. The function then calculates and returns the cross product of the two
Expand Down

0 comments on commit f17fb1b

Please sign in to comment.