Merge branch 'main' of https://github.com/SteveDoyle2/pyNastran

Author: SteveDoyle2

pyNastran/f06/dev/flutter/utils.py

@@ -15,7 +15,8 @@
 def load_f06_op2(f06_filename: str, log: SimpleLogger,
                  in_units: str,
                  out_units: str,
-                 use_rhoref: bool) -> tuple[OP2, dict[int, FlutterResponse]]:
+                 use_rhoref: bool,
+                 make_alt: bool=False) -> tuple[OP2, dict[int, FlutterResponse]]:
     """
     load a Vg-Vf plot from:
      - OP2 / F06
@@ -43,6 +44,7 @@ def load_f06_op2(f06_filename: str, log: SimpleLogger,
                 f06_units=in_units_dict,
                 out_units=out_units_dict,
                 use_rhoref=use_rhoref,
+                make_alt=make_alt,
                 log=log)
         except Exception as e:
             log.error(str(e))

pyNastran/f06/flutter_response.py

@@ -575,10 +575,6 @@ def _set_pknl_results(self,
 
         vel_units_in = in_units['velocity']
         q_units_in = in_units['dynamic_pressure']
-        if _is_q_units_consistent(density_units_in, vel_units_in, q_units_in):
-            q = 0.5 * rho * vel**2
-        else:  # pragma: no cover
-            raise NotImplementedError((density_units_in, vel_units_in, q_units_in))
 
         # eas  = (2 * q / rho_ref)**0.5
         # eas = V * sqrt(rho / rhoSL)
@@ -590,27 +586,35 @@ def _set_pknl_results(self,
         altitude_units = in_units['altitude']
 
         # print('density_units_in=%r density_units2=%r' % (density_units_in, density_units2))
-        kdensityi = convert_density(1., density_units_in, 'slug/ft^3')
-
         resultsi = results[:, :, :9]
         assert resultsi.shape[2] == 9, resultsi.shape
 
-        rho_in_slug_ft3 = rho * kdensityi
+        rho_in_slug_ft3 = rho_ref * rho
         ft_to_alt_unit = convert_altitude(1., 'ft', altitude_units)
         if self.make_alt:
+            nrho = len(rho)
             alt_ft = []
             for idensity, densityi in enumerate(rho_in_slug_ft3.ravel()):
                 try:
                     alt_fti = get_alt_for_density(densityi, density_units='slug/ft^3',
-                                                  alt_units='ft', nmax=20)
+                                                  alt_units='ft', nmax=50)
                 except Exception:
-                    raise RuntimeError(f'failed to find altitude for density[{idensity}]='
-                                       f'{rho.ravel()[idensity]:g}; density_units_in={density_units_in!r}')
+                    raise RuntimeError(f'Case {idensity+1}/{nrho}: failed to find altitude for density='
+                                       f'{rho.ravel()[idensity]:g}; density_units_in={density_units_in!r} ({densityi:g} slug/ft^3)\n'
+                                       '  output_rho_range = rho_ref * input_rho_range\n'
+                                       f'  rho_ref:             {rho_ref:g} {density_units_in}\n'
+                                       f'  input_rho_range:     [{rho.min():g}, {rho.max()}]\n'
+                                       f'  => output_rho_range: [{rho_in_slug_ft3.min():g}, {rho_in_slug_ft3.max()}] slug/ft^3')
                 alt_ft.append(alt_fti)
             alt = np.array(alt_ft, dtype=rho.dtype).reshape(vel.shape) * ft_to_alt_unit
         else:
             alt = np.full(vel.shape, np.nan, dtype=vel.dtype)
 
+        # get dynamic pressure
+        if _is_q_units_consistent(density_units_in, vel_units_in, q_units_in):
+            q = 0.5 * rho * vel**2
+        else:  # pragma: no cover
+            raise NotImplementedError((density_units_in, vel_units_in, q_units_in))
         results2 = np.dstack([resultsi, eas, q, alt])
         # results2[:, :, self.idensity] = rho
         return results2
@@ -809,6 +813,9 @@ def get_flutter_crossings(self,
             easi, dampi, freqi = remove_eas_range(
                 easi, dampi, freqi, eas_range)
 
+            if len(freqi) == 0:
+                continue
+
             for damping_targeti, damping_requiredi in damping_crossings_dict.items():
                 if dampi.max() < damping_requiredi:
                     continue
@@ -1796,6 +1803,8 @@ def plot_vg_vf(self, fig=None, damp_axes=None, freq_axes=None,
             vel_calc, damping_calc, freq_calc = remove_eas_range(
                 vel, damping, freq, eas_range)
 
+            if len(freq) == 0:
+                continue
             jcolor, color, linestyle2, symbol2, texti, is_removedi = _increment_jcolor(
                 mode, jcolor, color, linestyle, symbol,
                 freq, damping, freq_tol=freq_tol, freq_tol_remove=freq_tol_remove,
@@ -3653,14 +3662,22 @@ def _get_divergence(self,
         dampi = self.results[imode, :, self.idamping].flatten()
         easi = self.results[imode, :, self.ieas].flatten()
         freqi = self.results[imode, :, self.ifreq].flatten()
-
         easi, dampi, freqi = remove_excluded_points(
             easi, dampi, freqi, point_removal)
         easi, dampi, freqi = remove_eas_range(
             easi, dampi, freqi, eas_range)
 
+        if len(freqi) == 0:
+            continue
+
         for freq_targeti, freq_requiredi in freq_crossings_dict.items():
-            if freqi.min() > freq_requiredi:
+            try:
+                freqi_min = freqi.min()
+            except:
+                warnings.warn(f'mode={mode}; freqi={freqi}')
+                raise
+
+            if freqi_min > freq_requiredi:
                 continue
 
             # dfreq sign is flipped because get_zero_crossings

pyNastran/f06/parse_flutter.py

@@ -39,6 +39,7 @@ def make_flutter_response(f06_filename: PathLike,
                           f06_units=None, out_units=None,
                           use_rhoref: bool=False,
                           read_flutter: bool=True,
+                          make_alt: bool=True,
                           log: Optional[SimpleLogger]=None) -> tuple[dict[int, FlutterResponse],
                                                                      dict[str, Any]]:
     """
