USatm1976Comp now accepts altitude as either geopotential or geodetic (…

…#735)

* USatm1976Comp now accepts altitude as either geopotential (the default and previous implementation) or as geodetic.

* used some timings to confirm that an if statement is still faster than the lambda way of doing this.

docs/dymos_book/examples/min_time_climb/min_time_climb.ipynb

@@ -128,7 +128,7 @@
     "        nn = self.options['num_nodes']\n",
     "\n",
     "        self.add_subsystem(name='atmos',\n",
-    "                           subsys=USatm1976Comp(num_nodes=nn),\n",
+    "                           subsys=USatm1976Comp(num_nodes=nn, h_def='geodetic'),\n",
     "                           promotes_inputs=['h'])\n",
     "\n",
     "        self.add_subsystem(name='aero',\n",

dymos/examples/min_time_climb/min_time_climb_ode.py

@@ -14,7 +14,7 @@ def setup(self):
         nn = self.options['num_nodes']
 
         self.add_subsystem(name='atmos',
-                           subsys=USatm1976Comp(num_nodes=nn),
+                           subsys=USatm1976Comp(num_nodes=nn, h_def='geodetic'),
                            promotes_inputs=['h'])
 
         self.add_subsystem(name='aero',

dymos/models/atmosphere/atmos_1976.py

@@ -1237,6 +1237,8 @@ class USatm1976Comp(om.ExplicitComponent):
     Component model for the United States standard atmosphere 1976 tables.
 
     Data for the model was obtained from http://www.digitaldutch.com/atmoscalc/index.htm.
+    Based on the original model documented in https://www.ngdc.noaa.gov/stp/space-weather/online-publications/
+    miscellaneous/us-standard-atmosphere-1976/us-standard-atmosphere_st76-1562_noaa.pdf
 
     Parameters
     ----------
@@ -1255,11 +1257,19 @@ def initialize(self):
         gas_c = 1716.49  # Gas constant (ft lbf)/(slug R)
         self._K = gamma * gas_c
 
+        self.options.declare('h_def', values=('geopotential', 'geodetic'), default='geopotential',
+                             desc='The definition of altitude provided as input to the component.  If "geodetic",'
+                                  'it will be converted to geopotential based on Equation 19 in the original standard.')
+
     def setup(self):
         """
         Add component inputs and outputs.
         """
         nn = self.options['num_nodes']
+
+        self._geodetic = self.options['h_def'] == 'geodetic'
+        self._R0 = 6_356_766 / 0.3048  # Value of R0 from the original standard (m -> ft)
+
         self.add_input('h', val=1. * np.ones(nn), units='ft')
 
         self.add_output('temp', val=1. * np.ones(nn), units='degR')
@@ -1269,7 +1279,7 @@ def setup(self):
         self.add_output('drhos_dh', val=1. * np.ones(nn), units='slug/ft**4')
         self.add_output('sos', val=1 * np.ones(nn), units='ft/s')
 
-        arange = np.arange(nn)
+        arange = np.arange(nn, dtype=int)
         self.declare_partials(['temp', 'pres', 'rho', 'viscosity', 'drhos_dh', 'sos'], 'h',
                               rows=arange, cols=arange)
 
@@ -1287,6 +1297,11 @@ def compute(self, inputs, outputs):
         table_points = USatm1976Data.alt
         h = inputs['h']
 
+        if self._geodetic:
+            h = h / (self._R0 + h) * self._R0  # Equation 19 from the original standard.
+
+        # From this point forward, h is geopotential altitude (z in the original reference).
+
         idx = np.searchsorted(table_points, h, side='left')
         h_bin_left = np.hstack((table_points[0], table_points))
         dx = h - h_bin_left[idx]
@@ -1319,6 +1334,13 @@ def compute_partials(self, inputs, partials):
         """
         table_points = USatm1976Data.alt
         h = inputs['h']
+        dz_dh = 1.0
+
+        if self._geodetic:
+            dz_dh = (self._R0 / (self._R0 + h)) ** 2
+            h = h / (self._R0 + h) * self._R0  # Equation 19 from the original standard.
+
+        # From this point forawrd, h is geopotential altitude (z in the original reference).
 
         idx = np.searchsorted(table_points, h, side='left')
         h_index = np.hstack((table_points[0], table_points))
@@ -1340,12 +1362,20 @@ def compute_partials(self, inputs, partials):
         coeffs = USatm1976Data.akima_drho[idx]
         d2rho_dh2 = coeffs[:, 1] + dx * (2.0 * coeffs[:, 2] + 3.0 * coeffs[:, 3] * dx)
 
-        partials['temp', 'h'] = dT_dh.ravel()
-        partials['pres', 'h'] = dP_dh.ravel()
-        partials['rho', 'h'] = drho_dh.ravel()
-        partials['viscosity', 'h'] = dvisc_dh.ravel()
-        partials['drhos_dh', 'h'] = d2rho_dh2.ravel()
-        partials['sos', 'h'][...] = 0.5 / np.sqrt(self._K * T) * partials['temp', 'h'] * self._K
+        partials['temp', 'h'][...] = dT_dh.ravel()
+        partials['pres', 'h'][...] = dP_dh.ravel()
+        partials['rho', 'h'][...] = drho_dh.ravel()
+        partials['viscosity', 'h'][...] = dvisc_dh.ravel()
+        partials['drhos_dh', 'h'][...] = d2rho_dh2.ravel()
+        partials['sos', 'h'][...] = (0.5 / np.sqrt(self._K * T) * partials['temp', 'h'] * self._K)
+
+        if self._geodetic:
+            partials['sos', 'h'][...] *= dz_dh
+            partials['temp', 'h'][...] *= dz_dh
+            partials['viscosity', 'h'][...] *= dz_dh
+            partials['rho', 'h'][...] *= dz_dh
+            partials['pres', 'h'][...] *= dz_dh
+            partials['drhos_dh', 'h'][...] *= dz_dh ** 2
 
 
 if __name__ == "__main__":

dymos/models/atmosphere/test/test_atmos.py

@@ -9,7 +9,7 @@
 
 class TestAtmosphere(unittest.TestCase):
 
-    def test_temperature_comp(self):
+    def test_atmos_comp_geopotential(self):
         n = USatm1976Data.alt.size
 
         p = om.Problem(model=om.Group())
@@ -28,14 +28,47 @@ def test_temperature_comp(self):
         rho = p.get_val('atmos.rho', units='slug/ft**3')
         sos = p.get_val('atmos.sos', units='ft/s')
 
-        assert_near_equal(T, USatm1976Data.T, tolerance=1.0E-2)
-        assert_near_equal(P, USatm1976Data.P, tolerance=1.0E-2)
-        assert_near_equal(rho, USatm1976Data.rho, tolerance=1.0E-2)
-        assert_near_equal(sos, USatm1976Data.a, tolerance=1.0E-2)
+        assert_near_equal(T, USatm1976Data.T, tolerance=1.0E-4)
+        assert_near_equal(P, USatm1976Data.P, tolerance=1.0E-4)
+        assert_near_equal(rho, USatm1976Data.rho, tolerance=1.0E-4)
+        assert_near_equal(sos, USatm1976Data.a, tolerance=1.0E-4)
 
         cpd = p.check_partials(method='cs', out_stream=None)
         assert_check_partials(cpd)
 
+    def test_atmos_comp_geodetic(self):
+        n = USatm1976Data.alt.size
+
+        p = om.Problem(model=om.Group())
+
+        ivc = p.model.add_subsystem('ivc', subsys=om.IndepVarComp(), promotes_outputs=['*'])
+        ivc.add_output(name='alt', val=USatm1976Data.alt, units='ft')
+
+        p.model.add_subsystem('atmos', subsys=USatm1976Comp(num_nodes=n, h_def='geodetic'))
+        p.model.connect('alt', 'atmos.h')
+
+        p.setup(force_alloc_complex=True)
+
+        h = USatm1976Data.alt * 0.3048  # altitude data in meters
+        R0 = 6_356_766  # US 1976 std atm R0 in m
+        p.set_val('alt', R0 / (R0 - h) * h, units='m')  # US 1976 std atm geopotential altitude to geodetic (m)
+
+        p.run_model()
+
+        T = p.get_val('atmos.temp', units='degR')
+        P = p.get_val('atmos.pres', units='psi')
+        rho = p.get_val('atmos.rho', units='slug/ft**3')
+        sos = p.get_val('atmos.sos', units='ft/s')
+
+        assert_near_equal(T, USatm1976Data.T, tolerance=1.0E-4)
+        assert_near_equal(P, USatm1976Data.P, tolerance=1.0E-4)
+        assert_near_equal(rho, USatm1976Data.rho, tolerance=1.0E-4)
+        assert_near_equal(sos, USatm1976Data.a, tolerance=1.0E-4)
+
+        with np.printoptions(linewidth=100000):
+            cpd = p.check_partials(method='cs')
+        assert_check_partials(cpd)
+
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()