Merge pull request #124 from ImperialCollegeLondon/dev_cfl_linearuvlm

Update linear UVLM to account for CFL not equal to one in the wake convection

sharpy/linear/assembler/linearaeroelastic.py

@@ -1,11 +1,13 @@
-import sharpy.linear.utils.ss_interface as ss_interface
 import numpy as np
-import sharpy.linear.src.libss as libss
 import scipy.linalg as sclalg
 import warnings
+
+import sharpy.linear.utils.ss_interface as ss_interface
+import sharpy.linear.src.libss as libss
 import sharpy.utils.settings as settings
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.algebra as algebra
+import sharpy.utils.generator_interface as gi
 
 
 @ss_interface.linear_system
@@ -98,8 +100,44 @@ def initialise(self, data):
         # Create Linear UVLM
         self.uvlm = ss_interface.initialise_system('LinearUVLM')
         self.uvlm.initialise(data, custom_settings=self.settings['aero_settings'])
+
+        # Look for the aerodynamic solver
+        if 'StaticUvlm' in data.settings:
+            aero_solver_name = 'StaticUvlm'
+            aero_solver_settings = data.settings['StaticUvlm']
+        elif 'StaticCoupled' in data.settings:
+            aero_solver_name = data.settings['StaticCoupled']['aero_solver']
+            aero_solver_settings = data.settings['StaticCoupled']['aero_solver_settings']
+        elif 'StaticCoupledRBM' in data.settings:
+            aero_solver_name = data.settings['StaticCoupledRBM']['aero_solver']
+            aero_solver_settings = data.settings['StaticCoupledRBM']['aero_solver_settings']
+        elif 'DynamicCoupled' in data.settings:
+            aero_solver_name = data.settings['DynamicCoupled']['aero_solver']
+            aero_solver_settings = data.settings['DynamicCoupled']['aero_solver_settings']
+        elif 'StepUvlm' in data.settings:
+            aero_solver_name = 'StepUvlm'
+            aero_solver_settings = data.settings['StepUvlm']
+        else:
+            raise RuntimeError("ERROR: aerodynamic solver not found")
+
+        # Velocity generator
+        vel_gen_name = aero_solver_settings['velocity_field_generator']
+        vel_gen_settings = aero_solver_settings['velocity_field_input']
+
+        # Get the minimum parameters needed to define the wake
+        vel_gen_type = gi.generator_from_string(vel_gen_name)
+        vel_gen = vel_gen_type()
+        vel_gen.initialise(vel_gen_settings) 
+
+        wake_prop_settings = {'dt': self.settings['aero_settings']['dt'],
+                              'ts': data.ts,
+                              't': data.ts*self.settings['aero_settings']['dt'],
+                              'for_pos': data.structure.timestep_info[-1].for_pos,
+                              'cfl1': self.settings['aero_settings']['cfl1'],
+                              'vel_gen': vel_gen}
+
         if self.settings['uvlm_filename'] == '':
-            self.uvlm.assemble(track_body=self.settings['track_body'])
+            self.uvlm.assemble(track_body=self.settings['track_body'], wake_prop_settings=wake_prop_settings)
         else:
             self.load_uvlm_from_file = True
 

sharpy/linear/assembler/linearuvlm.py

@@ -128,6 +128,10 @@ class LinearUVLM(ss_interface.BaseElement):
     settings_default['vortex_radius'] = vortex_radius_def
     settings_description['vortex_radius'] = 'Distance below which inductions are not computed'
 
+    settings_types['cfl1'] = 'bool'
+    settings_default['cfl1'] = True
+    settings_description['cfl1'] = 'If it is ``True``, it assumes that the discretisation complies with CFL=1'
+
     settings_table = settings.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
@@ -238,7 +242,7 @@ def initialise(self, data, custom_settings=None):
             self.gust_assembler = lineargust.LinearGustGenerator()
             self.gust_assembler.initialise(data.aero)
 
-    def assemble(self, track_body=False):
+    def assemble(self, track_body=False, wake_prop_settings=None):
         r"""
         Assembles the linearised UVLM system, removes the desired inputs and adds linearised control surfaces
         (if present).
@@ -252,7 +256,7 @@ def assemble(self, track_body=False):
         .. math:: [\delta_1, \delta_2, \dots, \dot{\delta}_1, \dot{\delta_2}]
         """
 
-        self.sys.assemble_ss()
+        self.sys.assemble_ss(wake_prop_settings=wake_prop_settings)
 
         if self.scaled:
             self.sys.nondimss()

sharpy/linear/src/assembly.py

@@ -24,6 +24,7 @@
 import sharpy.linear.src.lib_dbiot as dbiot
 import sharpy.linear.src.lib_ucdncdzeta as lib_ucdncdzeta
 import sharpy.utils.algebra as algebra
+import sharpy.utils.cout_utils as cout
 
 # local indiced panel/vertices as per self.maps
 dmver = [0, 1, 1, 0]  # delta to go from (m,n) panel to (m,n) vertices
@@ -1155,26 +1156,42 @@ def dfunstdgamma_dot(Surfs):
     return DerList
 
 