@@ -57,7 +58,10 @@ def make_flutter_response(f06_filename: PathLike,
         False: assume the density in the table is absolute density
         True: assume the density should be defined by sea level density,
               so density is a density ratio
-
+    read_flutter : bool; default=True
+        don't parse the flutter to extract eigenvalues/mass
+    make_alt : bool; default=True
+        make altitude an output variable
 
     Returns
     -------
@@ -225,7 +229,8 @@ def make_flutter_response(f06_filename: PathLike,
                         in_units=f06_units,
                         use_rhoref=use_rhoref,
                         eigenvector=eigenvectors_array,
-                        eigr_eigi_velocity=eigr_eigi_velocity)
+                        eigr_eigi_velocity=eigr_eigi_velocity,
+                        make_alt=make_alt)
                     response.set_out_units(out_units)
                     #_remove_neutrinos(response, log)
                     flutters[subcase] = response
@@ -426,7 +431,8 @@ def make_flutter_response(f06_filename: PathLike,
                 in_units=f06_units,
                 use_rhoref=use_rhoref,
                 eigenvector=eigenvectors_array,
-                eigr_eigi_velocity=eigr_eigi_velocity)
+                eigr_eigi_velocity=eigr_eigi_velocity,
+                make_alt=make_alt)
             response.set_out_units(out_units)
             flutters[subcase] = response
 

pyNastran/utils/atmosphere.py

