Enable partials. Add save. Add scenario name.

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

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)

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)

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
 

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

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

tests/agard_grad.py

+#!/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()