-def wake_prop(Surfs, Surfs_star, use_sparse=False, sparse_format='lil'):
+def wake_prop(MS, use_sparse=False, sparse_format='lil', settings=None):
     """
     Assembly of wake propagation matrices, in sparse or dense matrices format
 
-    Note: wake propagation matrices are very sparse. Nonetheless, allocation
-    in dense format (from numpy.zeros) or sparse does not have important
-    differences in terms of cpu time and memory used as numpy.zeros does
-    not allocate memory until this is accessed.
+    Note:
+        Wake propagation matrices are very sparse. Nonetheless, allocation
+        in dense format (from numpy.zeros) or sparse does not have important
+        differences in terms of cpu time and memory used as numpy.zeros does
+        not allocate memory until this is accessed
+
+    Args:
+        MS (MultiSurface): MultiSurface instance
+        use_sparse (bool (optional)): Use sparse matrices
+        sparse_format (str (optional)): Use either ``csc`` or ``lil`` format
+        settings (dict (optional)): Dictionary with aerodynamic settings containing:
+            cfl1 (bool): Defines if the wake shape complies with CFL=1
+            dt (float): time step
     """
 
-    n_surf = len(Surfs)
-    assert len(Surfs_star) == n_surf, 'No. of wake and bound surfaces not matching!'
+    try:
+        cfl1 = settings['cfl1']
+    except (KeyError, TypeError):
+        # In case the key does not exist or settings=None
+        cfl1 = True
+    cout.cout_wrap("Computing wake propagation matrix with CFL1={}".format(cfl1), 1)
+
+    n_surf = len(MS.Surfs)
+    assert len(MS.Surfs_star) == n_surf, 'No. of wake and bound surfaces not matching!'
 
     dimensions = [None]*n_surf
     dimensions_star = [None]*n_surf
 
     for ss in range(n_surf):
 
-        Surf = Surfs[ss]
-        Surf_star = Surfs_star[ss]
+        Surf = MS.Surfs[ss]
+        Surf_star = MS.Surfs_star[ss]
 
         N, M, K = Surf.maps.N, Surf.maps.M, Surf.maps.K
         M_star, K_star = Surf_star.maps.M, Surf_star.maps.K
@@ -1184,10 +1201,64 @@ def wake_prop(Surfs, Surfs_star, use_sparse=False, sparse_format='lil'):
         dimensions[ss] = [M, N, K]
         dimensions_star[ss] = [M_star, N, K_star]
 
-    C_list, Cstar_list = wake_prop_from_dimensions(dimensions,
-                                                   dimensions_star,
-                                                   use_sparse=use_sparse,
-                                                   sparse_format=sparse_format)
+        if not cfl1:
+            # allocate...
+            if use_sparse:
+                if sparse_format == 'csc':
+                    C = libsp.csc_matrix((K_star, K))
+                    C_star = libsp.csc_matrix((K_star, K_star))
+                elif sparse_format == 'lil':
+                    C = sparse.lil_matrix((K_star, K))
+                    C_star = sparse.lil_matrix((K_star, K_star))
+            else:
+                C = np.zeros((K_star, K))
+                C_star = np.zeros((K_star, K_star))
+
+            C_list = []
+            Cstar_list = []
+
+            # Compute flow velocity at wake
+            uext = [np.zeros((3,
+                             dimensions_star[ss][0],
+                             dimensions_star[ss][1]))]
+
+            try:
+                Surf_star.zetac
+            except AttributeError:
+                Surf_star.generate_collocations()
+            # Compute induced velocities in the wake
+            Surf_star.u_ind_coll = np.zeros((3, M_star, N))
+
+            for iin in range(N):
+                # propagation from trailing edge
+                conv_vec = Surf_star.zetac[:, 0, iin] - Surf.zetac[:, -1, iin]
+                dist = np.linalg.norm(conv_vec)
+                conv_dir_te = conv_vec/dist
+                vel = Surf.u_input_coll[:, -1, iin]
+                vel_value = np.dot(vel, conv_dir_te)
+                cfl = settings['dt']*vel_value/dist
+
+                C[iin, N * (M - 1) + iin] = cfl
+                C_star[iin, iin] = 1.0 - cfl
+
+                # wake propagation
+                for mm in range(1, M_star):
+                    conv_vec = Surf_star.zetac[:, mm, iin] - Surf_star.zetac[:, mm - 1, iin]
+                    dist = np.linalg.norm(conv_vec)
+                    conv_dir = conv_vec/dist
+                    cfl = settings['dt']*vel_value/dist
+
+                    C_star[mm * N + iin, (mm - 1) * N + iin] = cfl
+                    C_star[mm * N + iin, mm * N + iin] = 1.0 - cfl
+
+            C_list.append(C)
+            Cstar_list.append(C_star)
+
+    if cfl1:
+        C_list, Cstar_list = wake_prop_from_dimensions(dimensions,
+                                                       dimensions_star,
+                                                       use_sparse=use_sparse,
+                                                       sparse_format=sparse_format)
 
     return C_list, Cstar_list
 

sharpy/linear/src/lin_aeroelastic.py

@@ -203,7 +203,7 @@ def reshape_struct_input(self):
 
     # self.dq[-3:]=-0.5*(wa*tsdata.quat[0]+np.cross(wa,tsdata.quat[1:]))
 
-    def assemble_ss(self, beam_num_modes=None):
+    def assemble_ss(self, beam_num_modes=None, wake_prop_settings=None):
         """Assemble State Space formulation"""
         data = self.data
 
@@ -213,7 +213,7 @@ def assemble_ss(self, beam_num_modes=None):
         tsstr = self.tsstr
 
         ### assemble linear uvlm
-        self.linuvlm.assemble_ss()
+        self.linuvlm.assemble_ss(wake_prop_settings=wake_prop_settings)
         SSaero = self.linuvlm.SS
 
         ### assemble gains and stiffening term due to non-zero forces