Add ground effect / AIC component refactoring (#348)

* semi-working, off from AVL w GE by a few percent

* Groud effect, integration and component tests, and docs

* Delete test_aero_analysis_ground_effect.py

This file shouldn't be here, not sure how it ended up committed

* Revert mangled file

* Fix duplicate jac index performance issue

* add comments; total deriv check output

* version bump

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = '2.2.2'
+__version__ = '2.3.0'

openaerostruct/aerodynamics/aero_groups.py

@@ -60,12 +60,21 @@ def setup(self):
         # While other components only depends on a single surface,
         # this component requires information from all surfaces because
         # each surface interacts with the others.
+
+        # check for ground effect and if so, promote
+        ground_effect = False
+        for surface in surfaces:
+            if surface.get('groundplane', False):
+                ground_effect = True
+
         if self.options['compressible'] == True:
             aero_states = CompressibleVLMStates(surfaces=surfaces, rotational=rotational)
             prom_in = ['v', 'alpha', 'beta', 'rho', 'Mach_number']
         else:
             aero_states = VLMStates(surfaces=surfaces, rotational=rotational)
             prom_in = ['v', 'alpha', 'beta', 'rho']
+        if ground_effect:
+            prom_in.append('height_agl')
 
         aero_states.linear_solver = om.LinearRunOnce()
 

openaerostruct/aerodynamics/compressible_states.py

@@ -48,8 +48,8 @@ def setup(self):
         num_collocation_points = 0
         for surface in surfaces:
             mesh=surface['mesh']
-            nx = self.nx = mesh.shape[0]
-            ny = self.ny = mesh.shape[1]
+            nx = mesh.shape[0]
+            ny = mesh.shape[1]
             num_collocation_points += (ny - 1) * (nx - 1)
 
         num_force_points = num_collocation_points

openaerostruct/aerodynamics/get_vectors.py

@@ -45,29 +45,37 @@ def setup(self):
             nx = mesh.shape[0]
             ny = mesh.shape[1]
             name = surface['name']
+
+            ground_effect = surface.get('groundplane', False)
+
             vectors_name = '{}_{}_vectors'.format(name, eval_name)
 
             # This is where we handle the symmetry in the VLM method.
             # If it's symmetric, we need to effectively mirror the mesh by
             # accounting for the ghost mesh. We do this by using an artificially
             # larger mesh here.
+
             if surface['symmetry']:
                 actual_ny_size = ny * 2 - 1
             else:
                 actual_ny_size = ny
+            if ground_effect:
+                actual_nx_size = nx * 2
+            else:
+                actual_nx_size = nx
 
-            self.add_input(name + '_vortex_mesh', val=np.zeros((nx, actual_ny_size, 3)), units='m')
-            self.add_output(vectors_name, val=np.ones((num_eval_points, nx, actual_ny_size, 3)), units='m')
+            self.add_input(name + '_vortex_mesh', val=np.zeros((actual_nx_size, actual_ny_size, 3)), units='m')
+            self.add_output(vectors_name, val=np.ones((num_eval_points, actual_nx_size, actual_ny_size, 3)), units='m')
 
             # Set up indices so we can get the rows and cols for the delcare
-            vector_indices = np.arange(num_eval_points * nx * actual_ny_size * 3)
+            vector_indices = np.arange(num_eval_points * actual_nx_size * actual_ny_size * 3)
             mesh_indices = np.outer(
                 np.ones(num_eval_points, int),
-                np.arange(nx * actual_ny_size * 3),
+                np.arange(actual_nx_size * actual_ny_size * 3),
             ).flatten()
             eval_indices = np.einsum('il,jk->ijkl',
                 np.arange(num_eval_points * 3).reshape((num_eval_points, 3)),
-                np.ones((nx, actual_ny_size), int),
+                np.ones((actual_nx_size, actual_ny_size), int),
             ).flatten()
 
             self.declare_partials(vectors_name, name + '_vortex_mesh', val=-1., rows=vector_indices, cols=mesh_indices)
@@ -91,9 +99,15 @@ def compute(self, inputs, outputs):
             mesh_name = name + '_vortex_mesh'
             vectors_name = '{}_{}_vectors'.format(name, eval_name)
 
-            mesh_reshaped = np.einsum('i,jkl->ijkl', np.ones(num_eval_points), inputs[mesh_name])
+            ground_effect = surface.get('groundplane', False)
 
-            if surface['symmetry']:
+            mesh_reshaped = np.einsum('i,jkl->ijkl', np.ones(num_eval_points), inputs[mesh_name])
+            if surface['symmetry'] and ground_effect:
+                eval_points_reshaped = np.einsum('il,jk->ijkl',
+                    inputs[eval_name],
+                    np.ones((2*nx, 2*ny-1)),
+                )
+            elif surface['symmetry']:
                 eval_points_reshaped = np.einsum('il,jk->ijkl',
                     inputs[eval_name],
                     np.ones((nx, 2*ny-1)),

openaerostruct/aerodynamics/states.py

@@ -31,8 +31,8 @@ def setup(self):
         num_collocation_points = 0
         for surface in surfaces:
             mesh=surface['mesh']
-            nx = self.nx = mesh.shape[0]
-            ny = self.ny = mesh.shape[1]
+            nx = mesh.shape[0]
+            ny = mesh.shape[1]
             num_collocation_points += (ny - 1) * (nx - 1)
 
         num_force_points = num_collocation_points

openaerostruct/aerodynamics/tests/test_eval_mtx.py

@@ -3,7 +3,7 @@
 from openmdao.utils.assert_utils import assert_check_partials
 
 from openaerostruct.aerodynamics.eval_mtx import EvalVelMtx
-from openaerostruct.utils.testing import run_test, get_default_surfaces
+from openaerostruct.utils.testing import run_test, get_default_surfaces, get_ground_effect_surfaces
 
 
 class Test(unittest.TestCase):
@@ -33,5 +33,32 @@ def test_assembled_jac(self):
         data = prob.check_partials(compact_print=True, out_stream=None, method='cs', step=1e-40)
         assert_check_partials(data, atol=1e20, rtol=1e-6)
 
+class GroundEffectTest(unittest.TestCase):
+
+    def test(self):
+        surfaces = get_ground_effect_surfaces()
+
+        comp = EvalVelMtx(surfaces=surfaces, num_eval_points=2, eval_name='test_name')
+
+        run_test(self, comp, complex_flag=True)
+
+    def test_assembled_jac(self):
+        surfaces = get_ground_effect_surfaces()
+
+        comp = EvalVelMtx(surfaces=surfaces, num_eval_points=2, eval_name='test_name')
+
+        prob = om.Problem()
+        prob.model.add_subsystem('comp', comp)
+
+        prob.model.linear_solver = om.DirectSolver(assemble_jac=True)
+        prob.model.options['assembled_jac_type'] = 'csc'
+
+        prob.setup(force_alloc_complex=True)
+
+        prob.run_model()
+
+        data = prob.check_partials(compact_print=True, out_stream=None, method='cs', step=1e-40)
+        assert_check_partials(data, atol=1e20, rtol=1e-6)
+
 if __name__ == '__main__':
     unittest.main()

openaerostruct/aerodynamics/tests/test_get_vectors.py

@@ -1,7 +1,7 @@
 import unittest
 
 from openaerostruct.aerodynamics.get_vectors import GetVectors
-from openaerostruct.utils.testing import run_test, get_default_surfaces
+from openaerostruct.utils.testing import run_test, get_default_surfaces, get_ground_effect_surfaces
 
 
 class Test(unittest.TestCase):
@@ -13,6 +13,13 @@ def test(self):
 
         run_test(self, comp)
 
+    def test_groundplane(self):
+        surfaces = get_ground_effect_surfaces()
+
+        comp = GetVectors(surfaces=surfaces, num_eval_points=10, eval_name='test_name')
+
+        run_test(self, comp)
+
 
 if __name__ == '__main__':
     unittest.main()

openaerostruct/aerodynamics/tests/test_vortex_mesh.py

@@ -1,7 +1,7 @@
 import unittest
 
 from openaerostruct.aerodynamics.vortex_mesh import VortexMesh
-from openaerostruct.utils.testing import run_test, get_default_surfaces
+from openaerostruct.utils.testing import run_test, get_default_surfaces, get_ground_effect_surfaces
 
 import numpy as np
 np.set_printoptions(linewidth=200)
@@ -16,6 +16,12 @@ def test(self):
 
         run_test(self, comp)
 
+    def test_groundplane(self):
+        surfaces = get_ground_effect_surfaces()
+
+        comp = VortexMesh(surfaces=surfaces)
+
+        run_test(self, comp, atol=1e6)
 
 if __name__ == '__main__':
     unittest.main()