add regression test for gust vane generator

Introduce a 6‑timestep gust vane simulation test to ensure backward
compatibility. The test compares the simulation output against the
current reference results to detect unintended changes in gust vane
behavior.
Adapt pazy related test modules to integrate gust vane generator functionalities.

Fix missing import in writevariables time.

Fix airfoil polar code.

Author: sduess

sharpy/aero/utils/airfoilpolars.py

@@ -61,7 +61,7 @@ def get_coefs(self, aoa_deg):
         cd = self.cd_interp(aoa_deg)
         cm = self.cm_interp(aoa_deg)
 
-        return cl[0], cd[0], cm[0]
+        return cl, cd, cm
 
     def get_aoa_deg_from_cl_2pi(self, cl):
 
@@ -115,7 +115,7 @@ def get_cdcm_from_cl(self, cl):
                 cd = np.interp(cl, self.table[i:i+2, 1], self.table[i:i+2, 2])
                 cm = np.interp(cl, self.table[i:i+2, 1], self.table[i:i+2, 3])
 
-        return float(cd), float(cm)
+        return cd, cm
 
 
 def interpolate(polar1, polar2, coef=0.5):

sharpy/postproc/writevariablestime.py

@@ -2,7 +2,7 @@
 import numpy as np
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.settings as settings_utils
-
+import sharpy.aero.utils.uvlmlib as uvlmlib
 
 @solver
 class WriteVariablesTime(BaseSolver):

tests/coupled/dynamic/test_gust_vanes.py

@@ -0,0 +1,49 @@
+import unittest
+import numpy as np
+import os
+from tests.coupled.static.test_pazy_static import TestPazyCoupled
+
+class TestPazyCoupledDynamicWithGustVanes(TestPazyCoupled, unittest.TestCase):
+    """
+    Test gust vanes with pazy wing in a dynamic coupled simulation with a free wake convection scheme. The induced velocity field
+    and wing root loads and moments are compared after 20 timesteps with the reference case produced with SHARPy v2.0 for backward compability.
+    """
+
+    def test_dynamic_aoa_symmetry_with_gust_vanes(self):
+        self.setup_test_folders('gust_vanes')
+        cs_deflection_file = self.route_test_dir + '/gust_vanes/gust_vane_deflections.csv'
+         
+        self.num_chordwise_panels //= 2
+        self.num_spanwise_nodes //= 2   
+        self.n_tsteps = 6
+        
+        self.run_test(True, dynamic=True, gust_vanes=True, cs_deflection_file=cs_deflection_file)
+
+    def evaluate_output(self):  
+        self.assert_induced_velocity_matches_reference(
+            time_step=6,
+            ref_u_ind_z=-0.4026807
+        )
+
+    def assert_induced_velocity_matches_reference(self, time_step, ref_u_ind_z):
+        """
+        Check if the vertical induced velocity at a given time step and node matches the reference.
+
+        Args:
+            time_step (int): Index of the time step to check (e.g., 6)
+             ref_u_ind_z (float): Reference vertical induced velocity [m/s]
+        """
+        file_path = os.path.join(self.output_folder, self.case_name, 'WriteVariablesTime', 'vel_field_uind_point0.dat')
+        induced_velocities = np.loadtxt(file_path)
+
+        actual_value = induced_velocities[time_step, -1]
+
+        np.testing.assert_almost_equal(
+            actual_value, ref_u_ind_z, decimal=3,
+            err_msg='Induced vertical velocity differs more than 0.1% from reference value.',
+            verbose=True
+        )
+        
+if __name__ == '__main__':
+    unittest.main()
+

tests/coupled/dynamic/test_pazy_dynamic.py

@@ -3,36 +3,51 @@
 import os
 from tests.coupled.static.test_pazy_static import TestPazyCoupled
 
-class TestPazyCoupledDynamic(TestPazyCoupled):
+class TestPazyCoupledDynamic(TestPazyCoupled, unittest.TestCase):
     """
     Test Pazy wing dynamic coupled case with a free wake convection scheme by compparing wing root loads and 
-    moments after 20 timesteps with the reference case produced with SHARPy v2.0 for backward compability.
+    moments after 2 timesteps with the reference case produced with SHARPy v2.0 for backward compability.
     Further, symmetry condition is checked for a dynamic free wake as well.
     """
 
-    def test_dynamic_aoa(self):        
-        self.run_test(False, dynamic = True)
-
     def test_dynamic_aoa_symmetry(self):
