First PR - working JAX pendulum comparison

Also cleaned up a lot of unneeded comments in some files, and prevents one of the tests from leaving undeleted .vtu files.

docs/source/content/example_notebooks/cantilever/model_static_cantilever.py

@@ -16,8 +16,7 @@ def generate_fem_file(route, case_name, num_elem, deadforce=600e3, followerforce
     length = 5
     num_node_elem=3
     num_node = (num_node_elem - 1)*num_elem + 1
-    # import pdb; pdb.set_trace()
-    angle = 0*np.pi/180.0       # Angle of the beam reference line within the x-y plane of the B frame.
+    angle = 0. * np.pi / 180.       # Angle of the beam reference line within the x-y plane of the B frame.
     x = (np.linspace(0, length, num_node))*np.cos(angle)
     y = (np.linspace(0, length, num_node))*np.sin(angle)
     z = np.zeros((num_node,))
@@ -41,15 +40,13 @@ def generate_fem_file(route, case_name, num_elem, deadforce=600e3, followerforce
                           + [0, 2, 1])
 
     # stiffness array
-    # import pdb; pdb.set_trace()
     num_stiffness = 1
     ea = 4.8e8
     ga = 3.231e8
     gj = 1.0e6
     ei = 9.346e6
     base_stiffness = np.diag([ea, ga, ga, gj, ei, ei])
     stiffness = np.zeros((num_stiffness, 6, 6))
-    # import pdb; pdb.set_trace()
     for i in range(num_stiffness):
         stiffness[i, :, :] = base_stiffness
 

scripts/optimiser/optimiser.py

@@ -265,8 +265,6 @@ def optimiser(in_dict, previous_x, previous_y):
 
     print('FINISHED')
 
-    import pdb; pdb.set_trace()
-
 def local_optimisation(opt, yaml_dict=None, min_method='Powell'):
     x_in = opt.X
     y_in = opt.Y
@@ -280,7 +278,6 @@ def local_optimisation(opt, yaml_dict=None, min_method='Powell'):
                'gtol': 1e-3}
     # scipy.optimize
     local_opt = optimize.minimize(
-                                  # lambda x: rbf_constrained(x, rbf, yaml_dict, opt),
                                   lambda x: gp_constrained(x, opt, yaml_dict),
                                   x0=opt.x_opt,
                                   method=min_method,

setup.py

@@ -145,7 +145,8 @@ def run(self):
         "openpyxl>=3.0.10",
         "lxml>=4.4.1",
         "PySocks",
-        "PyYAML"
+        "PyYAML",
+        "jax",
     ],
     extras_require={
         "docs": [

sharpy/aero/utils/airfoilpolars.py

@@ -11,6 +11,9 @@ class Polar:
 
     def __init__(self):
 
+        self.cm_interp = None
+        self.cd_interp = None
+        self.cl_interp = None
         self.table = None
         self.aoa_cl0_deg = None
 
@@ -34,15 +37,11 @@ def initialise(self, table):
         for ipoint in range(npoints - 1):
             if self.table[ipoint, 1] == 0.:
                 matches.append(self.table[ipoint, 0])
-            elif (self.table[ipoint, 1] < 0. and self.table[ipoint + 1, 1] > 0):
-            # elif ((self.table[ipoint, 1] < 0. and self.table[ipoint + 1, 1] > 0) or
-            #       (self.table[ipoint, 1] > 0. and self.table[ipoint + 1, 1] < 0)):
-                if (self.table[ipoint, 0] <= 0.):
+            elif self.table[ipoint, 1] < 0. and self.table[ipoint + 1, 1] > 0:
+                if self.table[ipoint, 0] <= 0.:
                     matches.append(np.interp(0,
                                              self.table[ipoint:ipoint+2, 1],
                                              self.table[ipoint:ipoint+2, 0]))
-                # else:
-                #     print("WARNING: Be careful negative camber airfoil not supported")
 
         iaoacl0 = 0
         aux = np.abs(matches[0])
@@ -62,7 +61,7 @@ def get_coefs(self, aoa_deg):
         cd = self.cd_interp(aoa_deg)
         cm = self.cm_interp(aoa_deg)
 
-        return cl, cd, cm
+        return cl[0], cd[0], cm[0]
 
     def get_aoa_deg_from_cl_2pi(self, cl):
 
@@ -116,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 cd, cm
+        return float(cd), float(cm)
 
 
 def interpolate(polar1, polar2, coef=0.5):

sharpy/controllers/multibodycontroller.py

@@ -0,0 +1,180 @@
+import numpy as np
+import os
+
+import sharpy.utils.controller_interface as controller_interface
+import sharpy.utils.settings as settings
+
+
+@controller_interface.controller
+class MultibodyController(controller_interface.BaseController):
+    r""" """
+
+    controller_id = "MultibodyController"
+
+    settings_types = dict()
+    settings_default = dict()
+    settings_description = dict()
+    settings_options = dict()
+
+    settings_types["ang_history_input_file"] = "str"
+    settings_default["ang_history_input_file"] = None
+    settings_description["ang_history_input_file"] = "Route and file name of the time history of desired CRV rotation"
+
+    settings_types["ang_vel_history_input_file"] = "str"
+    settings_default["ang_vel_history_input_file"] = ""
+    settings_description["ang_vel_history_input_file"] = ("Route and file name of the time history of desired CRV "
+                                                          "velocity")
+
+    settings_types["psi_dot_init"] = "list(float)"
+    settings_default["psi_dot_init"] = [0., 0., 0.]
+    settings_description["psi_dot_init"] = "Initial rotation velocity of hinge"
+
+    settings_types["dt"] = "float"
+    settings_default["dt"] = None
+    settings_description["dt"] = "Time step of the simulation"
+
+    settings_types["write_controller_log"] = "bool"
+    settings_default["write_controller_log"] = True
+    settings_description["write_controller_log"] = (
+        "Write a time history of input, required input, " + "and control"
+    )
+
+    settings_table = settings.SettingsTable()
+    __doc__ += settings_table.generate(
+        settings_types, settings_default, settings_description, settings_options
+    )
+
+    def __init__(self):
+        self.in_dict = None
+        self.data = None
+        self.settings = None
+
+        self.prescribed_ang_time_history = None
+        self.prescribed_ang_vel_time_history = None
+
+        # Time histories are ordered such that the [i]th element of each
+        # is the state of the controller at the time of returning.
+        # That means that for the timestep i,
+        # state_input_history[i] == input_time_history_file[i] + error[i]
+
+        self.real_state_input_history = list()
+        self.control_history = list()
+
+        self.controller_implementation = None
+        self.log = None
+
+    def initialise(self, data, in_dict, controller_id=None, restart=False):
+        self.in_dict = in_dict
+        settings.to_custom_types(
+            self.in_dict, self.settings_types, self.settings_default
+        )
+
+        self.settings = self.in_dict
+        self.controller_id = controller_id
+
+        # whilst PID control is not here implemented, I have left the remains for if it gets implemented in future
+        if self.settings["write_controller_log"]:
+            folder = data.output_folder + "/controllers/"
+            if not os.path.exists(folder):
+                os.makedirs(folder)
+            self.log = open(folder + self.controller_id + ".log.csv", "w+")
+            self.log.write(
+                ("#" + 1 * "{:>2}," + 6 * "{:>12}," + "{:>12}\n").format(
+                    "tstep",
+                    "time",
+                    "Ref. state",
+                    "state",
+                    "Pcontrol",
+                    "Icontrol",
+                    "Dcontrol",
+                    "control",
+                )
+            )
+            self.log.flush()
+
+        # save input time history
+        try:
+            self.prescribed_ang_time_history = np.loadtxt(
+                self.settings["ang_history_input_file"], delimiter=","
+            )
+        except:
+            try:
+                self.prescribed_ang_time_history = np.load(
+                    self.settings["ang_history_input_file"]
+                )
+            except:
+                raise OSError(
+                    "File {} not found in Controller".format(
+                        self.settings["ang_history_input_file"]
+                    )
+                )
+
+        if self.settings["ang_vel_history_input_file"]:
+            try:
+                self.prescribed_ang_vel_time_history = np.loadtxt(
+                    self.settings["ang_vel_history_input_file"], delimiter=","
+                )
+            except:
+                try:
+                    self.prescribed_ang_vel_time_history = np.load(
+                        self.settings["ang_vel_history_input_file"]
+                    )
+                except:
+                    raise OSError(
+                        "File {} not found in Controller".format(
+                            self.settings["ang_vel_history_input_file"]
+                        )
+                    )
+
+    def control(self, data, controlled_state):
+        r"""
+        Main routine of the controller.
+        Input is `data` (the self.data in the solver), and
+        `currrent_state` which is a dictionary with ['structural', 'aero']
+        time steps for the current iteration.
+
+        :param data: problem data containing all the information.
+        :param controlled_state: `dict` with two vars: `structural` and `aero`
+            containing the `timestep_info` that will be returned with the
+            control variables.
+
+        :returns: A `dict` with `structural` and `aero` time steps and control
+            input included.
+        """
+
+        control_command = self.prescribed_ang_time_history[data.ts - 1, :]
+
+        if self.prescribed_ang_vel_time_history is None:
+            if data.ts == 1:
+                psi_dot = self.settings["psi_dot_init"]
+            else:
+                psi_dot = (
+                    self.prescribed_ang_time_history[data.ts - 1, :]
+                    - self.prescribed_ang_time_history[data.ts - 2, :]
+                ) / self.settings["dt"]
+        else:
+            psi_dot = self.prescribed_ang_vel_time_history[data.ts - 1, :]
+
+        if controlled_state["structural"].mb_prescribed_dict is None:
+            controlled_state["structural"].mb_prescribed_dict = dict()
+        controlled_state["structural"].mb_prescribed_dict[self.controller_id] = {
+            "psi": control_command,
+            "psi_dot": psi_dot,
+        }
+        controlled_state["structural"].mb_prescribed_dict[self.controller_id].update(
+            {"delta_psi": control_command - self.prescribed_ang_time_history[0, :]}
+        )
+
+        return controlled_state, control_command
+
+    def controller_wrapper(
+        self, required_input, current_input, control_param, i_current
+    ):
+        self.controller_implementation.set_point(required_input[i_current - 1])
+        control_param, detailed_control_param = self.controller_implementation(
+            current_input[-1]
+        )
+        return control_param, detailed_control_param
+
+    def __exit__(self, *args):
+        self.log.close()

sharpy/controllers/takeofftrajectorycontroller.py

@@ -143,9 +143,6 @@ def control(self, data, controlled_state):
                     mb_dict[self.settings['controlled_constraint']]
         except KeyError:
             return controlled_state
-        except TypeError:
-            import pdb
-            pdb.set_trace()
 
         if self.controlled_body is None or self.controlled_node is None:
             self.controlled_body = constraint['body_number']

sharpy/generators/polaraeroforces.py

@@ -102,6 +102,9 @@ class PolarCorrection(generator_interface.BaseGenerator):
                                        header_line='This generator takes in the following settings.')
 
     def __init__(self):
+        self.folder = None
+        self.cd_from_cl = None
+        self.list_aoa_cl0 = None
         self.settings = None
 
         self.aero = None
@@ -202,11 +205,9 @@ def generate(self, **params):
                     c_bs = local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord))
                     forces_s = c_bs.T.dot(struct_forces[inode, :3])
                     moment_s = c_bs.T.dot(struct_forces[inode, 3:])
