Merge branch 'adri' into 'master'

v1.1 - Dynamic aeroelasticity and composites

See merge request !1

.gitattributes

 *.stp filter=lfs diff=lfs merge=lfs -text
+*.bdf filter=lfs diff=lfs merge=lfs -text

README.md

 # DARTFlo - MPhys
-This repository contains some examples of analysis and optimization cases with [DART](https://gitlab.uliege.be/am-dept/dartflo) and [OpenMDAO](https://openmdao.org/).
+This repository contains some examples of analysis and optimization cases with [DART](https://gitlab.uliege.be/am-dept/dartflo) and [OpenMDAO](https://openmdao.org/)/[MPhys](https://github.com/OpenMDAO/mphys).
 
 ## List of cases
-- acrm: aerostructural analysis of the NASA CRM.
+- crm: aerostructural analysis of the NASA CRM.
 - n12: aerodynamic shape optimiation of a NACA 0012 airfoil.
 - m6: aaerodynamic shape optimiation of the ONERA M6 wing.
-- rae: aerostructural optimization of a benchmark wing.
+- rae_alu: aerostructural optimization of a metalic benchmark wing.
+- rae_cfrp: aerostructural optimization of a composite benchmark wing.
 
 ## Version of solvers
-- [OpenMDAO](https://openmdao.org/): v3.29
+- [OpenMDAO](https://openmdao.org/): v3.30
 - [pyOptSparse](https://github.com/mdolab/pyoptsparse): v2.10.1
-- [MPhys](https://github.com/OpenMDAO/mphys): v1.2.0
+- [MPhys](https://github.com/OpenMDAO/mphys): v1.3.0
+- [OMFlut](https://gitlab.uliege.be/am-dept/omflut): v1.1.0
 - [pySpline](https://github.com/mdolab/pyspline): v1.5.1
 - [pyGeo](https://github.com/mdolab/pygeo): v1.13.0
 - [Gmsh](https://gmsh.info/): v4.10.5
 - [DART](https://gitlab.uliege.be/am-dept/dartflo): v1.2.0
 - [MELD](https://github.com/smdogroup/funtofem): git revision f2b39ef (06/28/2022)
-- [TACS](https://github.com/smdogroup/tacs): v3.3.1
+- [TACS](https://github.com/smdogroup/tacs): v3.7.1
+- [SDPM](https://gitlab.uliege.be/am-dept/sdpm): v1.2.0
+- [PyPk](https://gitlab.uliege.be/am-dept/pypk): v1.1.0

acrm/clean.sh

-rm crm_* cruise_* maneuver_*
\ No newline at end of file

crm/clean.sh

+set -e
+cd workspace
+rm aeroloads_* rae_cruise* rae_maneuver* ffd_* cruise_* maneuver_*

rae/flight.py → crm/flight.py

+from setup import flight_setup
+
 import numpy as np
 import openmdao.api as om
 
+class Flight(om.Group):
+    def setup(self):
+        # mission profile parameters
+        cfg = flight_setup.get()
+        # trim components
+        for i, scn in enumerate(['cruise', 'maneuver']):
+            self.add_subsystem(f'trim_{scn}', Trim(wf=cfg['fixed_weight'], fi=cfg['init_fuel'], q=cfg['dynamic_pressure'][i], s=cfg['ref_area'], g=cfg['gravity']))
+        # fuel volume component
+        self.add_subsystem('volume_match', VolumeMatch(fr=cfg['res_fuel'], fi=cfg['init_fuel']))
+        # fuel burn component
+        self.add_subsystem('fuel_burn', FuelBurn(wf=cfg['fixed_weight'], r=cfg['range'], v=cfg['v_inf'], ct=cfg['sfc'], cd0=cfg['cd_0']))
+        # fuel match component
+        self.add_subsystem('fuel_match', FuelMatch(fi=cfg['init_fuel']))
+
 class Trim(om.ExplicitComponent):
     """Vertical trim equation (mid-cruise): n = L / (Wland+Wfuel/2)
     """

acrm/run.py → crm/run.py

+from setup.struct_dv_setup import struct_comps
+from static import Static
+from struct_dv import StructDvMapper, SmoothnessEvaluatorGrid
+from flight import Flight
+
 import numpy as np
 import openmdao.api as om
 
-from mphys import MultipointParallel
-from dartflo.dart.api.mphys import DartBuilder
-from tacs.mphys import TacsBuilder
-from funtofem.mphys.mphys_meld import MeldBuilder
-from geometry import PyGeoBuilder
-from mphys.scenario_aerostructural import ScenarioAeroStructural
-
-import pygeo_setup
-import dart_setup
-import tacs_setup
-from struct_dv import struct_comps, StructDvMapper, SmoothnessEvaluatorGrid
-from flight import Trim, VolumeMatch, FuelBurn, FuelMatch
-
-_scenarios = ['cruise', 'maneuver']
-
-class Multi(MultipointParallel):
-    def setup(self):
-        # Create the scenarios
-        for i, scn in enumerate(_scenarios):
-            # create the builders
-            args = parseargs()
-            geom = PyGeoBuilder(pygeo_setup.get(i))
-            dart = DartBuilder(dart_setup.get(args.k, args.v, i), scenario='aerostructural', task='optimization', saveAll=True, raiseError=False)
-            tacs = TacsBuilder(tacs_setup.get())
-            meld = MeldBuilder(dart, tacs, isym=1)
-            # create the scenarios
-            csn = self.mphys_add_scenario(scn, ScenarioAeroStructural(aero_builder=dart, struct_builder=tacs, ldxfer_builder=meld, geometry_builder=geom, in_MultipointParallel=True), coupling_nonlinear_solver=om.NonlinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75), coupling_linear_solver=om.LinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75))
-
-    def configure(self):
-        # Configure the multipoint group
-        super().configure()
-        # Configure the scnenarios
-        for i, ssys in enumerate(self.system_iter(recurse=False)):
-            # set surface meshes
-            ssys.geometry.setConstrainedSurface(ssys.aero_mesh.get_triangulated_surface())
-            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val('x_aero0'))
-            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val('x_struct0_masked'))
-            # save sizes
-            sizes = ssys.geometry.getSizes()
-            self.geo_ntwist = sizes['twist']
-            self.geo_nshape = sizes['shape']
-
 class Top(om.Group):
     def setup(self):
         # MPhys
@@ -51,65 +14,46 @@ class Top(om.Group):
         # create mapper for structural DVs
         self.add_subsystem('struct_mapper', StructDvMapper(), promotes=['*'])
         # create multipoint group
-        self.add_subsystem('multi', Multi())
+        self.add_subsystem('static', Static())
         # structural DV smoothness component
         self.add_subsystem('spar_smoothness', SmoothnessEvaluatorGrid(columns=struct_comps['sizes']['spar'], rows=1))
         self.add_subsystem('skin_smoothness', SmoothnessEvaluatorGrid(columns=struct_comps['sizes']['skin'], rows=1))
         # mission profile parameters
-        self.add_subsystem('flight', om.Group())
-        fixed_weight = 170000 # fixed mass, ie landing mass minus wing mass (includes payload and reserve fuel), kg
-        dynamic_pressure = [10990.5, 23568.3] # dynamic pressure (FL370/M0.85, FL200/M0.85), Pa
-        ref_area = 383.7 # reference area, m^2
-        gravity = 9.81 # gravity, m/s^2
-        res_fuel = 15000 # reserve fuel, kg
-        init_fuel = 135000 # baseline fuel weight capacity, kg
-        range = 7725 # range, nm
-        v_inf = 487.5 # KTAS @ FL370, M0.85
-        sfc = 0.53 # SFC, 1/h
-        cd_0 = 0.0137 # CD0 for full aircraft (based on kenway2014)
-        # trim components
-        for i, scn in enumerate(['cruise', 'maneuver']):
-            self.flight.add_subsystem(f'trim_{scn}', Trim(wf=fixed_weight, fi=init_fuel, q=dynamic_pressure[i], s=ref_area, g=gravity))
-        # fuel volume component
-        self.flight.add_subsystem('volume_match', VolumeMatch(fr=res_fuel, fi=init_fuel))
-        # fuel burn component
-        self.flight.add_subsystem('fuel_burn', FuelBurn(wf=fixed_weight, r=range, v=v_inf, ct=sfc, cd0=cd_0))
-        # fuel match component
-        self.flight.add_subsystem('fuel_match', FuelMatch(fi=init_fuel))
+        self.add_subsystem('flight', Flight())
 
     def configure(self):
         # DV - fuel
-        self.dvs.add_output('av_fuel', val=135.0/135.0, desc='fuel available for cruise normalized by baseline fuel capacity')
+        self.dvs.add_output('av_fuel', val=110.0/135.0, desc='fuel available for cruise normalized by baseline fuel capacity')
         """self.add_design_var('av_fuel', lower=0.5, upper=1.5, scaler=1.0)"""
         # DV - aoa
-        aoas = [1.0, 4.0]
-        for i, scn in enumerate(_scenarios):
+        aoas = [3.0, 4.0]
+        for i, scn in enumerate(['cruise', 'maneuver']):
             self.dvs.add_output(f'aoa_{scn}', val=aoas[i], units='deg')
-            self.connect(f'aoa_{scn}', f'multi.{scn}.aoa')
+            self.connect(f'aoa_{scn}', f'static.{scn}.aoa')
             self.add_design_var(f'aoa_{scn}', lower=0.0, upper=6.0, scaler=1.0, units='deg')
         # DV - struct
         for name, value in struct_comps['values'].items():
             self.dvs.add_output(name, val=value)
         self.connect('spar','spar_smoothness.thickness')
         self.connect('skin','skin_smoothness.thickness')
-        for scn in _scenarios:
-            self.connect('dv_struct', f'multi.{scn}.dv_struct')
+        for scn in ['cruise', 'maneuver']:
+            self.connect('dv_struct', f'static.{scn}.dv_struct')
         """for name, value in struct_comps['values'].items():
             self.dvs.add_design_var(name, lower=0.001, upper=0.2, scaler=np.reciprocal(value))"""
         # DV - shape (global)
-        self.dvs.add_output('twist', val=np.array([0] * self.multi.geo_ntwist))
-        for scn in _scenarios:
-          self.connect('twist', f'multi.{scn}.twist')
+        self.dvs.add_output('twist', val=np.array([0] * self.static.geo_ntwist))
+        for scn in ['cruise', 'maneuver']:
+          self.connect('twist', f'static.{scn}.twist')
         """self.add_design_var('twist', lower=-5.0, upper=5.0, scaler=1.0)"""
         # DV - shape (local)
-        self.dvs.add_output('shape', val=np.array([0] * self.multi.geo_nshape))
-        for scn in _scenarios:
-            self.connect('shape', f'multi.{scn}.shape')
+        self.dvs.add_output('shape', val=np.array([0] * self.static.geo_nshape))
+        for scn in ['cruise', 'maneuver']:
+            self.connect('shape', f'static.{scn}.shape')
         """self.add_design_var('shape', lower=-0.1, upper=0.1, scaler=10.0)"""
         # CON - trim
-        for scn in _scenarios:
-            self.connect(f'multi.{scn}.cl', f'flight.trim_{scn}.cl')
-            self.connect(f'multi.cruise.mass', f'flight.trim_{scn}.ww')
+        for scn in ['cruise', 'maneuver']:
+            self.connect(f'static.{scn}.cl', f'flight.trim_{scn}.cl')
+            self.connect(f'static.cruise.mass', f'flight.trim_{scn}.ww')
             self.connect('av_fuel', f'flight.trim_{scn}.fa')
         self.add_constraint('flight.trim_cruise.n', equals=1.0, scaler=1.0)
         self.add_constraint('flight.trim_maneuver.n', equals=2.5, scaler=0.4)
@@ -118,50 +62,35 @@ class Top(om.Group):
         self.connect('flight.fuel_burn.fb', 'flight.fuel_match.fb')
         self.add_constraint('flight.fuel_match.fm', lower=1.0, upper=1.05, scaler=1.0)
         # CON - struct
-        """for scn in _scenarios:
-          self.add_constraint(f'multi.{scn}.ks_fail', upper=1.0, scaler=1.0)"""
-        self.add_constraint(f'multi.maneuver.ks_fail', upper=1.0, scaler=1.0)
+        self.add_constraint(f'static.maneuver.ks_fail', upper=1.0, scaler=1.0)
         """self.add_constraint('spar_smoothness.diff', ref=5e-3, upper=0.0, linear=True)
         self.add_constraint('skin_smoothness.diff', ref=5e-3, upper=0.0, linear=True)"""
         # CON - shape
-        """self.add_constraint('multi.cruise.thckle_con', lower=1.0, scaler=1.0)
-        self.add_constraint('multi.cruise.thckte_con', lower=1.0, scaler=1.0)
-        self.add_constraint('multi.cruise.thcksp_con', lower=0.8, scaler=1.0)
-        self.add_constraint('multi.cruise.le_con', equals=0.0, scaler=1.0, linear=True)
-        self.add_constraint('multi.cruise.te_con', equals=0.0, scaler=1.0, linear=True)"""
+        """self.add_constraint('static.cruise.thckle_con', lower=1.0, scaler=1.0)
+        self.add_constraint('static.cruise.thckte_con', lower=1.0, scaler=1.0)
+        self.add_constraint('static.cruise.thcksp_con', lower=0.8, scaler=1.0)
+        self.add_constraint('static.cruise.le_con', equals=0.0, scaler=1.0, linear=True)
+        self.add_constraint('static.cruise.te_con', equals=0.0, scaler=1.0, linear=True)"""
         # CON - shape (fuel tank volume)
-        self.connect('multi.cruise.vr', 'flight.volume_match.vr', src_indices=[0])
+        self.connect('static.cruise.vr', 'flight.volume_match.vr', src_indices=[0])
         self.connect('av_fuel', 'flight.volume_match.fa')
         """self.add_constraint('flight.volume_match.vc', lower=1.0, scaler=1.0)"""
         # OBJ
-        self.connect('multi.cruise.mass', 'flight.fuel_burn.ww')
-        self.connect('multi.cruise.cl', 'flight.fuel_burn.cl')
-        self.connect('multi.cruise.cd', 'flight.fuel_burn.cd')
+        self.connect('static.cruise.mass', 'flight.fuel_burn.ww')
+        self.connect('static.cruise.cl', 'flight.fuel_burn.cl')
+        self.connect('static.cruise.cd', 'flight.fuel_burn.cd')
         self.add_objective('flight.fuel_burn.fb', scaler=1e-5)
 
-def parseargs():
-    from argparse import ArgumentParser
-    parser = ArgumentParser()
-    parser.add_argument('-v', help='increase output verbosity', action='count', default=0)
-    parser.add_argument('-k', help='nb of threads', type=int, default=1)
-    return parser.parse_args()
-
-def print_time(selapsed):
-    days, rem = divmod(selapsed, 24*60*60)
-    hours, rem = divmod(rem, 60*60)
-    minutes, seconds = divmod(rem, 60)
-    print('Wall-clock time: {:0>2}-{:0>2}:{:0>2}:{:0>2}'.format(int(days),int(hours),int(minutes),int(seconds)))
-
 def main():
+    from utils import change_workdir, print_time
     import time
     from mpi4py import MPI
+    # Change to working directory
+    change_workdir()
     # Setup
-    prob = om.Problem(model=Top(), driver=om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-5, disp=True, maxiter=3, debug_print=['desvars', 'nl_cons', 'ln_cons']))
-    #prob = om.Problem(model=Top(), driver=om.pyOptSparseDriver(optimizer='SLSQP'))
-    #prob.driver.opt_settings['ACC'] = 1e-4
-    #prob.driver.opt_settings['MAXIT'] = 200
+    prob = om.Problem(model=Top(), driver=om.ScipyOptimizeDriver(optimizer='SLSQP', tol=1e-4, disp=True, maxiter=3, debug_print=['desvars', 'nl_cons', 'ln_cons']))
     if MPI.COMM_WORLD.size == 1:
-        prob.model.recording_options['includes'] = ['av_fuel', 'aoa_cruise', 'aoa_maneuver', 'twist', 'dihedral', 'taper', 'sweep', 'shape', 'ribs', 'le_spar', 'te_spar','up_skin', 'lo_skin', 'up_stringer', 'lo_stringer', 'multi.cruise.cl', 'multi.cruise.cd', 'multi.maneuver.cl', 'multi.maneuver.cd', 'multi.cruise.ks_fail', 'multi.maneuver.ks_fail', 'multi.cruise.mass', 'flight.trim_cruise.n', 'flight.trim_maneuver.n', 'flight.volume_match.vc', 'flight.fuel_burn.fb', 'flight.fuel_match.fm']
+        prob.model.recording_options['includes'] = ['av_fuel', 'aoa_cruise', 'aoa_maneuver', 'twist', 'shape', 'ribs', 'le_spar', 'te_spar','up_skin', 'lo_skin', 'up_stringer', 'lo_stringer', 'static.cruise.cl', 'static.cruise.cd', 'static.maneuver.cl', 'static.maneuver.cd', 'static.cruise.ks_fail', 'static.maneuver.ks_fail', 'static.cruise.mass', 'flight.trim_cruise.n', 'flight.trim_maneuver.n', 'flight.volume_match.vc', 'flight.fuel_burn.fb', 'flight.fuel_match.fm']
         prob.model.add_recorder(om.SqliteRecorder('cases.sql'))
     prob.setup(mode='rev')
     om.n2(prob, show_browser=False, outfile='mphys_dart_tacs.html')
@@ -178,16 +107,16 @@ def main():
     ns = []
     kss = []
     dzs = []
-    for scn in _scenarios:
+    for scn in ['cruise', 'maneuver']:
         aoas.append(prob.get_val(f'aoa_{scn}', get_remote=True)[0])
-        cls.append(prob.get_val(f'multi.{scn}.cl', get_remote=True)[0])
-        cds.append(prob.get_val(f'multi.{scn}.cd', get_remote=True)[0])
+        cls.append(prob.get_val(f'static.{scn}.cl', get_remote=True)[0])
+        cds.append(prob.get_val(f'static.{scn}.cd', get_remote=True)[0])
         ns.append(prob.get_val(f'flight.trim_{scn}.n', get_remote=True)[0])
-        kss.append(prob.get_val(f'multi.{scn}.ks_fail', get_remote=True)[0])
-        dzs.append(max(prob.get_val(f'multi.{scn}.u_struct', get_remote=True)[2::3]))
+        kss.append(prob.get_val(f'static.{scn}.ks_fail', get_remote=True)[0])
+        dzs.append(max(prob.get_val(f'static.{scn}.u_struct', get_remote=True)[2::3]))
     # Display
     if MPI.COMM_WORLD.rank == 0:
-        for i, scn in enumerate(_scenarios):
+        for i, scn in enumerate(['cruise', 'maneuver']):
             print(f'{scn}:'.capitalize())
             print('  aoa =', aoas[i], 'deg')
             print('  cl =', cls[i])
@@ -196,7 +125,7 @@ def main():
             print('  KSFailure =', kss[i])
             print('  z-deflection =', dzs[i], 'm')
         print('Geo:')
-        print('  mass =', prob['multi.cruise.mass'][0], 'kg')
+        print('  mass =', prob['static.cruise.mass'][0], 'kg')
         print('  twist =', prob['twist'])
         print('  volume matching =', prob['flight.volume_match.vc'][0])
         print('Perfo:')
@@ -211,13 +140,13 @@ def main():
             matrix[i,0] = i
             matrix[i,1] = case.get_design_vars(scaled=False)['aoa_cruise'][0]
             matrix[i,2] = case.get_design_vars(scaled=False)['aoa_maneuver'][0]
-            matrix[i,3] = case.get_val('multi.cruise.cl')[0]
-            matrix[i,4] = case.get_val('multi.maneuver.cl')[0]
-            matrix[i,5] = case.get_val('multi.cruise.cd')[0]
-            matrix[i,6] = case.get_val('multi.cruise.mass')[0]
+            matrix[i,3] = case.get_val('static.cruise.cl')[0]
+            matrix[i,4] = case.get_val('static.maneuver.cl')[0]
+            matrix[i,5] = case.get_val('static.cruise.cd')[0]
+            matrix[i,6] = case.get_val('static.cruise.mass')[0]
             matrix[i,7] = case.get_constraints(scaled=False)['flight.trim_cruise.n'][0]
             matrix[i,8] = case.get_constraints(scaled=False)['flight.trim_maneuver.n'][0]
-            matrix[i,9] = case.get_constraints(scaled=False)['multi.maneuver.ks_fail'][0]
+            matrix[i,9] = case.get_constraints(scaled=False)['static.maneuver.ks_fail'][0]
             matrix[i,10] = case.get_constraints(scaled=False)['flight.volume_match.vc'][0]
             matrix[i,11] = case.get_constraints(scaled=False)['flight.fuel_match.fm'][0]
             matrix[i,12] = case.get_objectives(scaled=False)['flight.fuel_burn.fb'][0]

acrm/crm.geo → crm/setup/crm.geo

@@ -16,7 +16,7 @@
 
 /* NASA CRM - NASA Common Research Model
  * CAD created by Guillem Battle I Capa
- * Valid mesh obtained with gmsh 4.8.4
+ * Valid mesh obtained with Gmsh 4.10.5
  */
 
 /* Model information
@@ -103,7 +103,7 @@ Characteristic Length {46,11} = msFus; // Fuselage-wake intersection
 // Wing
 Characteristic Length {24, 25,26,27,14} = msWingR; // wing root
 Characteristic Length {22, 20,21,23,9} = msWingY; // wing kink
-Characteristic Length {18,10,16,19} = 2 * msWingT; // wing tip TE
+Characteristic Length {18,10,16,19} = 3 * msWingT; // wing tip TE
 Characteristic Length {17} = msWingT; // wing tip LE
 
 // Tail
@@ -112,8 +112,8 @@ Characteristic Length {1,2,105} = 2 * msTailT; // tail tip TE
 Characteristic Length {104} = msTailT; // tail tip LE
 
 // Mesh algos
-Mesh.Algorithm = 1;
-Mesh.Algorithm3D = 2;
+Mesh.Algorithm = 6;
+Mesh.Algorithm3D = 10;
 Mesh.Optimize = 1;
 Mesh.Smoothing = 10;
 Mesh.SmoothNormals = 1;

acrm/dart_setup.py → crm/setup/dart_setup.py

+from utils import get_path
+
 def get(nthrd, verb, iscn):
     # common parameters
     cfg = {
@@ -5,7 +7,7 @@ def get(nthrd, verb, iscn):
         'Threads' : nthrd, # number of threads
         'Verb' : verb, # verbosity
         # Model (geometry or mesh)
-        'File' : 'crm.geo', # Input file containing the model
+        'File' : get_path('crm.geo', 'setup'), # Input file containing the model
         'Pars' : {'msWingR': 0.121, 'msWingY': 0.073, 'msWingT': 0.028,
                   'msTailR': 0.055, 'msTailT': 0.023,
                   'msFus': 0.628, 'msFar': 36.8}, # parameters for input file model
@@ -43,7 +45,7 @@ def get(nthrd, verb, iscn):
     # define freestream
     if iscn == 0:
         cfg['M_inf'] = 0.85 # cruise Mach
-        cfg['Q_inf'] = 10990.5 # FL370
+        cfg['Q_inf'] = 12091.9 # FL350
     elif iscn == 1:
         cfg['M_inf'] = 0.85 # maneuver Mach
         cfg['Q_inf'] = 23568.3 # FL200

crm/setup/flight_setup.py

+# Mass breakdown for the CRM (based on B777-200ER)
+# B777-200ER data
+#   MTOW     OEW     MZFW     MPAY     FUEL
+#   298t    138t     195t      57t     135t
+# Assumed data for design mission
+#   WING     PAY     RES
+#    20t     35t     15t
+def get():
+    return {
+        'fixed_weight': 168000, # fixed mass, ie landing mass minus wing mass (includes payload and reserve fuel), kg
+        'dynamic_pressure': [12091.9, 23568.3], # dynamic pressure (FL350/M0.85, FL200/M0.85), Pa
+        'ref_area': 383.7, # reference area, m^2
+        'gravity': 9.81, # gravity, m/s^2
+        'res_fuel': 15000, # reserve fuel, kg
+        'init_fuel': 135000, # baseline fuel weight capacity, kg
+        'range': 7500, # range, nm
+        'v_inf': 490.1, # KTAS @ FL370, M0.85
+        'sfc': 0.53, # SFC, 1/h
+        'cd_0': 0.0137 # CD0 for full aircraft (based on kenway2014)
+    }

acrm/pygeo_setup.py → crm/setup/pygeo_setup.py

+from utils import get_path
 import numpy as np
+
 def get(iscn):
     # Parameters that should be retrieved from the setup
     nTwist = 3
@@ -26,7 +28,7 @@ def get(iscn):
     # Configuration
     p = {
         # Input file
-        'geoFile' : 'ffd.xyz',
+        'geoFile' : get_path('ffd.xyz', 'setup'),
         # Disciplines
         'disciplines' : ['aero', 'struct'],
         # Design variables

crm/setup/struct_dv_setup.py

+# Structural components
+# sizes: dict of design variable sizes
+# indices: dict of list mapping a value of a DV to its index in NASTRAN/TACS
+# values: dict of list of dv values
+struct_comps = {'sizes': {'spar': 6,
+                          'ribs': 3,
+                          'stng': 2,
+                          'skin': 6},
+                'indices': {'spar': [0, 1, 2, 3, 4, 5],
+                            'ribs': [6, 16, 7],
+                            'stng': [8, 9],
+                            'skin': [10, 11, 12, 13, 14, 15]},
+                'values': {'spar': [19.05e-3, 15.24e-3, 12.7e-3, 9.525e-3, 7.112e-3, 2.54e-3],
+                           'ribs': [3.048e-3, 2.159e-3, 1.588e-3],
+                           'stng': [4.445e-3, 2.54e-3],
+                           'skin': [6.35e-3, 5.334e-3, 4.572e-3, 3.302e-3, 2.286e-3, 1.27e-3]}
+               }

acrm/tacs_setup.py → crm/setup/tacs_setup.py

+from utils import get_path
 from tacs import elements, constitutive, functions
 import numpy as np
 
 def get():
-    # define elements
-    """def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
-        # Material properties
-        rho = 2780.0  # density, kg/m^3
-        E = 73.1e9  # elastic modulus, Pa
-        nu = 0.33  # poisson's ratio
-        ys = 324.0e6  # yield stress, Pa
-        thickness = 0.003
-        min_thickness = 0.001
-        max_thickness = 0.02
-        # Setup (isotropic) property and constitutive objects
-        prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys)
-        # Set one thickness dv for every component
-        con = constitutive.IsoShellConstitutive(prop, t=thickness, tNum=dvNum, tlb=min_thickness, tub=max_thickness)
-        # For each element type in this component, pass back the appropriate tacs element object
-        transform = None
-        elem = elements.Quad4Shell(transform, con)
-        return elem"""
     # define problem
     def problem_setup(scenario_name, fea_assembler, problem):
         # Mass
@@ -38,14 +21,9 @@ def get():
         # Engine loads
         We = g * 78.0e3 / 12.0 # N per node
         problem.addLoadToNodes([125, 38, 73, 96, 35, 110, 75, 40, 124, 277, 39, 212], np.append(We, [0., 0., 0.]), nastranOrdering=True)
-        # Fuel loads
-        # TODO requires a bdf that "describes" components by ID
-        """compIds = fea_assembler.selectCompIDs(include='SKIN_LOWER') # TODO should select only part of lower skin
-        p = 90000 * g / (0.75 * 383) # N/m^2, TODO crudely approximated by applying W_fuel/(0.75*S_ref) on lower skin, should be variable
-        problem.addPressureToComponents(compIds, p)"""
 
     return {
         #'element_callback': element_callback,
         'problem_setup': problem_setup,
-        'mesh_file': 'crm.bdf'
+        'mesh_file': get_path('crm.bdf', 'setup')
         }

crm/static.py

+from utils import parse_args
+from setup import pygeo_setup, dart_setup, tacs_setup
+
+from dartflo.dart.api.mphys import DartBuilder
+from tacs.mphys import TacsBuilder
+from funtofem.mphys.mphys_meld import MeldBuilder
+from geometry import PyGeoBuilder
+from mphys.scenario_aerostructural import ScenarioAeroStructural
+from mphys import MultipointParallel
+import openmdao.api as om
+
+class Static(MultipointParallel):
+    def setup(self):
+        # Create the scenarios
+        args = parse_args()
+        for i, scn in enumerate(['cruise', 'maneuver']):
+            # create the builders
+            geom = PyGeoBuilder(pygeo_setup.get(i))
+            dart = DartBuilder(dart_setup.get(args.k, args.v, i), scenario='aerostructural', task='optimization', saveAll=False, raiseError=False)
+            tacs = TacsBuilder(**tacs_setup.get())
+            meld = MeldBuilder(dart, tacs, isym=1)
+            # create the scenarios
+            csn = self.mphys_add_scenario(scn, ScenarioAeroStructural(aero_builder=dart, struct_builder=tacs, ldxfer_builder=meld, geometry_builder=geom, in_MultipointParallel=True), coupling_nonlinear_solver=om.NonlinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75), coupling_linear_solver=om.LinearBlockGS(maxiter=10, iprint=2, use_aitken=True, rtol=1e-6, atol=1e-8, aitken_initial_factor=0.75))
+
+    def configure(self):
+        # Configure the multipoint group
+        super().configure()
+        # Configure the scnenarios
+        for i, ssys in enumerate(self.system_iter(recurse=False)):
+            # set surface meshes
+            ssys.geometry.setConstrainedSurface(ssys.aero_mesh.get_triangulated_surface())
+            ssys.geometry.addDisciplinePoints('aero', ssys.aero_mesh.get_val('x_aero0'))
+            ssys.geometry.addDisciplinePoints('struct', ssys.struct_mesh.masker.get_val('x_struct0_masked'))
+            # save sizes
+            sizes = ssys.geometry.getSizes()
+            self.geo_ntwist = sizes['twist']
+            self.geo_nshape = sizes['shape']

acrm/struct_dv.py → crm/struct_dv.py

+from setup.struct_dv_setup import struct_comps 
 import openmdao.api as om
 
-##### OLD BACKUP
-"""struct_comps = { 'names' : ['spar_Z1', 'spar_Z2', 'spar_Z3', 'spar_Z4', 'spar_Z5', 'spar_Z6', 'ribs_IB', 'ribs_MD', 'ribs_OB', 'stng_IB', 'stng_OB', 'skin_Z1', 'skin_Z2', 'skin_Z3', 'skin_Z4', 'skin_Z5', 'skin_Z6'],
-                 'indices' : [0, 1, 2, 3, 4, 5, 6, 16, 7, 8, 9, 10, 11, 12, 13, 14, 15],
-                 'values' : [19.05e-3, 15.24e-3, 12.7e-3, 9.525e-3, 7.112e-3, 2.54e-3, 3.048e-3, 2.159e-3, 1.588e-3, 4.445e-3, 2.54e-3, 6.35e-3, 5.334e-3, 4.572e-3, 3.302e-3, 2.286e-3, 1.27e-3]
-               }
-class StructDvMapper(om.ExplicitComponent):
-    def setup(self):
-        self.ndvs = len(struct_comps['names'])
-        for name in struct_comps['names']:
-            self.add_input(name, shape=1)
-        self.add_output('dv_struct', shape=self.ndvs)
-
-    def compute(self, inputs, outputs):
-        for i in range(self.ndvs):
-            outputs['dv_struct'][struct_comps['indices'][i]] = inputs[struct_comps['names'][i]][0]
-
-    def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
-        if mode == 'fwd':
-            raise NotImplementedError('StructDvMapper - forward mode not implemented!\n')
-        if mode == 'rev':
-            if 'dv_struct' in d_outputs:
-                for i in range(self.ndvs):
-                    d_inputs[struct_comps['names'][i]][0] += d_outputs['dv_struct'][struct_comps['indices'][i]]"""
-#####
-
-# Structural components
-# names: list of design variable names
-# sizes: list of design variable sizes
-# indices: list of list mapping a value of a DV to its index in NASTRAN/TACS
-# values: list of list of dv values
-#struct_comps = {'names': ['spars', 'ribs', 'stng', 'skin'],
-#                'sizes': [6, 3, 2, 6],
-#                'indices': [[0, 1, 2, 3, 4, 5], [6, 16, 7], [8, 9], [10, 11, 12, 13, 14, 15]],
-#                'values': [[19.05e-3, 15.24e-3, 12.7e-3, 9.525e-3, 7.112e-3, 2.54e-3],
-#                           [3.048e-3, 2.159e-3, 1.588e-3],
-#                           [4.445e-3, 2.54e-3],
-#                           [6.35e-3, 5.334e-3, 4.572e-3, 3.302e-3, 2.286e-3, 1.27e-3]]
-#               }
-struct_comps = {'sizes': {'spar': 6,
-                          'ribs': 3,
-                          'stng': 2,
-                          'skin': 6},
-                'indices': {'spar': [0, 1, 2, 3, 4, 5],
-                            'ribs': [6, 16, 7],
-                            'stng': [8, 9],
-                            'skin': [10, 11, 12, 13, 14, 15]},
-#                'values': {'spar': [19.00e-3, 19.00e-3, 19.00e-3, 19.00e-3, 19.00e-3, 19.00e-3],
-#                           'ribs': [10.00e-3, 10.00e-3, 10.00e-3],
-#                           'stng': [10.00e-3, 10.00e-3],
-#                           'skin': [10.00e-3, 10.00e-3, 10.00e-3, 8.00e-3, 8.00e-3, 8.00e-3]}
-                'values': {'spar': [19.05e-3, 15.24e-3, 12.7e-3, 9.525e-3, 7.112e-3, 2.54e-3],
-                           'ribs': [3.048e-3, 2.159e-3, 1.588e-3],
-                           'stng': [4.445e-3, 2.54e-3],
-                           'skin': [6.35e-3, 5.334e-3, 4.572e-3, 3.302e-3, 2.286e-3, 1.27e-3]}
-               }
-
 class StructDvMapper(om.ExplicitComponent):
     """OpenMDAO component to map the user-defined DV to those defined in NASTRAN/TACS
     """
@@ -139,7 +84,7 @@ class SmoothnessEvaluatorGrid(om.ExplicitComponent):
                             d_inputs['thickness'][i]  -= d_outputs['diff'][j+1]
                             j += 2
 
-if __name__ == "__main__":
+if __name__ == '__main__':
     from openmdao.api import Problem
     prob = Problem()
     prob.model.add_subsystem('StructDvDistributor',StructDvMapper())

crm/utils.py

+def parse_args():
+    """Parse command line arguments
+    """
+    from argparse import ArgumentParser
+    parser = ArgumentParser()
+    parser.add_argument('-v', help='increase output verbosity', action='count', default=0)
+    parser.add_argument('-k', help='nb of threads', type=int, default=1)
+    return parser.parse_args()
+
+def get_path(fname, dir_name=''):
+    """Append fname and dir_name to base path (where run.py is)
+    """
+    import os.path
+    return os.path.join(os.path.abspath(os.path.dirname(__file__)), dir_name, fname)
+
+def change_workdir():
+    """Create workspace directory and change to it
+    """
+    import os, os.path
+    # build the name of the workspace directory
+    wdir = os.path.join(os.getcwd(), 'workspace')
+    # create the directory
+    if not os.path.isdir(wdir):
+        print('- creating', wdir)
+        os.makedirs(wdir)
+    # change dir
+    print('- changing to', wdir)
+    os.chdir(wdir)
+
+def print_time(selapsed):
+    """Print elapsed time nicely
+    """
+    days, rem = divmod(selapsed, 24*60*60)
+    hours, rem = divmod(rem, 60*60)
+    minutes, seconds = divmod(rem, 60)
+    print('Wall-clock time: {:0>2}-{:0>2}:{:0>2}:{:0>2}'.format(int(days),int(hours),int(minutes),int(seconds)))