Skip to content
Snippets Groups Projects
Commit 4e58da94 authored by Adrien Crovato's avatar Adrien Crovato
Browse files

Merge branch 'adri' into 'master' 

Vesion 1.1

See merge request !1
parents 90e1daf6 be888a37
Branches master
Tags v1.1.0
1 merge request!1Vesion 1.1
Pipeline #50489 passed
Stage: build
    Stage: test
      Hide whitespace changes
      Inline Side-by-side
      .gitlab-ci.yml
      + 1
      − 1
      View file @ 4e58da94
      default: default:
      image: rboman/waves-py3:2020.3 image: rboman/waves-py3:2022.0
      .global_tag: &global_tag_def .global_tag: &global_tag_def
      tags: tags:
      ......
      pypk/__init__.py
      + 1
      − 1
      View file @ 4e58da94
      ...@@ -15,6 +15,6 @@ ...@@ -15,6 +15,6 @@
      # See the License for the specific language governing permissions and # See the License for the specific language governing permissions and
      # limitations under the License. # limitations under the License.
      __version__ = '1.0.0' __version__ = '1.1.0'
      from .api_core import init_pypk from .api_core import init_pypk
      pypk/api_core.py
      + 6
      − 3
      View file @ 4e58da94
      ...@@ -23,6 +23,8 @@ from .nipk import Nipk ...@@ -23,6 +23,8 @@ from .nipk import Nipk
      def init_pypk(cfg): def init_pypk(cfg):
      """Initialize PyPK """Initialize PyPK
      """ """
      # Case name
      name = cfg['name'] if 'name' in cfg else 'noname'
      # Fluid type # Fluid type
      if cfg['fluid'] == 'unmatched': if cfg['fluid'] == 'unmatched':
      n_speed = len(cfg['u_idx']) # number of airspeed test points n_speed = len(cfg['u_idx']) # number of airspeed test points
      ...@@ -33,14 +35,15 @@ def init_pypk(cfg): ...@@ -33,14 +35,15 @@ def init_pypk(cfg):
      else: else:
      raise RuntimeError('"fluid" entry must be "unmatched" or "matched_isa"!\n') raise RuntimeError('"fluid" entry must be "unmatched" or "matched_isa"!\n')
      # Eigen solver # Eigen solver
      eig = EigenSolver(cfg['l_ref'], cfg['n_modes']) g = cfg['g_struct'] if 'g_struct' in cfg else 0.0
      eig = EigenSolver(cfg['l_ref'], cfg['n_modes'], g)
      # Aggregator # Aggregator
      agg = KSAggregator(cfg['rho_ks']) agg = KSAggregator(cfg['rho_ks'])
      # p-k method # p-k method
      if cfg['method'] == 'pk': if cfg['method'] == 'pk':
      sol = Pk(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': elif cfg['method'] == 'nipk':
      sol = Nipk(flu, eig, agg, cfg['mach'], cfg['k_ref'], cfg['l_ref'], cfg['vrb'], cfg['n_modes'], n_speed) sol = Nipk(name, flu, eig, agg, cfg['mach'], cfg['k_ref'], cfg['l_ref'], cfg['vrb'], cfg['n_modes'], n_speed)
      else: else:
      raise RuntimeError('"method" entry must be "pk" or "nipk"!\n') raise RuntimeError('"method" entry must be "pk" or "nipk"!\n')
      return sol return sol
      pypk/api_om.py
      + 13
      − 3
      View file @ 4e58da94
      ...@@ -23,11 +23,13 @@ class PypkSolver(om.ExplicitComponent): ...@@ -23,11 +23,13 @@ class PypkSolver(om.ExplicitComponent):
      """OpenMDAO flutter solution wrapper """OpenMDAO flutter solution wrapper
      """ """
      def initialize(self): def initialize(self):
      self.options.declare('scn_name', desc='name of current scenario')
      self.options.declare('nm', desc='number of modes') self.options.declare('nm', desc='number of modes')
      self.options.declare('nu', desc='number of airspeed test points') self.options.declare('nu', desc='number of airspeed test points')
      self.options.declare('sol', desc='PyPK solver', recordable=False) self.options.declare('sol', desc='PyPK solver', recordable=False)
      def setup(self): def setup(self):
      self.scn_name = self.options['scn_name']
      self.nm = self.options['nm'] self.nm = self.options['nm']
      self.nu = self.options['nu'] self.nu = self.options['nu']
      self.sol = self.options['sol'] self.sol = self.options['sol']
      ...@@ -41,18 +43,26 @@ class PypkSolver(om.ExplicitComponent): ...@@ -41,18 +43,26 @@ class PypkSolver(om.ExplicitComponent):
      self.add_output('f_speed', val=1., desc='flutter speed') self.add_output('f_speed', val=1., desc='flutter speed')
      self.add_output('f_mode', val=1, desc='flutter mode') self.add_output('f_mode', val=1, desc='flutter mode')
      # Partials # Partials
      # self.declare_partials(of='g_ks', wrt=['M', 'K', 'Q_re', 'Q_im'], method='exact') self.declare_partials(of='g_agg', wrt=['M', 'K', 'Q_re', 'Q_im'], method='exact')
      def compute(self, inputs, outputs): def compute(self, inputs, outputs):
      self.sol.set_matrices(inputs['M'], inputs['K'], inputs['Q_re'] + 1j * inputs['Q_im']) self.sol.set_matrices(inputs['M'], inputs['K'], inputs['Q_re'] + 1j * inputs['Q_im'])
      self.sol.compute() self.sol.compute()
      self.sol.find_flutter() self.sol.find_flutter()
      self.sol.save('_' + self.scn_name)
      outputs['f'] = self.sol.freq outputs['f'] = self.sol.freq
      outputs['g'] = self.sol.damp outputs['g'] = self.sol.damp
      outputs['g_agg'] = self.sol.damp_agg outputs['g_agg'] = self.sol.damp_agg
      outputs['f_speed'] = self.sol.f_speed outputs['f_speed'] = self.sol.f_speed
      outputs['f_mode'] = self.sol.f_mode outputs['f_mode'] = self.sol.f_mode
      def compute_partials(self, inputs, partials):
      self.sol.compute_gradients()
      partials['g_agg', 'M'] = self.sol.damp_m
      partials['g_agg', 'K'] = self.sol.damp_k
      partials['g_agg', 'Q_re'] = self.sol.damp_qre
      partials['g_agg', 'Q_im'] = self.sol.damp_qim
      class PypkBuilder(Builder): class PypkBuilder(Builder):
      """PyPK builder for OpenMDAO """PyPK builder for OpenMDAO
      """ """
      ...@@ -66,7 +76,7 @@ class PypkBuilder(Builder): ...@@ -66,7 +76,7 @@ class PypkBuilder(Builder):
      """ """
      self.__sol = init_pypk(cfg) self.__sol = init_pypk(cfg)
      def get_solver(self): def get_solver(self, scenario_name=''):
      """Return OpenMDAO component containing the solver """Return OpenMDAO component containing the solver
      """ """
      return PypkSolver(nm=self.__sol.n_modes, nu=self.__sol.n_speed, sol=self.__sol) return PypkSolver(scn_name=scenario_name, nm=self.__sol.n_modes, nu=self.__sol.n_speed, sol=self.__sol)
      pypk/eigensolver.py
      + 32
      − 9
      View file @ 4e58da94
      ...@@ -20,26 +20,36 @@ import scipy.linalg as spla ...@@ -20,26 +20,36 @@ import scipy.linalg as spla
      class EigenSolver: class EigenSolver:
      """Eigen solver """Eigen solver
      """ """
      def __init__(self, l, n): def __init__(self, l, n, g=0):
      self._l = l # reference length self._l = l # reference length
      self._n = n # number of modes 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 """Compute eigenvalues and eigenvectors
      rho : fuid density rho : fuid density
      u : freestream velocity u : freestream velocity
      m : mass matrix m : mass matrix
      k : stiffness matrix k : stiffness matrix
      q : aerodynamic matrix q : aerodynamic matrix
      v0 : converged mode shapes from previous iteration (airspeed)
      """ """
      p, v = spla.eig(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, 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) p = np.sqrt(p)
      # Sort # Sort
      p = np.where(np.imag(p) < 0, -p, p) # change sign of p if imag(p) < 0 p = np.where(np.imag(p) < 0, -p, p) # change sign of p if imag(p) < 0
      srt = np.imag(p).argsort() if v0 is None:
      # Normalize (optional since scipy does it) srt = np.imag(p).argsort()
      for j in range(self._n): else:
      v[:,j] /= spla.norm(v[:,j]) srt = self._sort(v, v0)
      return p[srt].T, v[:,srt] return p[srt].T, v[:,srt]
      def compute_gradients(self, rho, u, m, k, q, p, v, dm, dk, dq): def compute_gradients(self, rho, u, m, k, q, p, v, dm, dk, dq):
      ...@@ -58,10 +68,23 @@ class EigenSolver: ...@@ -58,10 +68,23 @@ class EigenSolver:
      a = np.zeros((self._n + 1, self._n + 1), dtype=complex) a = np.zeros((self._n + 1, self._n + 1), dtype=complex)
      b = np.zeros((self._n + 1, 1), dtype=complex) b = np.zeros((self._n + 1, 1), dtype=complex)
      a[:-1,0] = -2 * p * u**2 / self._l**2 * m @ v a[:-1,0] = -2 * p * u**2 / self._l**2 * m @ v
      a[:-1,1:] = p**2 * u**2 / self._l**2 * m + k - 0.5 * rho * u**2 * q a[:-1,1:] = p**2 * u**2 / self._l**2 * m + self._g * k - 0.5 * rho * u**2 * q
      a[-1,0] = 0 a[-1,0] = 0
      a[-1,1:] = v.T a[-1,1:] = v.T
      b[:-1,0] = (p**2 * u**2 / self._l**2 * dm + dk - 0.5 * rho * u**2 * dq) @ v b[:-1,0] = (p**2 * u**2 / self._l**2 * dm + self._g * dk - 0.5 * rho * u**2 * dq) @ v
      b[-1,0] = 0 b[-1,0] = 0
      x = spla.lu_solve((spla.lu_factor(a)), b) x = spla.lu_solve((spla.lu_factor(a)), b)
      return x[0,0] 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
      + 14
      − 8
      View file @ 4e58da94
      ...@@ -20,8 +20,8 @@ import numpy as np ...@@ -20,8 +20,8 @@ import numpy as np
      class Nipk(Solution): class Nipk(Solution):
      """Non-iterative p-k method for flutter solution """Non-iterative p-k method for flutter solution
      """ """
      def __init__(self, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed): def __init__(self, name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed):
      super().__init__(flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed) super().__init__(name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed)
      self.name = 'NIPK' # solver name self.name = 'NIPK' # solver name
      self._pk = np.zeros((self.n_speed, self.n_modes, 2), dtype=complex) # eigenvalues (at reference reduced frequencies) self._pk = np.zeros((self.n_speed, self.n_modes, 2), dtype=complex) # eigenvalues (at reference reduced frequencies)
      self._vk = np.zeros((self.n_speed, self.n_modes, 2, self.n_modes), dtype=complex) # eigenmodes (at reference reduced frequencies) self._vk = np.zeros((self.n_speed, self.n_modes, 2, self.n_modes), dtype=complex) # eigenmodes (at reference reduced frequencies)
      ...@@ -33,21 +33,21 @@ class Nipk(Solution): ...@@ -33,21 +33,21 @@ class Nipk(Solution):
      def compute(self): def compute(self):
      """Compute frequencies and dampings """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) self.rho, self.speed = self._flu.compute(self._mach)
      for ivel in range(self.n_speed): for ivel in range(self.n_speed):
      # Compute eigenvalues at zero airspeed # Compute eigenvalues at zero airspeed
      if self.speed[ivel] == 0: if self.speed[ivel] == 0:
      for imod in range(self.n_modes): self.freq[ivel, :] = omega_0
      self.freq[ivel, imod] = 0.
      continue continue
      # Compute eigenvalues for each reference reduced frequency # Compute eigenvalues for each reference reduced frequency
      p_k = np.zeros((self.n_freqs, self.n_modes), dtype=complex) 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) v_k = np.zeros((self.n_freqs, self.n_modes, self.n_modes), dtype=complex)
      for ifrq in range(self.n_freqs): 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 # Match reduced frequency by linear interpolation for each mode
      for imod in range(self.n_modes): 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]) p_re, p_im = self.__interp(self._kref[idx], delta[idx], p_k[idx, imod])
      # Compute frequency and damping # Compute frequency and damping
      self.freq[ivel, imod] = p_im * self.speed[ivel] / self._lref / (2 * np.pi) self.freq[ivel, imod] = p_im * self.speed[ivel] / self._lref / (2 * np.pi)
      ...@@ -59,6 +59,11 @@ class Nipk(Solution): ...@@ -59,6 +59,11 @@ class Nipk(Solution):
      self._pk[ivel, imod, :] = p_k[idx, imod] self._pk[ivel, imod, :] = p_k[idx, imod]
      self._vk[ivel, imod, :, :] = v_k[idx, :, imod] self._vk[ivel, imod, :, :] = v_k[idx, :, imod]
      self._p[ivel, imod] = p_re + 1j * p_im 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 # Compute aggregated damping over modes then over dynamic pressures
      self._ksm = np.zeros(self.n_speed) self._ksm = np.zeros(self.n_speed)
      for i in range(self.n_speed): for i in range(self.n_speed):
      ...@@ -87,9 +92,10 @@ class Nipk(Solution): ...@@ -87,9 +92,10 @@ class Nipk(Solution):
      dq[np.unravel_index(i, dq.shape)] = 1. dq[np.unravel_index(i, dq.shape)] = 1.
      self.damp_qim[i] = self.__compute_dg(z, z, zk+1j*dq) 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 """Find the indices between which the frequency must be interpolated
      u : freestream velocity u : freestream velocity
      m : mode number
      p : eigenvalue solutions p : eigenvalue solutions
      k : reference reduced frequencies k : reference reduced frequencies
      """ """
      ...@@ -109,7 +115,7 @@ class Nipk(Solution): ...@@ -109,7 +115,7 @@ class Nipk(Solution):
      elif ai == len(delta) - 1: elif ai == len(delta) - 1:
      idx = [ai-1, ai] idx = [ai-1, ai]
      else: 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 return delta, idx
      def __interp(self, k, delta, p): def __interp(self, k, delta, p):
      ......
      pypk/pk.py
      + 17
      − 19
      View file @ 4e58da94
      ...@@ -17,41 +17,48 @@ ...@@ -17,41 +17,48 @@
      from .solution import Solution from .solution import Solution
      import numpy as np import numpy as np
      from scipy import interpolate from scipy import interpolate
      import scipy.linalg as spla
      class Pk(Solution): class Pk(Solution):
      """Standard p-k method for flutter solution """Standard p-k method for flutter solution
      """ """
      def __init__(self, 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__(flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed) super().__init__(name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed)
      self.name = 'PK' # solver name self.name = 'PK' # solver name
      self._tol = tol # tolerance self._tol = tol # tolerance
      self._max_it = max_it # maximum number of iterations
      def compute(self): def compute(self):
      """Compute frequencies and dampings """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) self.rho, self.speed = self._flu.compute(self._mach)
      for ivel in range(self.n_speed): for ivel in range(self.n_speed):
      # Compute eigenvalues at zero airspeed # Compute eigenvalues at zero airspeed
      if self.speed[ivel] == 0: if self.speed[ivel] == 0:
      for imod in range(self.n_modes): self.freq[ivel, :] = omega_0
      self.freq[ivel, imod] = 0.
      continue continue
      v_tmp = np.zeros((self.n_modes, self.n_modes), dtype=complex) # temporary container for modes
      for imod in range(self.n_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: while True:
      # Interpolate loads and solve eigenvalue problem # Interpolate loads and solve eigenvalue problem
      q = self.__interp(k, self._kref, self._q) 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] p = p[imod]
      # Check if reduced frequency is converged # 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.freq[ivel, imod] = np.imag(p) * self.speed[ivel] / self._lref / (2 * np.pi)
      self.damp[ivel, imod] = np.real(p) / np.imag(p) 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 break
      else: else:
      k = np.imag(p) k = np.imag(p)
      it += 1
      # Update modes for mode tracking
      v_0 = v_tmp
      # Compute aggregated damping over modes then over dynamic pressures # Compute aggregated damping over modes then over dynamic pressures
      self._ksm = np.zeros(self.n_speed) self._ksm = np.zeros(self.n_speed)
      for i in range(self.n_speed): for i in range(self.n_speed):
      ...@@ -63,15 +70,6 @@ class Pk(Solution): ...@@ -63,15 +70,6 @@ class Pk(Solution):
      """ """
      raise NotImplementedError() 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): def __interp(self, k, k_ref, q_ref):
      """Match the frequency and interpolate the eigenvalues and eigenvectors """Match the frequency and interpolate the eigenvalues and eigenvectors
      k : interpolation reduced frequency k : interpolation reduced frequency
      ......
      pypk/solution.py
      + 11
      − 11
      View file @ 4e58da94
      ...@@ -19,8 +19,9 @@ import numpy as np ...@@ -19,8 +19,9 @@ import numpy as np
      class Solution: class Solution:
      """Base class for solving the aeroelastic equation """Base class for solving the aeroelastic equation
      """ """
      def __init__(self, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed): def __init__(self, name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed):
      self.name = None # solver name self.name = None # solver name
      self._name = name # case name
      self._flu = flu # fluid type self._flu = flu # fluid type
      self._eig = eig # eigen solver self._eig = eig # eigen solver
      self._agg = agg # aggregator self._agg = agg # aggregator
      ...@@ -61,15 +62,15 @@ class Solution: ...@@ -61,15 +62,15 @@ class Solution:
      """Find the flutter speed and display it """Find the flutter speed and display it
      """ """
      found = False found = False
      self.f_speed = np.inf self.f_speed = -1.
      self.f_mode = np.inf self.f_mode = -1
      for ivel in range(len(self.speed) - 1): for ivel in range(len(self.speed) - 1):
      if (found): if (found):
      break break
      for imod in range(self.n_modes): for imod in range(self.n_modes):
      if (self.damp[ivel, imod] > 0 or self.damp[ivel, imod] < 0 and self.damp[ivel+1, imod] > 0): if (self.damp[ivel, imod] > 0 or self.damp[ivel, imod] < 0 and self.damp[ivel+1, imod] > 0):
      if self._v > 0: if self._v > 0:
      print('{:s}: flutter found at Mach {:.2f} for mode #{:d} near u_inf = {:.0f}.'.format(self.name, self._mach, imod, self.speed[ivel])) print('{:s}: flutter found at Mach {:.2f} for mode #{:d} at {:.0f} Hz near u_inf = {:.0f}.'.format(self.name, self._mach, imod, self.freq[ivel, imod], self.speed[ivel]))
      self.f_speed = self.speed[ivel] self.f_speed = self.speed[ivel]
      self.f_mode = imod self.f_mode = imod
      found = True found = True
      ...@@ -77,21 +78,20 @@ class Solution: ...@@ -77,21 +78,20 @@ class Solution:
      if self._v > 0 and not found: if self._v > 0 and not found:
      print('{:s}: no flutter found at Mach {:.2f} for speed range [{:.0f}; {:.0f}].'.format(self.name, self._mach, self.speed[1], self.speed[-1])) print('{:s}: no flutter found at Mach {:.2f} for speed range [{:.0f}; {:.0f}].'.format(self.name, self._mach, self.speed[1], self.speed[-1]))
      def save(self, sfx): def save(self, sfx=''):
      """Write results to disk """Write results to disk
      sfx : filename suffix sfx : filename suffix
      """ """
      hdr = '{:>8s}, {:>10s}, {:>10s}, {:>10s}'.format('Mach', 'u_f', 'mode', 'AGG(g)') hdr = 'Mach={:.2f}, u_f={:.2f}, mode={:d}, AGG(g)={:.2f}\n'.format(self._mach, self.f_speed, self.f_mode, self.damp_agg)
      np.savetxt(f'results_{sfx}.dat', np.array([self._mach, self.f_speed, self.f_mode, self.damp_agg]).reshape(1,4), fmt='%10.4f', delimiter=', ', header=hdr) hdr += '{:>9s}'.format('u_inf')
      hdr = '{:>9s}'.format('u_inf')
      for j in range(self.n_modes): for j in range(self.n_modes):
      hdr += ', {:>11s}'.format(f'f_{j}') hdr += ', {:>11s}'.format(f'f_{j}')
      for j in range(self.n_modes): for j in range(self.n_modes):
      hdr += ', {:>11s}'.format(f'g_{j}') hdr += ', {:>11s}'.format(f'g_{j}')
      data = np.concatenate((self.speed.reshape(self.speed.shape[0], 1), self.freq, self.damp), axis=1) data = np.concatenate((self.speed.reshape(self.speed.shape[0], 1), self.freq, self.damp), axis=1)
      np.savetxt(f'vgf_{sfx}.dat', data, fmt='%+.4e', delimiter=', ', header=hdr) np.savetxt(f'vgf_{self._name}{sfx}.dat', data, fmt='%+.4e', delimiter=', ', header=hdr)
      def plot(self, sfx, show=True, format=''): def plot(self, sfx='', show=True, format=''):
      """Plot frequency and damping as a function of the velocity """Plot frequency and damping as a function of the velocity
      show : whether to display the plot show : whether to display the plot
      format : format to which the plot will be written to disk, if not empty format : format to which the plot will be written to disk, if not empty
      ...@@ -116,5 +116,5 @@ class Solution: ...@@ -116,5 +116,5 @@ class Solution:
      axs[1].spines['right'].set_visible(False) axs[1].spines['right'].set_visible(False)
      axs[1].plot(u, [0.] * self.n_speed, color='k', lw=1, ls='-.') # zero-damping line axs[1].plot(u, [0.] * self.n_speed, color='k', lw=1, ls='-.') # zero-damping line
      fig.tight_layout() fig.tight_layout()
      if format: plt.savefig(f'vgf_{sfx}.' + format) if format: plt.savefig(f'vgf_{self._name}{sfx}.' + format)
      if show: plt.show() if show: plt.show()
      setup.py
      + 1
      − 1
      View file @ 4e58da94
      ...@@ -31,6 +31,6 @@ setup( ...@@ -31,6 +31,6 @@ setup(
      url='https://gitlab.uliege.be/am-dept/pypk', url='https://gitlab.uliege.be/am-dept/pypk',
      license='Apache 2.0', license='Apache 2.0',
      packages=find_packages(include=['pypk*']), packages=find_packages(include=['pypk*']),
      install_requires=['numpy>=1.22', 'scipy>=1.8', 'matplotlib>=3.5'], install_requires=['numpy>=1.22,<2.0.0', 'scipy>=1.8', 'matplotlib>=3.5'],
      classifiers=['Operating System :: OS Independent', 'Programming Language :: Python'], classifiers=['Operating System :: OS Independent', 'Programming Language :: Python'],
      ) )
      tests/agard.py
      + 3
      − 2
      View file @ 4e58da94
      ...@@ -23,6 +23,7 @@ import numpy as np ...@@ -23,6 +23,7 @@ import numpy as np
      def get_cfg(): def get_cfg():
      cfg = { cfg = {
      'name': 'agard', # case name
      'k_ref': np.array([0.1, 0.3]), # reference reduced frequencies 'k_ref': np.array([0.1, 0.3]), # reference reduced frequencies
      'l_ref': 0.5 * 0.559, # reference length (half root chord) 'l_ref': 0.5 * 0.559, # reference length (half root chord)
      'mach': 0.901, # freesteam Mach number 'mach': 0.901, # freesteam Mach number
      ...@@ -60,8 +61,8 @@ def compute(sol): ...@@ -60,8 +61,8 @@ def compute(sol):
      sol.compute() sol.compute()
      sol.find_flutter() sol.find_flutter()
      # Save results # Save results
      sol.save('agard_' + sol.name.lower()) sol.save('_' + sol.name.lower())
      sol.plot('agard_' + sol.name.lower(), show=False, format='.pdf') sol.plot('_' + sol.name.lower(), show=False, format='pdf')
      def main(): def main():
      # Set up p-k and non-iterative p-k method # Set up p-k and non-iterative p-k method
      ......
      tests/agard_grad.py 0 → 100644
      + 95
      − 0
      View file @ 4e58da94
      #!/usr/bin/env python3
      # -*- coding: utf8 -*-
      # test encoding: à-é-è-ô-ï-€
      # Copyright 2024 University of Liège
      #
      # Licensed under the Apache License, Version 2.0 (the "License");
      # you may not use this file except in compliance with the License.
      # You may obtain a copy of the License at
      #
      # http://www.apache.org/licenses/LICENSE-2.0
      #
      # Unless required by applicable law or agreed to in writing, software
      # distributed under the License is distributed on an "AS IS" BASIS,
      # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      # See the License for the specific language governing permissions and
      # limitations under the License.
      # Compute flutter of AGARD 445.6 wing
      from pypk import init_pypk
      import numpy as np
      def get_cfg():
      cfg = {
      'name': 'agard', # case name
      'k_ref': np.array([0.1, 0.3]), # reference reduced frequencies
      'l_ref': 0.5 * 0.559, # reference length (half root chord)
      'mach': 0.901, # freesteam Mach number
      'rho_inf': 0.0995, # freesteam density
      'u_idx': np.linspace(0.2, 0.5, 31, endpoint=True), # velocity index range
      'mu': 143.920, # mass ratio
      'f_ref': 40.3511, # frequency of first torsional mode
      'n_modes': 4, # number of modes
      'g_struct': 0.02, # structural damping
      'rho_ks': 1e6, # KS aggregation parameter for KS
      'method': 'nipk', # method type
      'fluid': 'unmatched', # fluid type
      'vrb': 1 # verbosity level
      }
      return cfg
      def get_matrices():
      # Structural mass matrix
      m = np.diag([1.0, 1.0, 1.0, 1.0])
      # Structural stiffness matrix
      k = np.diag([3695.2, 63392.2, 99892.1, 368026.0])
      # Generalized aerodynamic force matrices (at k=0.1 and k=0.3)
      q = np.array([[[-7.747e-01-2.058e-01j, -6.014e+00-1.949e-01j, 3.962e+00+6.734e-02j, -3.403e+00-1.784e-01j],
      [ 9.520e-01-1.703e-02j, 4.900e+00-1.707e+00j, -4.416e+00+1.048e+00j, -5.951e+00-3.470e-01j],
      [-6.597e-02+6.036e-02j, 3.580e-01+5.907e-01j, -1.046e+00-6.893e-01j, -8.272e-02-6.273e-01j],
      [-1.238e-01+8.540e-02j, 3.588e-01+5.434e-01j, 8.221e-01-4.957e-01j, 2.922e+00-9.594e-01j]],
      [[-6.560e-01-6.408e-01j, -5.628e+00-1.028e+00j, 3.637e+00+5.146e-01j, -3.498e+00-6.270e-01j],
      [ 9.918e-01-1.579e-02j, 3.545e+00-3.832e+00j, -3.405e+00+2.284e+00j, -5.881e+00-6.028e-01j],
      [-3.243e-02+1.393e-01j, 9.077e-01+1.124e+00j, -1.434e+00-1.622e+00j, -6.827e-01-1.797e+00j],
      [ 3.645e-03+2.359e-01j, 1.521e+00+7.309e-01j, 1.015e-01-9.865e-01j, 2.600e+00-3.367e+00j]]]
      )
      return m, k, q
      def main():
      # Set up non-iterative p-k method
      nipk = init_pypk(get_cfg())
      m, k, q = get_matrices()
      nipk.set_matrices(m, k, q)
      # Compute solution
      nipk.compute()
      nipk.find_flutter()
      # Compute analytical gradients
      print('Computing NIPK gradients...')
      nipk.compute_gradients()
      np.savetxt('gradients_agard.dat', nipk.damp_k, fmt='%+.4e', delimiter=', ', header='dg/dk_i')
      # Compute gradients using finite differences
      eps = 1e-3
      dg = np.zeros(k.size)
      for i in range(k.size):
      idx = np.unravel_index(i, k.shape)
      k_ref = k[idx]
      k[idx] = k_ref + eps
      nipk.compute()
      gp = nipk.damp_agg
      k[idx] = k_ref - eps
      nipk.compute()
      gm = nipk.damp_agg
      dg[i] = (gp - gm) / (2 * eps)
      k[idx] = k_ref
      # Check
      for i in range(k.size):
      print(f'- checking gradient {i}/{k.size}')
      diff = abs(nipk.damp_k[i] - dg[i]) / nipk.damp_k[i]
      if diff > 1e-3:
      raise RuntimeError('Test failed: relative error {} between analtical derivative ({}) and finite difference derivative ({}) is too large!\n'.format(diff, nipk.damp_k[i], dg[i]))
      print('All tests ok.')
      if __name__ == '__main__':
      main()
      0% or .
      You are about to add 0 people to the discussion. Proceed with caution.
      Finish editing this message first!
      Please register or to comment