diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2f4f3ca4b36dcd3b2e7496634eaacbc1831681ac..a9fc51d82de9710ac7fe9abcec88b9a586a0c044 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,5 +1,5 @@
default:
- image: rboman/waves-py3:2020.3
+ image: rboman/waves-py3:2022.0
.global_tag: &global_tag_def
tags:
diff --git a/pypk/__init__.py b/pypk/__init__.py
index 873ccbb461dcad5bc403516eb1b447b9251aee7c..38e750a84e338ef497d4da65167bd2a5b0130c25 100644
--- a/pypk/__init__.py
+++ b/pypk/__init__.py
@@ -15,6 +15,6 @@
# See the License for the specific language governing permissions and
# limitations under the License.
-__version__ = '1.0.0'
+__version__ = '1.1.0'
from .api_core import init_pypk
diff --git a/pypk/api_core.py b/pypk/api_core.py
index ff57890cdc3241c909abf76a185bc067e483fde2..3d61916cd5e28ae8d1dee3d5acc33d148277c428 100644
--- a/pypk/api_core.py
+++ b/pypk/api_core.py
@@ -23,6 +23,8 @@ from .nipk import Nipk
def init_pypk(cfg):
"""Initialize PyPK
"""
+ # Case name
+ name = cfg['name'] if 'name' in cfg else 'noname'
# Fluid type
if cfg['fluid'] == 'unmatched':
n_speed = len(cfg['u_idx']) # number of airspeed test points
@@ -33,14 +35,15 @@ def init_pypk(cfg):
else:
raise RuntimeError('"fluid" entry must be "unmatched" or "matched_isa"!\n')
# 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
agg = KSAggregator(cfg['rho_ks'])
# p-k method
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':
- 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:
raise RuntimeError('"method" entry must be "pk" or "nipk"!\n')
return sol
diff --git a/pypk/api_om.py b/pypk/api_om.py
index 95861575dcc4c16c33f010ffc807eb98d45b3b9e..894bfe0a2968c69cdec27e6fed293aa73363bc28 100644
--- a/pypk/api_om.py
+++ b/pypk/api_om.py
@@ -23,11 +23,13 @@ class PypkSolver(om.ExplicitComponent):
"""OpenMDAO flutter solution wrapper
"""
def initialize(self):
+ self.options.declare('scn_name', desc='name of current scenario')
self.options.declare('nm', desc='number of modes')
self.options.declare('nu', desc='number of airspeed test points')
self.options.declare('sol', desc='PyPK solver', recordable=False)
def setup(self):
+ self.scn_name = self.options['scn_name']
self.nm = self.options['nm']
self.nu = self.options['nu']
self.sol = self.options['sol']
@@ -41,18 +43,26 @@ class PypkSolver(om.ExplicitComponent):
self.add_output('f_speed', val=1., desc='flutter speed')
self.add_output('f_mode', val=1, desc='flutter mode')
# 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):
self.sol.set_matrices(inputs['M'], inputs['K'], inputs['Q_re'] + 1j * inputs['Q_im'])
self.sol.compute()
self.sol.find_flutter()
+ self.sol.save('_' + self.scn_name)
outputs['f'] = self.sol.freq
outputs['g'] = self.sol.damp
outputs['g_agg'] = self.sol.damp_agg
outputs['f_speed'] = self.sol.f_speed
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):
"""PyPK builder for OpenMDAO
"""
@@ -66,7 +76,7 @@ class PypkBuilder(Builder):
"""
self.__sol = init_pypk(cfg)
- def get_solver(self):
+ def get_solver(self, scenario_name=''):
"""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)
diff --git a/pypk/eigensolver.py b/pypk/eigensolver.py
index d0fe1eb8f684f27461cf142f9fa3971ef1e0781b..6c2872158884f39021a5a432a937c2f5e0bbfc51 100644
--- a/pypk/eigensolver.py
+++ b/pypk/eigensolver.py
@@ -20,26 +20,36 @@ import scipy.linalg as spla
class EigenSolver:
"""Eigen solver
"""
- def __init__(self, l, n):
+ def __init__(self, l, n, g=0):
self._l = l # reference length
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(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)
# 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):
@@ -58,10 +68,23 @@ class EigenSolver:
a = np.zeros((self._n + 1, self._n + 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,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,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
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 7b5110633e9785d449b3c573e45fe4964196d6b1..c097f7407a7afba49341c14483010ddac45c186a 100644
--- a/pypk/nipk.py
+++ b/pypk/nipk.py
@@ -20,8 +20,8 @@ import numpy as np
class Nipk(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):
- super().__init__(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__(name, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed)
self.name = 'NIPK' # solver name
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)
@@ -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 e4aa9d28d6278d108c312ff664068ba4d936aa26..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, flu, eig, agg, mach, k_ref, l_ref, vrb, n_modes, n_speed, tol=1e-3):
- super().__init__(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, 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
diff --git a/pypk/solution.py b/pypk/solution.py
index a14acc25f015ade6a719887d9db2833a2650023e..90805758716e4a13b875ac6673b3d1f2bae597c7 100644
--- a/pypk/solution.py
+++ b/pypk/solution.py
@@ -19,8 +19,9 @@ import numpy as np
class Solution:
"""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 = name # case name
self._flu = flu # fluid type
self._eig = eig # eigen solver
self._agg = agg # aggregator
@@ -61,15 +62,15 @@ class Solution:
"""Find the flutter speed and display it
"""
found = False
- self.f_speed = np.inf
- self.f_mode = np.inf
+ self.f_speed = -1.
+ self.f_mode = -1
for ivel in range(len(self.speed) - 1):
if (found):
break
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._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_mode = imod
found = True
@@ -77,21 +78,20 @@ class Solution:
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]))
- def save(self, sfx):
+ def save(self, sfx=''):
"""Write results to disk
sfx : filename suffix
"""
- hdr = '{:>8s}, {:>10s}, {:>10s}, {:>10s}'.format('Mach', 'u_f', 'mode', 'AGG(g)')
- 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 = 'Mach={:.2f}, u_f={:.2f}, mode={:d}, AGG(g)={:.2f}\n'.format(self._mach, self.f_speed, self.f_mode, self.damp_agg)
+ hdr += '{:>9s}'.format('u_inf')
for j in range(self.n_modes):
hdr += ', {:>11s}'.format(f'f_{j}')
for j in range(self.n_modes):
hdr += ', {:>11s}'.format(f'g_{j}')
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
show : whether to display the plot
format : format to which the plot will be written to disk, if not empty
@@ -116,5 +116,5 @@ class Solution:
axs[1].spines['right'].set_visible(False)
axs[1].plot(u, [0.] * self.n_speed, color='k', lw=1, ls='-.') # zero-damping line
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()
diff --git a/setup.py b/setup.py
index f04087f82d41674f2bc65cd17ed37632b3995722..86ec6a775c8695c934be915e5b0507e98c4108cd 100644
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,6 @@ setup(
url='https://gitlab.uliege.be/am-dept/pypk',
license='Apache 2.0',
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'],
)
diff --git a/tests/agard.py b/tests/agard.py
index 489f3395e20adf2d7fa1715c1ccef55be4ba08ba..a62f30fb5bd5029657c5694eab46668054eee5da 100644
--- a/tests/agard.py
+++ b/tests/agard.py
@@ -23,6 +23,7 @@ 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
@@ -60,8 +61,8 @@ def compute(sol):
sol.compute()
sol.find_flutter()
# Save results
- sol.save('agard_' + sol.name.lower())
- sol.plot('agard_' + sol.name.lower(), show=False, format='.pdf')
+ sol.save('_' + sol.name.lower())
+ sol.plot('_' + sol.name.lower(), show=False, format='pdf')
def main():
# Set up p-k and non-iterative p-k method
diff --git a/tests/agard_grad.py b/tests/agard_grad.py
new file mode 100644
index 0000000000000000000000000000000000000000..77ddfdfdac2f83f626d19f3602b84849dd8f4416
--- /dev/null
+++ b/tests/agard_grad.py
@@ -0,0 +1,95 @@
+#!/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()