diff --git a/pypk/api_core.py b/pypk/api_core.py
index ff57890cdc3241c909abf76a185bc067e483fde2..d5bab2b2c88f8426f6b62935965d1bf79e609c3a 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
@@ -38,9 +40,9 @@ def init_pypk(cfg):
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)
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/nipk.py b/pypk/nipk.py
index 7b5110633e9785d449b3c573e45fe4964196d6b1..c6e875996d848ecdb42ced55046f3c7ce29dd3b5 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)
diff --git a/pypk/pk.py b/pypk/pk.py
index e4aa9d28d6278d108c312ff664068ba4d936aa26..695ba7f4c12401625407f40fe1073e9f4479fbd7 100644
--- a/pypk/pk.py
+++ b/pypk/pk.py
@@ -22,8 +22,8 @@ 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):
+ 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
diff --git a/pypk/solution.py b/pypk/solution.py
index a14acc25f015ade6a719887d9db2833a2650023e..d2c9c7cf1d127750b4acce63ecaf6ab9dc43fa86 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
@@ -69,7 +70,7 @@ class Solution:
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,21 @@ 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)
+ np.savetxt(f'results_{self._name}{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')
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 +117,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/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..e44e5c495abab311e0ba21c86532839b052275f3
--- /dev/null
+++ b/tests/agard_grad.py
@@ -0,0 +1,94 @@
+#!/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
+ '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()