Add CRM files. Add max number of GMRES iterations in API. Add CREDITS.

CREDITS.md

+DARTFlo (DART) is Copyright (C) of University of Liège.
+
+The main author of DART is Adrien Crovato.
+
+Direct and indirect contributions have been made to DART by the following people:
+- Romain Boman
+- Luc Papeleux
+- Amaury Bilocq, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195).
+- Paul Dechamps, [Improvement of the viscous-inviscid interaction method implemented in DARTFlo](http://hdl.handle.net/2268.2/13886).
+- Guillaume Brian, [Investigation of various transonic stabilisation techniques for full potential flow calculations in DARTFlo](http://hdl.handle.net/2268.2/14510).
+- Guillem Battle I Capa, Aerodynamic and aeroelastic computations of full aircraft configurations.

README.md

@@ -11,10 +11,10 @@ DART is currently capable of rapidly solving steady transonic flows on arbitrary
   - steady transonic flows
   - viscous-inviscid interaction (2D only)
 * Numerical methods
-  - linear algrebra using [Eigen](http://eigen.tuxfamily.org/) and [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html#gs.a3o4w5) or [openBLAS](https://www.openblas.net/)
+  - linear algrebra using [Eigen](http://eigen.tuxfamily.org/)
   - direct (forward) and adjoint (backward) modes
   - nonlinear Newton solver with Bank&Rose line search
-  - linear Intel PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) solvers
+  - linear [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html) PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) solvers
 * Interfaced with
   - [VTK](https://vtk.org/)
   - [CUPyDO](https://github.com/ulgltas/CUPyDO)
@@ -28,6 +28,8 @@ Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preli
 
 Crovato Adrien, [DARTFlo - Discrete Adjoint for Rapid Transonic Flows](http://hdl.handle.net/2268/264284), Technical note, University of Liège, 2021.
 
-Crovato A., Boman R., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
+Crovato A., et al. [A discrete adjoint full potential formulation for fast aerostructural optimization in preliminary aircraft design](), Submitted to Aerospace Science and Technology, 2022.
 
-Bilocq Amaury, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195), Master thesis, University of Liège, 2020.
+Crovato A., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
+
+Crovato A., et al., [Fast Full Potential Based Aerostructural Optimization Calculations for Preliminary Aircraft Design](http://hdl.handle.net/2268/292456), 19th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2022.

dart/api/core.py

@@ -238,7 +238,7 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
         gfil = cfg['G_fill'] if 'G_fill' in cfg else 2
         grst = cfg['G_restart'] if 'G_restart' in cfg else 50
         gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5
-        linsol = tbox.Gmres(gfil, 1e-6, grst, gtol)
+        linsol = tbox.Gmres(gfil, 1e-6, grst, gtol, 200)
     else:
         raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
     # initialize the nonlinear (outer) solver

dart/models/crm.geo

+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/* NASA CRM - NASA Common Research Model
+ * CAD created by Guillem Battle I Capa
+ * Valid mesh obtained with gmsh 4.8.4
+ */
+
+/* Model information
+ * Length of fuselage: 62.79 m
+ * Chord length at wing root: 12.09 m
+ * Chord length at Yehudi break: 7.30 m
+ * Chord length at wing tip: 2.76 m
+ * Chord length at tail root: 5.45 m
+ * Chord length at tail tip: 2.33 m
+ * Reference wing area: 1/2 * 383.69 m^2
+ * Reference chord length: 7.01 m
+ */
+
+/// Constants
+scl = 2; // scaling factor for lifting surfaces mesh size (testing purpose)
+DefineConstant[ msFar = {36.8, Name "Mesh size on farfield" },
+                msFus = {0.628, Name "Mesh size on fuselage" },
+                msWingR = {0.121 * scl, Name "Mesh size on wing root" },
+				msWingY = {0.073 * scl, Name "Mesh size on wing Yehudi break" },
+				msWingT = {0.028 * scl, Name "Mesh size on wing tip" },
+				msTailR = {0.055 * scl, Name "Mesh size on tail root" },
+                msTailT = {0.023 * scl, Name "Mesh size on tail tip" } ];
+
+/// Geometry definition
+// Import
+Geometry.OCCAutoFix = 1;
+Geometry.OCCFixDegenerated = 1;
+Geometry.OCCFixSmallEdges = 1;
+Geometry.OCCFixSmallFaces = 1;
+Geometry.OCCSewFaces = 1;
+Geometry.Tolerance = 1e-5; // 0.01 mm
+Geometry.OCCTargetUnit = "M";
+Merge "crm.stp";
+Coherence;
+
+// Re-orient
+ReverseMesh Surface{1:48}; //To reverse the direction of the mesh normal vectors
+
+// Volume
+Surface Loop(1) = {4:48};
+Volume(1) = {1};
+
+// Embedded
+Line{3} In Surface{16};
+Line{7} In Surface{16};
+Line{15} In Surface{16};
+Surface{1} In Volume{1};
+Surface{2} In Volume{1};
+Surface{3} In Volume{1};
+
+// Physical groups
+Physical Line("wakeExW") = {10,11,12,13};
+Physical Line("teW") = {9,14};
+Physical Line("wakeTipW") = {10};
+Physical Line("wakeExT") = {1,2,5};
+Physical Line("teT") = {6};
+Physical Line("wakeTipT") = {1,2};
+Physical Surface("symmetry") = {19,20,21,22};
+Physical Surface("downstream") = {16};
+Physical Surface("upstream") = {18};
+Physical Surface("farfield") = {14,15,17};
+Physical Surface("fuselage") = {13,23:45};
+Physical Surface("wing") = {4,5,6,11,12};
+Physical Surface("wing_") = {7,8,9,10};
+Physical Surface("tail") = {47,48};
+Physical Surface("tail_") = {46};
+Physical Surface("wakeW") = {2,3};
+Physical Surface("wakeT") = {1};
+Physical Volume("field") = {1};
+
+/// Mesh definition
+// Farfield
+Characteristic Length {34,35,36,37,38,39,40,43,44,45} = msFar; // farfield boundaries
+Characteristic Length {3,4,7,8,15,41,42} = msFar; // wake
+
+// Fuselage
+Characteristic Length {11,12,13,51,61,50,73,74,62,63,64,65,47,29,70,46} = msFus; // body
+Characteristic Length {51,52,53,54,55,56,59,60,89,85,86,90,91,92} = 0.5 * msFus; // nose and windshield
+Characteristic Length {57,58,87,88} = 0.2 * msFus; // surface between windshield
+Characteristic Length {49,72,28,33,32,24,31,30,75,71,48,76,69,47} = 0.75 * msFus; // wingbox
+Characteristic Length {5,66,67,68} = 0.25 * msFus; // APU
+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 {17} = msWingT; // wing tip LE
+
+// Tail
+Characteristic Length {79:84,93:99, 100:102,77,78,103,6} = msTailR; // Tail root
+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.Optimize = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+Mesh 1;
+Geometry.Tolerance = 1e-6;
+Coherence Mesh;

dart/validation/crm.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute the flow around the NASA CRM (Common Research Model) at CL = 0.5 and Mach 0.85
+# CAD created by Guillem Battle I Capa
+# Valid mesh obtained with gmsh 4.8.4
+
+try:
+   import tboxVtk
+   import vtk
+   hasVTK = True
+except:
+   hasVTK = False
+
+def cfg():
+    from fwk.wutils import parseargs
+    import os.path
+    # Mesh size
+    wrms = 0.121 # wing root trailing edge mesh size
+    wyms = 0.073 # wing root trailing edge mesh size
+    wtms = 0.028 # wing root trailing edge mesh size
+    trms = 0.055 # wing root trailing edge mesh size
+    ttms = 0.023 # wing root trailing edge mesh size
+    fums = 0.628 # fuselage mesh size
+    fams = 36.8 # farfield mesh size
+    # Parse
+    args = parseargs()
+    # Parameters
+    return {
+    # Options
+    'Threads' : args.k, # number of threads
+    'Verb' : args.verb + 2, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/crm.geo'), # Input file containing the model
+    'Pars' : {'msWingR' : wrms, 'msWingY' : wyms, 'msWingT' : wtms,
+              'msTailR' : trms, 'msTailT' : ttms,
+              'msFus' : fums, 'msFar' : fams}, # parameters for input file model
+    'Dim' : 3, # problem dimension (2 or 3)
+    'Format' : 'vtk', # save format (vtk or gmsh)
+    # Markers...
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wings' : ['wing', 'tail'], # LIST of names of physical groups containing the lifting surface boundary (first element will be the body of interest for aerostructural and optimization)
+    'Wakes' : ['wakeW', 'wakeT'], # LIST of names of physical group containing the wake
+    'WakeTips' : ['wakeTipW', 'wakeTipT'], # LIST of names of physical group containing the edge of the wake
+    'WakeExs' : ['wakeExW', 'wakeExT'], # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
+    'Tes' : ['teW', 'teT'], # LIST of names of physical group containing the trailing edge
+    'Fuselage' : 'fuselage', # name of physical group containing the fuselage boundaries
+    'Symmetry' : 'symmetry', # name of physical group containing the symmetry boundaries
+    # Freestream
+    'M_inf' : 0.85, # freestream Mach number
+    'AoA' : 1.2, # freestream angle of attack [deg] (optional, default=0)
+    'AoS' : 0.0, # freestream angle of sideslip [deg] (optional, default=0)
+    'Q_inf' : 0.0, # freesteam dynamic pressure (only required for aerostructural computations)
+    # Geometry
+    'S_ref' : 383.69/2, # reference surface length
+    'c_ref' : 7.005, # reference chord length
+    'x_ref' : 33.667, # x-coordinate of reference point for moment computation
+    'y_ref' : 11.906, # y-coordinate of reference point for moment computation
+    'z_ref' : 4.520, # z-coordinate of reference point for moment computation
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (PARDISO, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner (optional, default=2)
+    'G_tol' : 1e-5, # tolerance for GMRES (optional, default=1e-5)
+    'G_restart' : 50, # restart for GMRES (optional, default=50)
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 25 # maximum number of iterations for nonlinear solver
+    }
+
+def main():
+    from dart.api.core import initDart
+    from dart import utils
+    from fwk import Timers
+    from fwk.testing import CTest, CTests
+    import numpy as np
+    import matplotlib.pyplot as plt
+    # pre
+    tms = Timers()
+    tms['pre'].start()
+    _dart = initDart(cfg(), scenario='aerodynamic', task='analysis')
+    msh = _dart['msh']
+    sol = _dart['sol']
+    wrt = _dart['wrt']
+    tms['pre'].stop()
+    # solve
+    tms['sol'].start()
+    sol.run()
+    sol.save(wrt)
+    tms['sol'].stop()
+    # post
+    tms['pos'].start()
+    if hasVTK:
+        utils.writeSlices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
+        for i in range(10):
+            data = np.loadtxt(f'{msh.name}_slice_{i}.dat', delimiter=',', skiprows=1) # read
+            plt.plot(data[:,3], data[:,4], lw=2) # plot results
+        font = {'family': 'serif', 'color': 'darkred', 'weight': 'normal', 'size': 16}
+        plt.xlabel('$x/c$', fontdict=font)
+        plt.ylabel('$C_p$', fontdict=font)
+        plt.gca().invert_yaxis()
+        plt.show()
+    tms['pos'].stop()
+
+    # display timers
+    print('CPU statistics')
+    print(tms)
+
+    # check results
+    tests = CTests()
+    tests.add(CTest('CL', sol.Cl, 0.5105, 1e-1))
+    tests.add(CTest('CD', sol.Cd, 0.0143, 1e-1))
+    tests.add(CTest('CM', sol.Cm, -0.1214, 1e-1))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()