Add the automatic differentiation multibody solver based on JAX

docs/source/content/casefiles.rst

@@ -319,7 +319,99 @@ Item by item:
 
     This is an optional group to add if correcting the aerodynamic forces using airfoil polars is desired. A polar
     should be included for each airfoil defined. Each entry consists of a 4-column table. The first column corresponds
-    to the angle of attack (in radians) and then the ``C_L``, ``C_D`` and ``C_M``. 
+    to the angle of attack (in radians) and then the ``C_L``, ``C_D`` and ``C_M``.
+
+Multibody file
+--------------
+
+All the aerodynamic data is contained in ``case.mb.h5``.
+
+This file encapsulates both the initial conditions for the multiple bodies, and the constraints between them.
+
+Item by item:
+
+* ``num_bodies``: Number of bodies.
+
+* ``num_constraints``: Number of constraints between bodies.
+
+    The initial conditions for each body and the constraint definitions are defined in groups. The body groups are named
+    as ``body_xx``, where the xx is replaced with a two digit body number starting from 00, e.g. ``body_00``. Each of these
+    groups should have the following items:
+
+* ``FoR_acceleration [6]``: Frame of reference initial acceleration.
+
+    An array of the stacked linear and rotational initial accelerations in the inertial frame.
+
+* ``FoR_velocity [6]``: Frame of reference initial velocity.
+
+    An array of the stacked linear and rotational initial velocities in the inertial frame.
+
+* ``FoR_position [6]``: Frame of reference initial position.
+
+    An array of the stacked linear and rotational initial positions in the inertial frame.
+
+* ``quat [4]``: Frame of reference initial orientation.
+
+    A quaternion describing the initial rotation between the body attached and inertial frames.
+
+* ``body_number``: Body number.
+
+    An integer used to identify the body when creating constraints.
+
+* ``FoR_movement``: Type of frame of reference movement.
+
+    Use "free" to include rigid body motion, or "prescribed" for a clamped body.
+
+The constraint groups are named similarly to the bodies, using ``constraint_xx`` where xx is the two digit constraint number
+starting from 00. Each of these groups should have the following items:
+
+* ``scalingFactor``: Scaling factor.
+
+    This value scales the multibody equations, where generally settings to ``dt^2`` provides acceptable results.
+
+* ``behaviour``: Constraint behaviour.
+
+    This string defines the type of constraint applied to a single or multiple bodies. A wide range of standard lower-pair
+    kinematic joints are available, such as hinge and spherical joints, as well as prescribed rotation joints. A list of the
+    available joints can be found in ``sharpy/structure/utils/lagrangeconstraints.py`` and
+    ``sharpy/structure/utils/lagrangeconstraintsjax.py``, depending on the solver being used. Due to every constraint being
+    different, further parameters depend upon the constraint type used. These parameters are added as variables to the constraint
+    group. Some examples which may be included are listed below:
+
+* ``body_FoR``: Body frame of reference number.
+
+    The number of the body which is constrained by its body attached frame of reference. For example for a double pendulum,
+    this would be the lower link.
+
+* ``body``: Body number.
+
+    The number of the body which is constrained by one of its nodes. For example for a double pendulum,
+    this would be the upper link.
+
+* ``node_in_body``: Node in body.
+
+    This is paired to the ``body`` parameter, and this indicates which node within the body is to be constrained.
+
+* ``rot_axisB [3]``: Rotation axis in the B frame.
+
+    For a hinge constraint, this defines a vector for the hinge axis for the ``body``. This is defined in the material frame.
+
+* ``rot_axisA2 [3]``: Rotation axis in the A frame.
+
+    For a hinge constraint, this defines a vector for the hinge axis for the ``body_FoR``. This is defined in the body attached frame.
+
+* ``controller_id``: Controller ID for using an actuated constraint.
+
+    This should use the same ID as the ``MultibodyController`` used in the ``DynamicCoupled`` simulation, allowing the rotation to be controlled over time.
+
+* ``aerogrid_warp_factor [num_elem, 3]``: Aerodynamic grid warping factor.
+
+    For simulating wings with dynamic sweep by a local rotation around the z axis, the aerodynamic grid needs to warp to
+    account for this. To create a smooth warp, it can be preferred to gradually change the sweep (and therefore chord to
+    maintain the "same" geometry) at elements around the discontinuity. This parameter sets the sweep per node using this
+    parameter to scale the input angle. If not included, it will be ignored when generating the aerodynamic grid and have
+    no effect.
+
 
 Nonlifting Body file
 -----------------
@@ -354,7 +446,7 @@ Item by item:
 
 * ``z_0_ellipse [num_node]``: Vertical offset of the ellipse center from the beam node.
 
-    Is an array with the vertical offset of the center of the elliptical cross-sectoin from the fuselage node.
+    Is an array with the vertical offset of the center of the elliptical cross-section from the fuselage node.
 
 *  ``surface_m [num_surfaces]``: Radial panelling.
 

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/models/aerogrid.py

@@ -193,11 +193,6 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, setting
                 else:
                     global_node_in_surface[i_surf].append(i_global_node)
 
-                # master_elem, master_elem_node = beam.master[i_elem, i_local_node, :]
-                # if master_elem < 0:
-                    # master_elem = i_elem
-                    # master_elem_node = i_local_node
-
                 # find the i_surf and i_n data from the mapping
                 i_n = -1
                 ii_surf = -1
@@ -270,15 +265,28 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, setting
                             raise NotImplementedError(str(self.data_dict['control_surface_type'][i_control_surface]) +
                                 ' control surfaces are not yet implemented')
 
-
+                # add sweep for aerogrid warping in constraint defintition
+                # if no constraint_xx.aerogrid warp factor is provided, this will be ignored
+                ang_warp = 0.
+                if structure_tstep.mb_dict is not None:
+                    if structure_tstep.mb_prescribed_dict is not None:
+                        for i_constraint in range(structure_tstep.mb_dict['num_constraints']):
+                            try:
+                                cst_name = f"constraint_{i_constraint:02d}"
+                                ctrl_id = structure_tstep.mb_dict[cst_name]['controller_id'].decode('UTF-8')
+                                f_warp = structure_tstep.mb_dict[cst_name]['aerogrid_warp_factor'][i_elem, i_local_node]
+                                ang_z = structure_tstep.mb_prescribed_dict[ctrl_id]['delta_psi'][2]
+                                ang_warp += f_warp * ang_z
+                            except KeyError:
+                                continue
 
                 node_info = dict()
                 node_info['i_node'] = i_global_node
                 node_info['i_local_node'] = i_local_node
-                node_info['chord'] = self.data_dict['chord'][i_elem, i_local_node]
+                node_info['chord'] = self.data_dict['chord'][i_elem, i_local_node] / np.cos(ang_warp)
                 node_info['eaxis'] = self.data_dict['elastic_axis'][i_elem, i_local_node]
                 node_info['twist'] = self.data_dict['twist'][i_elem, i_local_node]
-                node_info['sweep'] = self.data_dict['sweep'][i_elem, i_local_node]
+                node_info['sweep'] = self.data_dict['sweep'][i_elem, i_local_node] + ang_warp
                 node_info['M'] = self.dimensions[i_surf, 0]
                 node_info['M_distribution'] = self.data_dict['m_distribution'].decode('ascii')
                 node_info['airfoil'] = self.data_dict['airfoil_distribution'][i_elem, i_local_node]

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