diff --git a/pypk/api_core.py b/pypk/api_core.py
index 4787eb10579d712050fa867719dad899a1688cf0..3d61916cd5e28ae8d1dee3d5acc33d148277c428 100644
--- a/pypk/api_core.py
+++ b/pypk/api_core.py
@@ -41,7 +41,7 @@ def init_pypk(cfg):
agg = KSAggregator(cfg['rho_ks'])
# p-k method
if cfg['method'] == 'pk':
- sol = Pk(name, flu, eig, agg, cfg['mach'], cfg['k_ref'], cfg['l_ref'], cfg['vrb'], cfg['n_modes'], n_speed, tol=1e-3)
+ sol = Pk(name, flu, eig, agg, cfg['mach'], cfg['k_ref'], cfg['l_ref'], cfg['vrb'], cfg['n_modes'], n_speed, tol=1e-3, max_it=10)
elif cfg['method'] == 'nipk':
sol = Nipk(name, flu, eig, agg, cfg['mach'], cfg['k_ref'], cfg['l_ref'], cfg['vrb'], cfg['n_modes'], n_speed)
else:
diff --git a/pypk/eigensolver.py b/pypk/eigensolver.py
index 9ff0e6e787ce202567799d10ab85244b85e5654a..6c2872158884f39021a5a432a937c2f5e0bbfc51 100644
--- a/pypk/eigensolver.py
+++ b/pypk/eigensolver.py
@@ -25,22 +25,31 @@ class EigenSolver:
self._n = n # number of modes
self._g = (1 + 1j * g) # structural damping (complex proportional stiffness)
- def compute(self, rho, u, m, k, q):
+ def compute_natural(self, m, k):
+ """Compute natural frequencies and mode shapes
+ m : mass matrix
+ k : stiffness matrix
+ """
+ l, v = spla.eig(k, m) # det(-omega^2 * M + K) = 0
+ return np.sqrt(np.real(l)), v
+
+ def compute(self, rho, u, m, k, q, v0=None):
"""Compute eigenvalues and eigenvectors
rho : fuid density
u : freestream velocity
m : mass matrix
k : stiffness matrix
q : aerodynamic matrix
+ v0 : converged mode shapes from previous iteration (airspeed)
"""
p, v = spla.eig(self._g * k - 0.5 * rho * u**2 * q, -u**2 / self._l**2 * m) # det((p*u/l)^2*M + K - q_inf*Q) = 0
p = np.sqrt(p)
# Sort
p = np.where(np.imag(p) < 0, -p, p) # change sign of p if imag(p) < 0
- srt = np.imag(p).argsort()
- # Normalize (optional since scipy does it)
- for j in range(self._n):
- v[:,j] /= spla.norm(v[:,j])
+ if v0 is None:
+ srt = np.imag(p).argsort()
+ else:
+ srt = self._sort(v, v0)
return p[srt].T, v[:,srt]
def compute_gradients(self, rho, u, m, k, q, p, v, dm, dk, dq):
@@ -66,3 +75,16 @@ class EigenSolver:
b[-1,0] = 0
x = spla.lu_solve((spla.lu_factor(a)), b)
return x[0,0]
+
+ def _sort(self, v, v0):
+ """Sort current modes according to their correlation to previous modes
+ v : mode shapes from current iteration
+ v0 : converged mode shapes from previous iteration
+ """
+ # Compute correlation matrix
+ mac = np.zeros((self._n, self._n))
+ for i in range(self._n):
+ for j in range(self._n):
+ mac[i,j] = np.linalg.norm(np.dot(v[:,i], np.conjugate(v0[:,j]))) / (np.linalg.norm(v[:,i]) * np.linalg.norm(v0[:,j]))
+ # Find maximum for each mode
+ return np.argmax(mac, axis=0)
diff --git a/pypk/nipk.py b/pypk/nipk.py
index c6e875996d848ecdb42ced55046f3c7ce29dd3b5..c097f7407a7afba49341c14483010ddac45c186a 100644
--- a/pypk/nipk.py
+++ b/pypk/nipk.py
@@ -33,21 +33,21 @@ class Nipk(Solution):
def compute(self):
"""Compute frequencies and dampings
"""
+ omega_0, v_0 = self._eig.compute_natural(self._m, self._k) # natural frequencies and modes
self.rho, self.speed = self._flu.compute(self._mach)
for ivel in range(self.n_speed):
# Compute eigenvalues at zero airspeed
if self.speed[ivel] == 0:
- for imod in range(self.n_modes):
- self.freq[ivel, imod] = 0.
+ self.freq[ivel, :] = omega_0
continue
# Compute eigenvalues for each reference reduced frequency
p_k = np.zeros((self.n_freqs, self.n_modes), dtype=complex)
v_k = np.zeros((self.n_freqs, self.n_modes, self.n_modes), dtype=complex)
for ifrq in range(self.n_freqs):
- p_k[ifrq,:], v_k[ifrq,:,:] = self._eig.compute(self.rho[ivel], self.speed[ivel], self._m, self._k, self._q[ifrq,:,:])
+ p_k[ifrq,:], v_k[ifrq,:,:] = self._eig.compute(self.rho[ivel], self.speed[ivel], self._m, self._k, self._q[ifrq,:,:], v_0)
# Match reduced frequency by linear interpolation for each mode
for imod in range(self.n_modes):
- delta, idx = self.__find_index(self.speed[ivel], p_k[:, imod], self._kref)
+ delta, idx = self.__find_index(self.speed[ivel], imod, p_k[:, imod], self._kref)
p_re, p_im = self.__interp(self._kref[idx], delta[idx], p_k[idx, imod])
# Compute frequency and damping
self.freq[ivel, imod] = p_im * self.speed[ivel] / self._lref / (2 * np.pi)
@@ -59,6 +59,11 @@ class Nipk(Solution):
self._pk[ivel, imod, :] = p_k[idx, imod]
self._vk[ivel, imod, :, :] = v_k[idx, :, imod]
self._p[ivel, imod] = p_re + 1j * p_im
+ # Interpolate mode shapes for mode tracking
+ v_0 = np.zeros((self.n_modes, self.n_modes), dtype=complex)
+ for imod in range(self.n_modes):
+ a = (self._p[ivel, imod].imag - self._kk[ivel, imod, 0]) / (self._kk[ivel, imod, 1] - self._kk[ivel, imod, 0])
+ v_0[:, imod] = (1 - a) * self._vk[ivel, imod, 0, :] + a * self._vk[ivel, imod, 1, :]
# Compute aggregated damping over modes then over dynamic pressures
self._ksm = np.zeros(self.n_speed)
for i in range(self.n_speed):
@@ -87,9 +92,10 @@ class Nipk(Solution):
dq[np.unravel_index(i, dq.shape)] = 1.
self.damp_qim[i] = self.__compute_dg(z, z, zk+1j*dq)
- def __find_index(self, u, p, k):
+ def __find_index(self, u, m, p, k):
"""Find the indices between which the frequency must be interpolated
u : freestream velocity
+ m : mode number
p : eigenvalue solutions
k : reference reduced frequencies
"""
@@ -109,7 +115,7 @@ class Nipk(Solution):
elif ai == len(delta) - 1:
idx = [ai-1, ai]
else:
- raise RuntimeError(f'NIPK: could not match frequency for mode {i} at velocity {u}!\n')
+ raise RuntimeError(f'NIPK: could not match frequency for mode {m} at velocity {u}!\n')
return delta, idx
def __interp(self, k, delta, p):
diff --git a/pypk/pk.py b/pypk/pk.py
index 695ba7f4c12401625407f40fe1073e9f4479fbd7..803cfd1b71b2b92b985f34c8a61a3f83b82cd0f4 100644
--- a/pypk/pk.py
+++ b/pypk/pk.py
@@ -17,41 +17,48 @@
from .solution import Solution
import numpy as np
from scipy import interpolate
-import scipy.linalg as spla
class Pk(Solution):
"""Standard p-k method for flutter solution
"""
- def __init__(self, name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed, tol=1e-3):
+ def __init__(self, name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed, tol=1e-3, max_it=10):
super().__init__(name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed)
self.name = 'PK' # solver name
self._tol = tol # tolerance
+ self._max_it = max_it # maximum number of iterations
def compute(self):
"""Compute frequencies and dampings
"""
- omega = self.__sfreq(self._m, self._k) # natural frequencies
+ omega_0, v_0 = self._eig.compute_natural(self._m, self._k) # natural frequencies and modes
self.rho, self.speed = self._flu.compute(self._mach)
for ivel in range(self.n_speed):
# Compute eigenvalues at zero airspeed
if self.speed[ivel] == 0:
- for imod in range(self.n_modes):
- self.freq[ivel, imod] = 0.
+ self.freq[ivel, :] = omega_0
continue
+ v_tmp = np.zeros((self.n_modes, self.n_modes), dtype=complex) # temporary container for modes
for imod in range(self.n_modes):
- k = omega[imod] * self._lref / self.speed[ivel] # guess reduced frequency
+ k = omega_0[imod] * self._lref / self.speed[ivel] # guess reduced frequency
+ it = 0
while True:
- # Interpolate loads and solve eigenvalue problem
+ # Interpolate loads and solve eigenvalue problem
q = self.__interp(k, self._kref, self._q)
- p, _ = self._eig.compute(self.rho[ivel], self.speed[ivel], self._m, self._k, q)
+ p, v = self._eig.compute(self.rho[ivel], self.speed[ivel], self._m, self._k, q, v_0)
p = p[imod]
# Check if reduced frequency is converged
- if abs(np.imag(p) - k) < self._tol:
+ if abs(np.imag(p) - k) < self._tol or it == self._max_it:
+ v_tmp[:, imod] = v[:, imod]
self.freq[ivel, imod] = np.imag(p) * self.speed[ivel] / self._lref / (2 * np.pi)
self.damp[ivel, imod] = np.real(p) / np.imag(p)
+ if it == self._max_it:
+ print(f'PK: could not converge for mode {imod} at velocity {self.speed[ivel]}, using k={np.imag(p)} with tolerance {abs(np.imag(p) - k)}.')
break
else:
k = np.imag(p)
+ it += 1
+ # Update modes for mode tracking
+ v_0 = v_tmp
# Compute aggregated damping over modes then over dynamic pressures
self._ksm = np.zeros(self.n_speed)
for i in range(self.n_speed):
@@ -63,15 +70,6 @@ class Pk(Solution):
"""
raise NotImplementedError()
- def __sfreq(self, m, k):
- """Compute the natural frequencies of the structural model
- TODO: should it be given directly by the structural solver?
- m : mass matrix
- k : stiffness matrix
- """
- om, _ = spla.eig(k, m) # det(-omega^2 * M + K) = 0
- return np.sqrt(np.real(om))
-
def __interp(self, k, k_ref, q_ref):
"""Match the frequency and interpolate the eigenvalues and eigenvectors
k : interpolation reduced frequency