@@ -38,7 +38,7 @@
 
 
 def get_alt_for_density(density: float, density_units: str='slug/ft^3',
-                        alt_units: str='ft', nmax: int=20, tol: float=5.) -> float:
+                        alt_units: str='ft', nmax: int=20, tol: float=0.1) -> float:
     """
     Gets the altitude associated with a given air density.
 
@@ -52,7 +52,7 @@ def get_alt_for_density(density: float, density_units: str='slug/ft^3',
         sets the units for the output altitude; ft, m, kft
     nmax : int; default=20
         max number of iterations for convergence
-    tol : float; default=5.
+    tol : float; default=0.1
         tolerance in alt_units
 
     Returns
@@ -71,15 +71,20 @@ def get_alt_for_density(density: float, density_units: str='slug/ft^3',
 
     # Newton's method
     alt1 = alt_old
+    log_density = np.log(density)
     while abs(alt_final - alt_old) > tol and n < nmax:
         alt_old = alt_final
         alt1 = alt_old
         alt2 = alt_old + dalt
         rho1: float = atm_density(alt1, density_units=density_units)
         rho2: float = atm_density(alt2, density_units=density_units)
-        m = dalt / (rho2 - rho1)
-        alt_final = m * (density - rho1) + alt1
+        # natural log
+        log_rho1 = np.log(rho1)
+        log_rho2 = np.log(rho2)
+        m = dalt / (log_rho2 - log_rho1)
+        alt_final = m * (log_density - log_rho1) + alt1
         n += 1
+
     if abs(alt_final - alt_old) > tol:
         raise RuntimeError(f'Did not converge; Check your units; n=nmax={nmax}\n'
                            f'target alt={alt_final} alt_current={alt1}')
@@ -91,7 +96,7 @@ def get_alt_for_eas_with_constant_mach(equivalent_airspeed: float,
                                        mach: float,
                                        velocity_units: str='ft/s',
                                        alt_units: str='ft',
-                                       nmax: int=20, tol: float=5.) -> float:
+                                       nmax: int=20, tol: float=0.1) -> float:
     """
     Gets the altitude associated with a equivalent airspeed.
 
