Add mode tracking. Add max iteration number for pk method. Bugfix.

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:

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)

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):

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