Merge pull request #23 from ImperialCollegeLondon/dev_linuvlm

Improved linear solver documentation and minor Krylov ROM fixes

sharpy/linear/assembler/lincontrolsurfacedeflector.py

@@ -61,7 +61,7 @@ def assemble(self):
         """
         pass
 
-    def generate(self, linuvlm=None, tsaero0=None, tsstruct0=None, aero=None, structure=None, track_body=False):
+    def generate(self, linuvlm=None, tsaero0=None, tsstruct0=None, aero=None, structure=None):
         """
         Generates a matrix mapping a linear control surface deflection onto the aerodynamic grid.
 
@@ -97,7 +97,6 @@ def generate(self, linuvlm=None, tsaero0=None, tsstruct0=None, aero=None, struct
         Kdisp = np.zeros((3 * linuvlm.Kzeta, n_control_surfaces))
         Kvel = np.zeros((3 * linuvlm.Kzeta, n_control_surfaces))
         Kmom = np.zeros((3 * linuvlm.Kzeta, n_control_surfaces))
-        Knew = np.zeros((3 * linuvlm.Kzeta, 3 * linuvlm.Kzeta + n_control_surfaces))
         zeta0 = np.concatenate([tsaero0.zeta[i_surf].reshape(-1, order='C') for i_surf in range(n_surf)])
 
         Cga = algebra.quat2rotation(tsstruct0.quat).T

sharpy/linear/assembler/linearaeroelastic.py

@@ -1,5 +1,4 @@
 import sharpy.linear.utils.ss_interface as ss_interface
-import scipy.sparse as scsp
 import numpy as np
 import sharpy.linear.src.libss as libss
 import scipy.linalg as sclalg
@@ -91,7 +90,7 @@ def initialise(self, data):
             self.settings = data.settings['LinearAssembler']['linear_system_settings']
         except KeyError:
             self.settings = None
-        settings.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
 
         if self.settings['use_euler']:
             self.settings['aero_settings']['use_euler'] = True
@@ -101,7 +100,7 @@ def initialise(self, data):
         self.uvlm = ss_interface.initialise_system('LinearUVLM')
         self.uvlm.initialise(data, custom_settings=self.settings['aero_settings'])
         if self.settings['uvlm_filename'] == '':
-            self.uvlm.assemble(track_body=self.settings['track_body'].value)
+            self.uvlm.assemble(track_body=self.settings['track_body'])
         else:
             self.load_uvlm_from_file = True
 
@@ -197,7 +196,7 @@ def assemble(self):
             self.couplings['Ksa'] = Ksa
             self.couplings['Kas'] = Kas
 
-            if self.settings['beam_settings']['modal_projection'].value == True and \
+            if self.settings['beam_settings']['modal_projection'] == True and \
                     self.settings['beam_settings']['inout_coords'] == 'modes':
                 # Project UVLM onto modal space
                 phi = beam.sys.U
@@ -232,7 +231,7 @@ def assemble(self):
                 Tas[:flex_nodes + 6, :flex_nodes + 6] /= uvlm.sys.ScalingFacts['length']
                 Tas[total_dof: total_dof + flex_nodes + 6] /= uvlm.sys.ScalingFacts['length']
             else:
-                if not self.settings['beam_settings']['modal_projection'].value:
+                if not self.settings['beam_settings']['modal_projection']:
                     Tas /= uvlm.sys.ScalingFacts['length']
 
         ss = libss.couple(ss01=uvlm.ss, ss02=beam.ss, K12=Tas, K21=Tsa)
@@ -255,7 +254,7 @@ def assemble(self):
         # Save zero force reference
         self.linearisation_vectors['forces_aero_beam_dof'] = Ksa.dot(self.linearisation_vectors['forces_aero'])
 
-        if self.settings['beam_settings']['modal_projection'].value == True and \
+        if self.settings['beam_settings']['modal_projection'] == True and \
                 self.settings['beam_settings']['inout_coords'] == 'modes':
             self.linearisation_vectors['forces_aero_beam_dof'] = out_mode_matrix.dot(self.linearisation_vectors['forces_aero_beam_dof'])
 
@@ -299,7 +298,7 @@ def runrom_rbm(self, uvlm):
         else:
             orient_dof = 4
         # Input side
-        if self.settings['beam_settings']['modal_projection'].value == True and \
+        if self.settings['beam_settings']['modal_projection'] == True and \
                 self.settings['beam_settings']['inout_coords'] == 'modes':
             rem_int_modes = np.zeros((ss.inputs, ss.inputs - rig_nodes))
             rem_int_modes[rig_nodes:, :] = np.eye(ss.inputs - rig_nodes)

sharpy/linear/assembler/linearbeam.py

@@ -108,7 +108,7 @@ class LinearBeam(BaseElement):
     settings_types['gravity'] = 'bool'
     settings_description['gravity'] = 'Linearise gravitational forces'
 
-    settings_types['remove_dofs'] = 'list'
+    settings_types['remove_dofs'] = 'list(str)'
     settings_default['remove_dofs'] = []
     settings_description['remove_dofs'] = 'Remove desired degrees of freedom: ``eta``, ``V``, ``W`` or ``orient``'
 
@@ -174,13 +174,13 @@ def assemble(self, t_ref=None):
         Returns:
 
         """
-        if self.settings['gravity'].value:
+        if self.settings['gravity']:
             self.sys.linearise_gravity_forces()
 
         if self.settings['remove_dofs']:
             self.trim_nodes(self.settings['remove_dofs'])
 
-        if self.settings['modal_projection'].value and self.settings['remove_sym_modes'].value and self.clamped:
+        if self.settings['modal_projection'] and self.settings['remove_sym_modes'] and self.clamped:
             self.remove_symmetric_modes()
 
         if t_ref is not None:

sharpy/linear/assembler/linearuvlm.py

@@ -13,17 +13,71 @@
 
 @ss_interface.linear_system
 class LinearUVLM(ss_interface.BaseElement):
-    """
+    r"""
     Linear UVLM System Assembler
 
+    Produces state-space model of the form
+
+        .. math::
+
+            \mathbf{x}_{n+1} &= \mathbf{A}\,\mathbf{x}_n + \mathbf{B} \mathbf{u}_{n+1} \\
+            \mathbf{y}_n &= \mathbf{C}\,\mathbf{x}_n + \mathbf{D} \mathbf{u}_n
+
+    where the state, inputs and outputs are:
+
+        .. math:: \mathbf{x}_n = \{ \delta \mathbf{\Gamma}_n,\, \delta \mathbf{\Gamma_{w_n}},\,
+            \Delta t\,\delta\mathbf{\Gamma}'_n,\, \delta\mathbf{\Gamma}_{n-1} \}
+
+        .. math:: \mathbf{u}_n = \{ \delta\mathbf{\zeta}_n,\, \delta\mathbf{\zeta}'_n,\,
+            \delta\mathbf{u}_{ext,n} \}
+
+        .. math:: \mathbf{y} = \{\delta\mathbf{f}\}
+
+    with :math:`\mathbf{\Gamma}\in\mathbb{R}^{MN}` being the vector of vortex circulations,
+    :math:`\mathbf{\zeta}\in\mathbb{R}^{3(M+1)(N+1)}` the vector of vortex lattice coordinates and
+    :math:`\mathbf{f}\in\mathbb{R}^{3(M+1)(N+1)}` the vector of aerodynamic forces and moments. Note
+    that :math:`(\bullet)'` denotes a derivative with respect to time.
+
+    Note that the input is atypically defined at time ``n+1``. If the setting
+    ``remove_predictor = True`` the predictor term ``u_{n+1}`` is eliminated through
+    the change of state[1]:
+
+        .. math::
+            \mathbf{h}_n &= \mathbf{x}_n - \mathbf{B}\,\mathbf{u}_n \\
+
+    such that:
+
+        .. math::
+            \mathbf{h}_{n+1} &= \mathbf{A}\,\mathbf{h}_n + \mathbf{A\,B}\,\mathbf{u}_n \\
+            \mathbf{y}_n &= \mathbf{C\,h}_n + (\mathbf{C\,B}+\mathbf{D})\,\mathbf{u}_n
+
+    which only modifies the equivalent :math:`\mathbf{B}` and :math:`\mathbf{D}` matrices.
+
+    The ``integr_order`` setting refers to the finite differencing scheme used to calculate the bound circulation
+    derivative with respect to time :math:`\dot{\mathbf{\Gamma}}`. A first order scheme is used when
+    ``integr_order == 1``
+
+    .. math:: \dot{\mathbf{\Gamma}}^{n+1} = \frac{\mathbf{\Gamma}^{n+1}-\mathbf{\Gamma}^n}{\Delta t}
+
+    If ``integr_order == 2`` a higher order scheme is used (but it isn't exactly second order accurate [1]).
+
+    .. math:: \dot{\mathbf{\Gamma}}^{n+1} = \frac{3\mathbf{\Gamma}^{n+1}-4\mathbf{\Gamma}^n + \mathbf{\Gamma}^{n-1}}
+        {2\Delta t}
+
+    References:
+        [1] Franklin, GF and Powell, JD. Digital Control of Dynamic Systems, Addison-Wesley Publishing Company, 1980
+
+        [2] Maraniello, S., & Palacios, R.. State-Space Realizations and Internal Balancing in Potential-Flow
+            Aerodynamics with
+            Arbitrary Kinematics. AIAA Journal, 57(6), 1–14. 2019. https://doi.org/10.2514/1.J058153
 
     """
     sys_id = 'LinearUVLM'
 
     settings_types = dict()
     settings_default = dict()
     settings_description = dict()
-    
+
     settings_types['dt'] = 'float'
     settings_default['dt'] = 0.1
     settings_description['dt'] = 'Time step'
@@ -48,7 +102,7 @@ class LinearUVLM(ss_interface.BaseElement):
     settings_default['density'] = 1.225
     settings_description['density'] = 'Air density'
 
-    settings_types['remove_inputs'] = 'list'
+    settings_types['remove_inputs'] = 'list(str)'
     settings_default['remove_inputs'] = []
     settings_description['remove_inputs'] = 'List of inputs to remove. ``u_gust`` to remove external velocity input.'
 
@@ -89,7 +143,7 @@ class LinearUVLM(ss_interface.BaseElement):
                                        scaling_settings_description)
 
     def __init__(self):