@@ -131,6 +136,7 @@ def get_alt_for_eas_with_constant_mach(equivalent_airspeed: float,
     k = np.sqrt(T0 / p0)
     #eas = a * mach * sqrt((p * T0) / (T * p0)) = a * mach * sqrt(p / T) * k
 
+    log_eas = np.log(equivalent_airspeed)
     # Newton's method
     while abs(alt_final - alt_old) > tol and n < nmax:
         alt_old = alt_final
@@ -144,8 +150,10 @@ def get_alt_for_eas_with_constant_mach(equivalent_airspeed: float,
         sos2 = np.sqrt(1.4 * R * T2)
         eas1 = sos1 * mach * np.sqrt(press1 / T1) * k
         eas2 = sos2 * mach * np.sqrt(press2 / T2) * k
-        m = dalt / (eas2 - eas1)
-        alt_final = m * (equivalent_airspeed - eas1) + alt1
+        log_eas1 = np.log(eas1)
+        log_eas2 = np.log(eas2)
+        m = dalt / (log_eas2 - log_eas1)
+        alt_final = m * (log_eas - log_eas1) + alt1
         n += 1
 
     if n > nmax - 1:
@@ -156,7 +164,7 @@ def get_alt_for_eas_with_constant_mach(equivalent_airspeed: float,
 
 def get_alt_for_q_with_constant_mach(q: float, mach: float,
                                      pressure_units: str='psf', alt_units: str='ft',
-                                     nmax: int=20, tol: float=5.) -> float:
+                                     nmax: int=20, tol: float=0.1) -> float:
     """
     Gets the altitude associated with a dynamic pressure.
 
@@ -172,7 +180,7 @@ def get_alt_for_q_with_constant_mach(q: float, mach: float,
         the altitude units; ft, kft, m
     nmax : int; default=20
         max number of iterations for convergence
-    tol : float; default=5.
+    tol : float; default=0.1
         tolerance in alt_units
 
     Returns
@@ -189,7 +197,7 @@ def get_alt_for_q_with_constant_mach(q: float, mach: float,
 
 def get_alt_for_pressure(pressure: float,
                          pressure_units: str='psf', alt_units: str='ft',
-                         nmax: int=20, tol: float=5.) -> float:
+                         nmax: int=20, tol: float=0.1) -> float:
     """
     Gets the altitude associated with a pressure.
 
@@ -203,7 +211,7 @@ def get_alt_for_pressure(pressure: float,
         the altitude units; ft, kft, m
     nmax : int; default=20
         max number of iterations for convergence
-    tol : float; default=5.
+    tol : float; default=0.1
         altitude tolerance in alt_units
 
     Returns
@@ -220,14 +228,17 @@ def get_alt_for_pressure(pressure: float,
     n = 0
 
     # Newton's method
+    log_pressure = np.log(pressure)
     while abs(alt_final - alt_old) > tol and n < nmax:
         alt_old = alt_final
         alt1 = alt_old
         alt2 = alt_old + dalt
         press1 = atm_pressure(alt1)
         press2 = atm_pressure(alt2)
-        m = dalt / (press2 - press1)
-        alt_final = m * (pressure - press1) + alt1
+        log_press1 = np.log(press1)
+        log_press2 = np.log(press2)
+        m = dalt / (log_press2 - log_press1)
+        alt_final = m * (log_pressure - log_press1) + alt1
         n += 1
 
     if n > nmax - 1:
@@ -243,7 +254,7 @@ def get_alt_for_pressure(pressure: float,
 def get_alt_for_mach_eas(mach: float,
                          eas: float, eas_units: str='knots',
                          alt_units: str='ft',
-                         nmax: int=20, tol: float=5.) -> float:
+                         nmax: int=20, tol: float=0.1) -> float:
     """
     Gets the altitude associated with an equivalent airspeed.
 
@@ -259,7 +270,7 @@ def get_alt_for_mach_eas(mach: float,
         the altitude units; ft, kft, m
     nmax : int; default=20
         max number of iterations for convergence
-    tol : float; default=5.
+    tol : float; default=0.1
         altitude tolerance in alt_units
 
     Returns
@@ -276,17 +287,20 @@ def get_alt_for_mach_eas(mach: float,
     n = 0
 
     # Newton's method
+    log_eas = np.log(eas)
     while abs(alt_final - alt_old) > tol and n < nmax:
         alt_old = alt_final
         alt1 = alt_old
         alt2 = alt_old + dalt
         eas1 = atm_equivalent_airspeed(alt1, mach, alt_units='ft', eas_units='ft/s')
         eas2 = atm_equivalent_airspeed(alt2, mach, alt_units='ft', eas_units='ft/s')
+        log_eas1 = np.log(eas1)
+        log_eas2 = np.log(eas2)
 
         # y = (y2-y1)/(x2-x1)*(x-x1) + y1
         # y = m * (x-x1) + y1
-        m = dalt / (eas2 - eas1)
-        alt_final = m * (eas - eas1) + alt1
+        m = dalt / (log_eas2 - log_eas1)
+        alt_final = m * (log_eas - log_eas1) + alt1
         n += 1
 
     if n > nmax - 1:
@@ -1236,6 +1250,7 @@ def make_flfacts_eas_sweep_constant_mach(mach: float,
     #
     # then pressure in a sane unit (e.g., psi/psf/Pa)
     # without a wacky conversion
+    assert eas.min() > 0, eas
     eas_fts = convert_velocity(eas, eas_units, 'ft/s')
     rho0_english: float = atm_density(
         0., R=1716., alt_units=alt_units,
@@ -1250,6 +1265,7 @@ def make_flfacts_eas_sweep_constant_mach(mach: float,
         alt_fti = get_alt_for_pressure(pressure_psfi, pressure_units='psf',
                                        alt_units='ft', nmax=30, tol=1.)
         rhoi = atm_density(alt_fti, R=1716., alt_units='ft', density_units=density_units)
+        assert np.isfinite(alt_ft[i]), alt_fti
         rho[i] = rhoi
         alt_ft[i] = alt_fti
     alt = convert_altitude(alt_ft, 'ft', alt_units)