Skip to content

Commit

Permalink
Merge pull request #1119 from jguarato/fix/plot-modes
Browse files Browse the repository at this point in the history 
Fix plotting of modes
  • Loading branch information
@raphaeltimbo
raphaeltimbo authored Oct 17, 2024
2 parents 92b522f + 411b649 commit d939735
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 84 deletions.
53 changes: 8 additions & 45 deletions ross/results.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,7 @@ def __init__(self, *, node, node_pos, ru_e, rv_e):
# kappa encodes the relation between the axis and the precession.
minor = np.sqrt(lam.min())
major = np.sqrt(lam.max())
kappa = minor / major

diff = nv - nu

# we need to evaluate if 0 < nv - nu < pi.
Expand All @@ -289,10 +289,11 @@ def __init__(self, *, node, node_pos, ru_e, rv_e):
# if nv = nu or nv = nu + pi then the response is a straight line.
if diff == 0 or diff == np.pi:
kappa = 0

# if 0 < nv - nu < pi, then a backward rotating mode exists.
elif 0 < diff < np.pi:
kappa *= -1
kappa = -minor / major
else:
kappa = minor / major

self.minor_axis = np.real(minor)
self.major_axis = np.real(major)
Expand Down Expand Up @@ -892,7 +893,6 @@ def __init__(
self.number_dof = number_dof
self.modes = self.evectors[: self.ndof]
self.shapes = []
kappa_modes = []
for mode in range(len(self.wn)):
self.shapes.append(
Shape(
Expand All @@ -904,42 +904,6 @@ def __init__(
number_dof=self.number_dof,
)
)
kappa_color = []
kappa_mode = self.kappa_mode(mode)
for kappa in kappa_mode:
kappa_color.append("blue" if kappa > 0 else "red")
kappa_modes.append(kappa_color)
self.kappa_modes = kappa_modes

@staticmethod
def whirl(kappa_mode):
"""Evaluate the whirl of a mode.
Parameters
----------
kappa_mode : list
A list with the value of kappa for each node related
to the mode/natural frequency of interest.
Returns
-------
whirldir : str
A string indicating the direction of precession related to the
kappa_mode.
Example
-------
>>> kappa_mode = [-5.06e-13, -3.09e-13, -2.91e-13, 0.011, -4.03e-13, -2.72e-13, -2.72e-13]
>>> ModalResults.whirl(kappa_mode)
'Forward'
"""
if all(kappa >= -1e-3 for kappa in kappa_mode):
whirldir = "Forward"
elif all(kappa <= 1e-3 for kappa in kappa_mode):
whirldir = "Backward"
else:
whirldir = "Mixed"
return whirldir

@staticmethod
@np.vectorize
Expand Down Expand Up @@ -1048,11 +1012,10 @@ def whirl_direction(self):
-------
whirl_w : array
An array of strings indicating the direction of precession related
to the kappa_mode. Backward, Mixed or Forward depending on values
of kappa_mode.
to the mode shape. Backward, Mixed or Forward.
"""
# whirl direction/values are methods because they are expensive.
whirl_w = [self.whirl(self.kappa_mode(wd)) for wd in range(len(self.wd))]
whirl_w = [self.shapes[wd].whirl for wd in range(len(self.wd))]

return np.array(whirl_w)

Expand Down Expand Up @@ -1211,7 +1174,7 @@ def plot_mode_3d(
f"{title}<br>"
f"Mode {mode} | "
f"{frequency['speed']} {frequency_units} | "
f"whirl: {self.whirl_direction()[mode]} | "
f"whirl: {shape.whirl} | "
f"{frequency[frequency_type]} {frequency_units} | "
f"{damping_name} = {damping_value:.2f}"
),
Expand Down Expand Up @@ -1306,7 +1269,7 @@ def plot_mode_2d(
f"{title}<br>"
f"Mode {mode} | "
f"{frequency['speed']} {frequency_units} | "
f"whirl: {self.whirl_direction()[mode]} | "
f"whirl: {shape.whirl} | "
f"{frequency[frequency_type]} {frequency_units} | "
f"{damping_name} = {damping_value:.2f}"
),
Expand Down
65 changes: 30 additions & 35 deletions ross/rotor_assembly.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,8 +373,6 @@ def flatten(l):
Ip_dsk = np.sum([disk.Ip for disk in self.disk_elements])
self.Ip = Ip_sh + Ip_dsk

self._v0 = None # used to call eigs

# number of dofs
half_ndof = self.number_dof / 2
self.ndof = int(
Expand Down Expand Up @@ -739,7 +737,7 @@ def run_modal(self, speed, num_modes=12, sparse=True, synchronous=False):
sparse : bool, optional
If True, ARPACK is used to calculate a desired number (according to
num_modes) or eigenvalues and eigenvectors.
If False, scipy.linalg.eig() is used to calculate all the eigenvalues and
If False, `scipy.linalg.eig()` is used to calculate all the eigenvalues and
eigenvectors.
Default is True.
synchronous : bool, optional
Expand Down Expand Up @@ -1400,7 +1398,7 @@ def _index(eigenvalues):
positive = [i for i in ind[len(a) // 2 :]]
negative = [i for i in ind[: len(a) // 2]]

idx = np.array([positive, negative]).flatten()
idx = np.array([*positive, *negative])

return idx

Expand All @@ -1412,7 +1410,7 @@ def _eigen(
frequency=None,
sorted_=True,
A=None,
sparse=True,
sparse=None,
synchronous=False,
):
"""Calculate eigenvalues and eigenvectors.
Expand All @@ -1436,14 +1434,16 @@ def _eigen(
frequency: float, pint.Quantity
Excitation frequency. Default units is rad/s.
sorted_ : bool, optional
Sort considering the imaginary part (wd)
Default is True
Sort considering the imaginary part (wd).
Default is True.
A : np.array, optional
Matrix for which eig will be calculated.
Defaul is the rotor A matrix.
sparse : bool, optional
If sparse, eigenvalues will be calculated with arpack.
Default is True.
If True, eigenvalues are computed using ARPACK. If False, they are
computed with `scipy.linalg.eig()`. When sparse is False, eigenvalues
are filtered to exclude rigid body modes. If sparse is None, no filtering
is applied. Default is None.
synchronous : bool, optional
If True a synchronous analysis is carried out.
Default is False.
Expand All @@ -1465,44 +1465,42 @@ def _eigen(
if A is None:
A = self.A(speed=speed, frequency=frequency, synchronous=synchronous)

filter_eigenpairs = lambda values, vectors, indices: (
values[indices],
vectors[:, indices],
)

if synchronous:
evalues, evectors = la.eig(A)

idx = np.where(np.imag(evalues) != 0)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
evalues, evectors = filter_eigenpairs(evalues, evectors, idx)
idx = np.where(np.abs(np.real(evalues) / np.imag(evalues)) < 1000)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
evalues, evectors = filter_eigenpairs(evalues, evectors, idx)
else:
if sparse is True:
if sparse:
try:
evalues, evectors = las.eigs(
A,
k=2 * num_modes,
k=min(2 * num_modes, max(num_modes, A.shape[0] - 2)),
sigma=1,
ncv=4 * num_modes,
which="LM",
v0=self._v0,
v0=np.ones(A.shape[0]),
)
# store v0 as a linear combination of the previously
# calculated eigenvectors to use in the next call to eigs
self._v0 = np.real(sum(evectors.T))

# Disregard rigid body modes:
idx = np.where(np.abs(evalues) > 0.1)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
except las.ArpackError:
evalues, evectors = la.eig(A)
else:
evalues, evectors = la.eig(A)

if sorted_ is False:
return evalues, evectors
if sparse is not None:
idx = np.where((np.imag(evalues) != 0) & (np.abs(evalues) > 0.1))[0]
evalues, evectors = filter_eigenpairs(evalues, evectors, idx)

idx = self._index(evalues)
if sorted_:
idx = self._index(evalues)
evalues, evectors = filter_eigenpairs(evalues, evectors, idx)

return evalues[idx], evectors[:, idx]
return evalues, evectors

def _lti(self, speed, frequency=None):
"""Continuous-time linear time invariant system.
Expand Down Expand Up @@ -1653,8 +1651,7 @@ def transfer_matrix(self, speed=None, frequency=None, modes=None):

# calculate eigenvalues and eigenvectors using la.eig to get
# left and right eigenvectors.

evals, psi = self._eigen(speed=speed, frequency=frequency, sparse=False)
evals, psi = self._eigen(speed=speed, frequency=frequency)

psi_inv = la.inv(psi)

Expand Down Expand Up @@ -2002,8 +1999,8 @@ def run_unbalance_response(
Unbalance magnitude (kg.m).
unbalance_phase : list, float, pint.Quantity
Unbalance phase (rad).
frequency : list, float, pint.Quantity
Array with the desired range of frequencies (rad/s).
frequency : list, pint.Quantity
List with the desired range of frequencies (rad/s).
modes : list, optional
Modes that will be used to calculate the frequency response
(all modes will be used if a list is not given).
Expand Down Expand Up @@ -4017,8 +4014,6 @@ def __init__(
Ip_dsk = np.sum([disk.Ip for disk in self.disk_elements])
self.Ip = Ip_sh + Ip_dsk

self._v0 = None # used to call eigs

# number of dofs
self.ndof = int(
4 * max([el.n for el in shaft_elements])
Expand Down
8 changes: 4 additions & 4 deletions ross/tests/test_rotor_assembly.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1835,11 +1835,11 @@ def rotor_6dof():

def test_modal_6dof(rotor_6dof):
modal = rotor_6dof.run_modal(speed=0, sparse=False)
wn = np.array([0.0, 47.62138, 91.79647, 96.28891, 274.56591, 296.5005])
wd = np.array([0.01079, 47.62156, 91.79656, 96.289, 274.56607, 296.50068])
wn = np.array([47.62138, 91.79647, 96.28891, 274.56591, 296.5005])
wd = np.array([47.62156, 91.79656, 96.289, 274.56607, 296.50068])

assert_almost_equal(modal.wn[:6], wn, decimal=2)
assert_almost_equal(modal.wd[:6], wd, decimal=2)
assert_almost_equal(modal.wn[:5], wn, decimal=2)
assert_almost_equal(modal.wd[:5], wd, decimal=2)


@pytest.fixture
Expand Down

0 comments on commit d939735

Please sign in to comment.