-        
+
         self.sys = None
         self.ss = None
         self.tsaero0 = None

sharpy/linear/src/lingebm.py

@@ -241,12 +241,6 @@ def num_modes(self):
 
     @num_modes.setter
     def num_modes(self, value):
-        # self.U = self.U[:, :value]
-        # self.freq_natural = self.freq_natural[:value]
-        # if self.freq_damp is not None:
-        #     self.freq_damp = self.freq_damp[:value]
-        # if self.V is not None:
-        #     self.V = self.V[:value, :]
         self.update_truncated_modes(value)
         self._num_modes = value
 
@@ -501,7 +495,7 @@ def linearise_gravity_forces(self, tsstr=None):
         if tsstr is None:
             tsstr = self.tsstruct0
 
-        if self.settings['print_info'].value:
+        if self.settings['print_info']:
             try:
                 cout.cout_wrap('\nLinearising gravity terms...')
             except ValueError:
@@ -543,7 +537,7 @@ def linearise_gravity_forces(self, tsstr=None):
         Xcg_A = -np.array([Mrr[2, 4], Mrr[0, 5], Mrr[1, 3]]) / Mrr[0, 0]
         Xcg_Askew = algebra.skew(Xcg_A)
 
-        if self.settings['print_info'].value:
+        if self.settings['print_info']:
             cout.cout_wrap('\tM = %.2f kg' % Mrr[0, 0], 1)
             cout.cout_wrap('\tX_CG A -> %.2f %.2f %.2f' %(Xcg_A[0], Xcg_A[1], Xcg_A[2]), 1)
 