+        self.setup_test_folders('pazy')
+        self.num_chordwise_panels //= 4
+        self.num_spanwise_nodes //= 4  
         self.run_test(True, dynamic = True)
 
-    def evaluate_output(self):       
-        ref_Fz = -1663.6376358639793
-        ref_My = -12.656888844293645
-
-        file = os.path.join(self.output_folder, self.case_name, 'beam/beam_loads_%i.csv' % (self.n_tsteps))
-        beam_loads_ts = np.loadtxt(file, delimiter=',')
-        error_Fz = (float(beam_loads_ts[0, 6])-ref_Fz)/ref_Fz
-        error_My = (float(beam_loads_ts[0, 8])-ref_My)/ref_My
-        
-        np.testing.assert_almost_equal(error_Fz, 0.,
-                                       decimal=3,
-                                       err_msg='Vertical load on wing root differs more than 0.1 %% from reference value.',
-                                       verbose=True)
-        np.testing.assert_almost_equal(error_My, 0.,
-                                       decimal=3,
-                                       err_msg='Pitching moment on wing root differs more than 0.1 %% from reference value.',
-                                       verbose=True)
+    def evaluate_output(self):                                                      
+        self.assert_root_forces_match_reference(ref_Fz=-2.5274941e+04, 
+                                                ref_My=-1.0703502e+01)
+
+    def assert_root_forces_match_reference(self, ref_Fz, ref_My):
+        """
+        Check if the vertical force and pitching moment at the wing root match reference values.
+
+        Args:
+            ref_Fz (float): Reference vertical force at wing root [N]
+            ref_My (float): Reference pitching moment at wing root [Nm]
+        """
+        file_path = os.path.join(self.output_folder, self.case_name, 'beam', f'beam_loads_{self.n_tsteps}.csv')
+        beam_loads_ts = np.loadtxt(file_path, delimiter=',')
+
+        actual_Fz = float(beam_loads_ts[0, 6])
+        actual_My = float(beam_loads_ts[0, 8])
+
+        error_Fz = (actual_Fz - ref_Fz) / ref_Fz
+        error_My = (actual_My - ref_My) / ref_My
+
+        np.testing.assert_almost_equal(
+            error_Fz, 0.0, decimal=3,
+            err_msg='Vertical load on wing root differs more than 0.1% from reference value.',
+            verbose=True
+        )
+        np.testing.assert_almost_equal(
+            error_My, 0.0, decimal=3,
+            err_msg='Pitching moment on wing root differs more than 0.1% from reference value.',
+            verbose=True
+        )    
+    
 
 if __name__ == '__main__':
     unittest.main()

tests/coupled/static/pazy/generate_pazy.py