-                    drag_force = forces_s[0]
                     lift_force = forces_s[2]
                     # Compute the associated lift
                     cl = np.sign(lift_force) * np.linalg.norm(lift_force) / coef
-                    cd_sharpy = np.linalg.norm(drag_force) / coef
 
                     if self.cd_from_cl:
                         # Compute the drag from the UVLM computed lift

sharpy/solvers/beamloader.py

@@ -136,9 +136,4 @@ def run(self, **kwargs):
         self.data.structure.generate(self.fem_data_dict, self.settings)
         self.data.structure.dyn_dict = self.dyn_data_dict
 
-        # Change the beam description to the local FoR for multibody
-        # if (self.data.structure.num_bodies > 1):
-        #     self.data.structure.ini_info.whole_structure_to_local_AFoR(self.data.structure)
-        #     self.data.structure.timestep_info[0].whole_structure_to_local_AFoR(self.data.structure)
-
         return self.data

sharpy/solvers/dynamiccoupled.py

@@ -543,7 +543,7 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
                 state = {'structural': structural_kstep,
                          'aero': aero_kstep}
                 for k, v in self.controllers.items():
-                    state = v.control(self.data, state)
+                    state, control = v.control(self.data, state)
                     # this takes care of the changes in options for the solver
                     structural_kstep, aero_kstep = self.process_controller_output(
                         state)
@@ -751,8 +751,6 @@ def convergence(self, k, tstep, previous_tstep,
         """
         # check for non-convergence
         if not all(np.isfinite(tstep.q)):
-            import pdb
-            pdb.set_trace()
             raise Exception(
                 '***Not converged! There is a NaN value in the forces!')