@@ -604,7 +598,7 @@ def linearise_gravity_forces(self, tsstr=None):
                 # Nodal CG in G frame - debug
                 Xcg_G_n = Pga.dot(Xcg_A_n)
 
-                if self.settings['print_info'].value:
+                if self.settings['print_info']:
                     cout.cout_wrap("Node %2d \t-> B %.3f %.3f %.3f" %(i_node, Xcg_B[0], Xcg_B[1], Xcg_B[2]), 2)
                     cout.cout_wrap("\t\t\t-> A %.3f %.3f %.3f" %(Xcg_A_n[0], Xcg_A_n[1], Xcg_A_n[2]), 2)
                     cout.cout_wrap("\t\t\t-> G %.3f %.3f %.3f" %(Xcg_G_n[0], Xcg_G_n[1], Xcg_G_n[2]), 2)
@@ -686,7 +680,7 @@ def linearise_gravity_forces(self, tsstr=None):
         if self.modal:
             self.Ccut = self.U.T.dot(self.Cstr.dot(self.U))
 
-        if self.settings['print_info'].value:
+        if self.settings['print_info']:
             cout.cout_wrap('\tUpdated the beam C, modal C and K matrices with the terms from the gravity linearisation\n')
 
     def assemble(self, Nmodes=None):
@@ -723,12 +717,12 @@ def assemble(self, Nmodes=None):
         assert self.inout_coords in ['modes', 'nodes'], \
             'inout_coords=%s not implemented!' % self.inout_coords
 
-        cond_mass_matrix = np.linalg.cond(self.Mstr)
-        if np.log10(cond_mass_matrix) >= 10.:
-            warnings.warn('Mass matrix is poorly conditioned (Cond = 10^%f). Inverse may not be correct.'
-                          % np.log10(cond_mass_matrix), 3)
-        else:
-            cout.cout_wrap('Mass matrix condition = %e' % cond_mass_matrix)
+        # cond_mass_matrix = np.linalg.cond(self.Mstr)
+        # if np.log10(cond_mass_matrix) >= 10.:
+        #     warnings.warn('Mass matrix is poorly conditioned (Cond = 10^%f). Inverse may not be correct.'
+        #                   % np.log10(cond_mass_matrix), 3)
+        # else:
+        #     cout.cout_wrap('Mass matrix condition = %e' % cond_mass_matrix)
 
         dlti = self.dlti
         modal = self.modal