@@ -9,9 +9,10 @@ def generate_pazy(u_inf, case_name, output_folder='/output/', cases_folder='', *
     num_modes = 16
     gravity_on = kwargs.get('gravity_on', True)
     symmetry_condition = kwargs.get('symmetry_condition', False)
+    gust_vanes = kwargs.get('gust_vanes', False)
     dynamic = kwargs.get('dynamic', False)
     n_tsteps = kwargs.get('n_tsteps', 1)
-
+    cs_deflection_file=kwargs.get('cs_deflection_file', None)
     # Lattice Discretisation
     M = kwargs.get('M', 4)
     N = kwargs.get('N', 32)
@@ -42,15 +43,19 @@ def generate_pazy(u_inf, case_name, output_folder='/output/', cases_folder='', *
         ws.reduce_model_to_symmetric_wing()
     ws.generate_aero_file()
     ws.generate_fem_file()
-    set_final_settings(ws, dynamic,
+    set_final_settings(ws, 
+                        dynamic=dynamic,
+                        surface_m=M,
                         output_folder=output_folder,
                         symmetry_condition = symmetry_condition,
+                        gust_vanes = gust_vanes,
                         gravity_on = gravity_on,
-                        n_tsteps=n_tsteps)
+                        n_tsteps=n_tsteps,
+                        cs_deflection_file=cs_deflection_file)
 
     sharpy.sharpy_main.main(['', ws.route + ws.case_name + '.sharpy'])
 
-def set_final_settings(ws, dynamic = False, output_folder='/output/', symmetry_condition = False, gravity_on = True, n_tsteps=1,flag_multiple_mstar_input=False):
+def set_final_settings(ws, dynamic = False, surface_m=8, output_folder='/output/', symmetry_condition = False, gravity_on = True, n_tsteps=1,flag_multiple_mstar_input=False, gust_vanes = False, cs_deflection_file=None):
     ws.config['SHARPy'] = {
         'flow':
             ['BeamLoader',
@@ -204,7 +209,80 @@ def set_final_settings(ws, dynamic = False, output_folder='/output/', symmetry_c
                                         },
             }
         ws.config['BeamLoads'] = {'csv_output': True}
-
+    if gust_vanes:
+        import numpy as np
+        ws.config = apply_gust_vane_settings(ws.config,
+                                             cs_deflection_file,
+                                             ws.dt,
+                                             ws.u_inf,
+                                             surface_m,
+                                             False,
+                                             symmetry_condition)
+        ws.config['DynamicCoupled']['postprocessors'].append('WriteVariablesTime')
+        ws.config['DynamicCoupled']['postprocessors_settings']['WriteVariablesTime'] = {
+            'vel_field_variables': ['uind'],
+            'vel_field_points': np.array([-1.25, 0.0, 0.25 ]),
+                                        }
     ws.config.write()
 
     sharpy.sharpy_main.main(['', ws.route + ws.case_name + '.sharpy'])
+
+def apply_gust_vane_settings(settings,
+                             cs_deflection_file,
+                             dt,
+                             u_inf,
+                             surface_m,
+                             vertical,
+                             symmetry_condition):
+    """
+    Updates the SHARPy configuration with gust vane definitions.
+
+    Args:
+        settings: SHARPy configuration object.
+        cs_deflection_file: Path to the gust vane deflection file.
+        dt: Time step size.
+        u_inf: Freestream velocity.
+        surface_m: Chordwise discretisation.
+        vertical: Orientation of the gust vanes (vertical or horizontal).
+        symmetry_condition: Whether the problem is symmetric.
+        only_gust_vanes: Whether gust vanes are the only gust source.
+
+    Returns:
+        cs_deflection_file. Modified the input settings dictionary.
+    """
+    # breakpoint()
+    wake_length_vanes = 5
+    gust_vane_parameters = {
+        'M': surface_m * 3,
+        'N': 40,
+        'M_star': int(wake_length_vanes / (dt * u_inf)),
+        'span': 5,
+        'chord': 0.3,
+        'control_surface_deflection_generator_settings': {
+            'dt': dt,
+            'deflection_file': cs_deflection_file
+        }
+    }
+
+    settings['AerogridLoader']['gust_vanes'] = True
+    settings['AerogridLoader']['gust_vanes_generator_settings'] = {
+        'n_vanes': 2,
+        'streamwise_position': [-1.5, -1.5],
+        'vertical_position': [-0.25, 0.25],
+        'symmetry_condition': symmetry_condition,
+        'vane_parameters': [gust_vane_parameters, gust_vane_parameters],
+        'vertical': vertical
+    }
+
+    # Override velocity field for gust vanes
+    stepuvlm_updates = {
+            'convection_scheme': 3,
+            'velocity_field_generator': 'SteadyVelocityField',
+            'velocity_field_input': {
+                'u_inf': u_inf,
+                'u_inf_direction': [1., 0., 0.],
+            }
+        }
+    settings['DynamicCoupled']['aero_solver_settings'].update(stepuvlm_updates)
+
+    return settings

tests/coupled/static/test_pazy_static.py

@@ -11,43 +11,46 @@ class TestPazyCoupled(unittest.TestCase):
     def setUp(self):
         self.u_inf = 50
         self.alpha = 7
-        self.M = 16
-        self.N= 64
+        self.num_chordwise_panels = 16
+        self.num_spanwise_nodes= 64
         self.Msf = 1
-
-        self.n_tsteps = 20
-
-        self.route_test_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
-        self.cases_folder = self.route_test_dir + '/pazy/cases/'
-        self.output_folder = self.route_test_dir + '/pazy/cases/'
-
-    def run_test(self, symmetry_condition, dynamic=False):
-        self.case_name = 'pazy_uinf{:04g}_alpha{:04g}_symmetry_{}'.format(self.u_inf * 10, self.alpha * 10, str(int(symmetry_condition)))
+        self.n_tsteps = 2
+        
+        self.setup_test_folders('pazy')
+       
+    def run_test(self, symmetry_condition, dynamic=False, gust_vanes=False, cs_deflection_file=None):
+        self.case_name = 'pazy_uinf{:04g}_alpha{:04g}_symmetry_{}_gustvanes_{}'.format(self.u_inf * 10, self.alpha * 10, str(int(symmetry_condition)), int(gust_vanes))
 
         gp.generate_pazy(self.u_inf, self.case_name, self.output_folder, self.cases_folder,
                          alpha=self.alpha,
-                         M=self.M,
-                         N=self.N,
+                         M=self.num_chordwise_panels,
+                         N=self.num_spanwise_nodes,
                          Msf=self.Msf,
                          symmetry_condition=symmetry_condition,
                          dynamic=dynamic,
-                         n_tsteps=self.n_tsteps)
-
+                         n_tsteps=self.n_tsteps,
+                         gust_vanes=gust_vanes,
+                         cs_deflection_file=cs_deflection_file)
+     
         self.evaluate_output()
 
     def evaluate_output(self):   
-        pass     
-    #     node_number = self.N / 2 # wing tip node
-
-    #     # Get results in A frame
-    #     tip_displacement = np.loadtxt(self.output_folder + '/' + self.case_name + '/WriteVariablesTime/struct_pos_node{:g}.dat'.format(node_number))
-    #     # current reference from Technion abstract
-    #     ref_displacement = 2.033291e-1  # m
-    #     print("delta z = ", tip_displacement[-1])
-    #     np.testing.assert_almost_equal(tip_displacement[-1], ref_displacement,
-    #                                    decimal=3,
-    #                                    err_msg='Wing tip displacement not within 3 decimal points of reference.',
-    #                                    verbose=True)
+        pass
+
+    def setup_test_folders(self, subfolder):
+        """
+        Set up the case and output directories for the current test,
+        based on the location of the test file that calls this method.
+
+        Args:
+            subfolder (str): Subfolder name like "pazy" or "gust_vanes"
+        """
+        import inspect 
+
+        caller_file = inspect.stack()[1].filename
+        self.route_test_dir = os.path.abspath(os.path.dirname(caller_file))
+        self.cases_folder = os.path.join(self.route_test_dir, subfolder, 'cases')
+        self.output_folder = self.cases_folder
 
     def tearDown(self):
         cases_folder = self.route_test_dir + '/pazy/cases/'
@@ -56,7 +59,7 @@ def tearDown(self):
             import shutil
             shutil.rmtree(cases_folder)
 
-class TestPazyCoupledStatic(TestPazyCoupled):
+class TestPazyCoupledStatic(TestPazyCoupled, unittest.TestCase):
     """
     Test Pazy wing static coupled case and compare against a benchmark result.
 
@@ -70,16 +73,29 @@ def test_static_aoa(self):
     def test_static_aoa_symmetry(self):
         self.run_test(True)
 
-    def evaluate_output(self):     
-        node_number = self.N / 2 # wing tip node
-        # Get results in A frame
-        tip_displacement = np.loadtxt(self.output_folder + '/' + self.case_name + '/WriteVariablesTime/struct_pos_node{:g}.dat'.format(node_number))
-        # current reference from Technion abstract
-        ref_displacement = 2.033291e-1  # m
-        np.testing.assert_almost_equal(tip_displacement[-1], ref_displacement,
-                                       decimal=3,
-                                       err_msg='Wing tip displacement not within 3 decimal points of reference.',
-                                       verbose=True)
+    def evaluate_output(self):  
+        self.assert_tip_displacement_matches_reference(
+            node_index=self.num_spanwise_nodes / 2,
+            ref_displacement=2.033291e-1
+        )
+        
+    def assert_tip_displacement_matches_reference(self, node_index: int, ref_displacement: float):
+        """
+        Check if the tip node displacement matches the reference value.
+
+        Args:
+            node_index (int): Node index along the span (e.g., N/2 for tip)
+            ref_displacement (float): Reference displacement [m]
+        """
+        file_path = os.path.join(self.output_folder, self.case_name,
+                                f'WriteVariablesTime/struct_pos_node{int(node_index)}.dat')
+        tip_displacement = np.loadtxt(file_path)
+
+        np.testing.assert_almost_equal(
+            tip_displacement[-1], ref_displacement, decimal=3,
+            err_msg='Wing tip displacement not within 0.001 m of reference.',
+            verbose=True
+        )   
 
 if __name__ == '__main__':
     unittest.main()