.clang-format

@@ -4,7 +4,25 @@
 BasedOnStyle: LLVM
 UseTab: Never
 IndentWidth: 4
-BreakBeforeBraces: Allman
+BreakBeforeBraces: Custom
+BraceWrapping:
+  AfterNamespace: true
+  AfterClass: true
+  AfterStruct: true
+  AfterFunction: true
+  BeforeLambdaBody: false
+  AfterControlStatement: Always
+  BeforeElse: true
+  AfterCaseLabel: true
+  BeforeWhile: false
+  BeforeCatch: true
+  AfterEnum: true
+  AfterUnion: true
+  AfterExternBlock: true
+  SplitEmptyNamespace: false
+  SplitEmptyFunction: false
+  SplitEmptyRecord: false
+  IndentBraces: false
 AccessModifierOffset: -4
 SortIncludes: false
 ColumnLimit: 0

.gitlab-ci.yml

 default:
-    image: rboman/waves-py3:2020.3
+    image: rboman/waves-py3:2023.0
     before_script:
-        - source /opt/intel/mkl/bin/mklvars.sh intel64
-        - source /opt/intel/tbb/bin/tbbvars.sh intel64
+        - source /opt/intel/oneapi/mkl/latest/env/vars.sh
+        - source /opt/intel/oneapi/tbb/latest/env/vars.sh
         - echo $(nproc)
         - printenv | sort
 
@@ -24,8 +24,8 @@ format:
     <<: *global_tag_def
     stage: build
     script:
-        - clang-format --version # we use clang-format-10 exclusively
-        - ./scripts/format_code.py
+        - clang-format --version
+        - ./ext/amfe/scripts/format_code.py
         - mkdir -p patches
         - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi
     artifacts:

CMakeLists.txt

@@ -15,11 +15,7 @@
 # ----------------------------------------------------------------------------
 PROJECT(DARTFlo)
 # ----------------------------------------------------------------------------
-CMAKE_MINIMUM_REQUIRED(VERSION 3.1)
-IF(${CMAKE_VERSION} VERSION_GREATER "3.14.0") # we might want to update the project and use the NEW behavior here
-    cmake_policy(SET CMP0078 OLD)
-    cmake_policy(SET CMP0086 OLD)
-ENDIF()
+CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
 
 # -- I/O
 # Lib/Exe dir
@@ -60,22 +56,17 @@ ENDIF()
 
 # -- DEPENDENCIES
 # Python
-IF (CMAKE_VERSION VERSION_LESS 3.12.0)
-    FIND_PACKAGE(PythonInterp 3.6 REQUIRED)
-    FIND_PACKAGE(PythonLibs 3.6 REQUIRED)
-ELSE()
-    FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
-    # use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
-    SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
-    SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
-    SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS}) 
-    SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})     
-ENDIF()
+# use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
+FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
+SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
+SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
+SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS}) 
+SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})     
 
 # SWIG
 FIND_PACKAGE(SWIG REQUIRED)
 IF(CMAKE_GENERATOR MATCHES "Visual Studio") # not MSVC because of nmake & jom
-    SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}/${CMAKE_BUILD_TYPE}/")
+    SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}/$(Configuration)/")
 ELSE()
     SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}")
 ENDIF()
@@ -108,8 +99,13 @@ ENABLE_TESTING()
 # -- INSTALL
 IF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
     EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -m site --user-site OUTPUT_VARIABLE PY_SITE OUTPUT_STRIP_TRAILING_WHITESPACE)
+    STRING(REGEX REPLACE "\\\\" "/" PY_SITE ${PY_SITE})
     SET(CMAKE_INSTALL_PREFIX "${PY_SITE}/dartflo" CACHE STRING "Install location" FORCE)
     SET(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT FALSE)
+ELSE()
+    IF(NOT(CMAKE_INSTALL_PREFIX MATCHES "dartflo$"))
+        SET(CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}/dartflo" CACHE STRING "Install location" FORCE)
+    ENDIF()
 ENDIF()
 IF(UNIX AND NOT APPLE)
     SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}")

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

 # DARTFlo
-DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
-DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [openMDAO](https://openmdao.org/) so that aeroelastic computations and optimization can be easily carried out.
+DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/Python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
+DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [OpenMDAO](https://openmdao.org/) so that aeroelastic analysis and optimization calculations can be carried out easily.
 
 ![](/dox/title.png)
 
 ## Main features
-* Cross platform (Windows and Unix) C++/python code
-* Physical model
-  - unstructured meshes for 2D and 3D geometries using [gmsh](https://gmsh.info/)
+* Code
+  - cross platform (Windows and Unix)
+  - C++ with Python interface
+* Inputs and outputs
+  - 2D and 3D unstructured meshes in [Gmsh](https://gmsh.info/) format, interface with [GmshCFD](https://github.com/acrovato/gmshcfd/)
+  - results in Gmsh or [VTK](https://vtk.org/) format
+* Flow modeling
   - steady transonic flows
-  - viscous-inviscid interaction (2D only)
+  - viscous-inviscid interaction (see [BLASTER](https://gitlab.uliege.be/am-dept/blaster))
 * 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/)
-  - 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
-* Interfaced with
-  - [VTK](https://vtk.org/)
+  - adjoint solver
+  - mesh morphing
+* Linear algebra
+  - [Eigen](http://eigen.tuxfamily.org/)
+  - [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/) linear solvers
+* Multi-disciplinary anlaysis and optimization
   - [CUPyDO](https://github.com/ulgltas/CUPyDO)
-  - [openMDAO](https://openmdao.org/)
+  - [OpenMDAO](https://openmdao.org/)
 
 ## Documentation
-Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).
-
-## References
-Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design](http://hdl.handle.net/2268/251906), PhD thesis, University of Liège, 2020.
-
-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.
-
-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.
+Detailed instructions for building and using the code, as well as references to theory manuals and scientific articles can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).

dart/CMakeLists.txt

@@ -24,5 +24,4 @@ INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
               ${CMAKE_CURRENT_LIST_DIR}/utils.py
         DESTINATION dart)
 INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
-                  ${CMAKE_CURRENT_LIST_DIR}/viscous
         DESTINATION dart)

dart/README.md

@@ -2,6 +2,5 @@
 * [tests](dart/tests): simple tests to be run in battery checking the basic functionnalities of the solver
 * [validation](dart/validation): large-scale tests used to validate the solver results
 * [benchmark](dart/benchmark): large-scale tests used to monitor the performance
-* [viscous](dart/viscous): viscous-inviscid interaction module
 * [api](dart/api): APIs for DART
 * [cases](dart/cases): sample cases meant to be run using the APIs

dart/_src/CMakeLists.txt

@@ -19,34 +19,17 @@ INCLUDE(${SWIG_USE_FILE})
 FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
 FILE(GLOB ISRCS *.i)
 
-SET(CMAKE_SWIG_FLAGS "")
 SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
-
-SET(SWINCFLAGS 
--I${PROJECT_SOURCE_DIR}/dart/src
--I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/src
--I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
--I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/src
--I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
-)
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
-
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
-    SWIG_ADD_MODULE(dartw python ${ISRCS} ${SRCS})
-else()
-    SWIG_ADD_LIBRARY(dartw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_dartw)
-
-TARGET_INCLUDE_DIRECTORIES(_dartw PRIVATE ${PROJECT_SOURCE_DIR}/dart/_src
-                                          ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
-                                          ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
-                                          ${PYTHON_INCLUDE_PATH}
-)
-
-SWIG_LINK_LIBRARIES(dartw 
-                    dart tbox fwk ${PYTHON_LIBRARIES}
-)
-
-INSTALL(FILES ${CMAKE_SWIG_OUTDIR}/dartw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
-INSTALL(TARGETS _dartw DESTINATION ${CMAKE_INSTALL_PREFIX})
+SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_dartw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(dartw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET dartw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(dartw)
+
+TARGET_INCLUDE_DIRECTORIES(dartw PRIVATE ${PROJECT_SOURCE_DIR}/dart/_src
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
+                                         ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(dartw PRIVATE dart tbox fwk ${PYTHON_LIBRARIES})
+
+INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/dartw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS dartw DESTINATION ${CMAKE_INSTALL_PREFIX})

dart/_src/dartw.i

@@ -70,8 +70,8 @@ threads="1"
 %immutable dart::Body::Cm;
 %immutable dart::Body::nLoads;
 %include "wBody.h"
-%include "wBlowing.h"
 %warnfilter(+509);
+%include "wBlowing.h"
 
 %immutable dart::Problem::msh; // read-only variable (no setter)
 %immutable dart::Problem::nDim;

dart/api/core.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.
-
-
-## Initialize DART
-# Adrien Crovato
-
-def initDart(params, scenario='aerodynamic', task='analysis'):
-    """Initlialize DART for various API
-
-    Parameters
-    ----------
-    params : dict
-        Dictionary of parameters to configure DART
-    scenario : string, optional
-        Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
-    task : string, optional
-        Type of task (available: analysis, optimization) (default: analysis)
-
-    Returns
-    -------
-    _dim : int
-        Dimension of the problem
-    _qinf : float
-        Freestream dynamic pressure
-    _msh : tbox.MshData object
-        Mesh data structure
-    _wrtr : tbox.MshExport object
-        Mesh writer
-    _mrf : tbox.MshDeform object
-        Mesh morpher
-    _pbl : dart.Probem object
-        Problem definition
-    _bnd : Dart.Body object
-        Body of interest
-    _sol : dart.Newton object
-        Direct (Newton) solver
-    _adj : dart.Adjoint object
-        Adjoint solver
-    """
-    # Imports
-    import math
-    import fwk.wutils as wu
-    import tbox
-    import tbox.gmsh as gmsh
-    import dart
-
-    # Basic checks and config
-    # dimension
-    if params['Dim'] != 2 and params['Dim'] != 3:
-        raise RuntimeError('Problem dimension should be 2 or 3, but ' + params['Dim'] + ' was given!\n')
-    else:
-        _dim = params['Dim']
-    # aoa and aos
-    if 'AoA' in params:
-        alpha = params['AoA'] * math.pi / 180
-    else:
-        alpha = 0
-    if _dim == 2:
-        if 'AoS' in params and params['AoS'] != 0:
-            raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
-        else:
-            beta = 0
-        if 'Symmetry' in params:
-            raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
-    else:
-        if 'AoS' in params:
-            if params['AoS'] != 0 and 'Symmetry' in params:
-                raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-            beta = params['AoS'] * math.pi / 180
-        else:
-            beta = 0
-    # save format
-    if params['Format'] == 'vtk':
-        try:
-            import tboxVtk
-            Writer = tboxVtk.VtkExport
-            print('VTK libraries found! Results will be saved in VTK format.\n')
-        except:
-            Writer = tbox.GmshExport
-            print('VTK libraries not found! Results will be saved in gmsh format.\n')
-    else:
-        Writer = tbox.GmshExport
-        print('Results will be saved in gmsh format.\n')
-    # number of threads
-    if 'Threads' in params:
-        nthrd = params['Threads']
-    else:
-        nthrd = 1
-    wu.initMKL(nthrd) # initialize threading layer and number of threads
-    # verbosity
-    if 'Verb' in params:
-        verb = params['Verb']
-    else:
-        verb = 0
-    # scenario and task type
-    if scenario != 'aerodynamic' and scenario != 'aerostructural':
-        raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
-    if task != 'analysis' and task != 'optimization':
-        raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n')
-    # dynamic pressure
-    if scenario == 'aerostructural':
-        _qinf = params['Q_inf']
-    else:
-        _qinf = 0
-
-    # Mesh creation
-    _msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
-    # add the wake
-    if _dim == 2:
-        mshCrck = tbox.MshCrack(_msh, _dim)
-        mshCrck.setCrack(params['Wake'])
-        mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wing']])
-        mshCrck.run()
-    else:
-        for i in range(len(params['Wakes'])):
-            mshCrck = tbox.MshCrack(_msh, _dim)
-            mshCrck.setCrack(params['Wakes'][i])
-            mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wings'][i]])
-            if 'Fuselage' in params:
-                mshCrck.addBoundaries([params['Fuselage']])
-            if 'Symmetry' in params:
-                mshCrck.addBoundaries([params['Symmetry']])
-            if 'WakeTips' in params:
-                mshCrck.setExcluded(params['WakeTips'][i]) # 2.5D computations
-            mshCrck.run()
-    # save the updated mesh in native (gmsh) format and then set the writer to the requested format
-    tbox.GmshExport(_msh).save(_msh.name)
-    del mshCrck
-    _wrtr = Writer(_msh)
-    print('\n')
-
-    # Mesh morpher creation
-    if scenario == 'aerostructural' or task == 'optimization':
-        _mrf = tbox.MshDeform(_msh, _dim, nthrds=nthrd)
-        _mrf.setField(params['Fluid'])
-        _mrf.addFixed(params['Farfield'])
-        if _dim == 2:
-            _mrf.addMoving([params['Wing']])
-            _mrf.addInternal([params['Wake'], params['Wake']+'_'])
-        else:
-            if 'Fuselage' in params:
-                _mrf.addFixed(params['Fuselage'])
-            _mrf.setSymmetry(params['Symmetry'], 1)
-            for i in range(len(params['Wings'])):
-                if i == 0:
-                    _mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
-                else:
-                    _mrf.addFixed([params['Wings'][i]])
-                _mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
-        _mrf.initialize()
-    else:
-        _mrf = None
-
-    # Problem creation
-    _pbl = dart.Problem(_msh, _dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
-    # add the field
-    _pbl.set(dart.Fluid(_msh, params['Fluid'], params['M_inf'], _dim, alpha, beta))
-    # add the initial condition
-    _pbl.set(dart.Initial(_msh, params['Fluid'], _dim, alpha, beta))
-    # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
-    for i in range (len(params['Farfield'])):
-        _pbl.add(dart.Freestream(_msh, params['Farfield'][i], _dim, alpha, beta))
-    # add solid boundaries
-    if _dim == 2:
-        bnd = dart.Body(_msh, [params['Wing'], params['Fluid']])
-        _bnd = bnd
-        _pbl.add(bnd)
-    else:
-        _bnd = None
-        for bd in params['Wings']:
-            bnd = dart.Body(_msh, [bd, params['Fluid']])
-            if _bnd is None:
-                _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
-            _pbl.add(bnd)
-        if 'Fuselage' in params:
-            bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']])
-            _pbl.add(bnd)
-    # add Wake/Kutta boundary conditions
-    if _dim == 2:
-        _pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
-        _pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
-    else:
-        for i in range(len(params['Wakes'])):
-            _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['TeTips'][i]]))
-
-    # Direct (forward) solver creation
-    # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
-    # initialize the linear (inner) solver
-    if params['LSolver'] == 'PARDISO':
-        linsol = tbox.Pardiso()
-    elif params['LSolver'] == 'MUMPS':
-        linsol = tbox.Mumps()
-    elif params['LSolver'] == 'GMRES':
-        linsol = tbox.Gmres()
-        linsol.setFillFactor(params['G_fill']) if 'G_fill' in params else linsol.setFillFactor(2)
-        linsol.setRestart(params['G_restart']) if 'G_restart' in params else linsol.setRestart(50)
-        linsol.setTolerance(params['G_tol']) if 'G_tol' in params else linsol.setTolerance(1e-5)
-    else:
-        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + params['LSolver'] + ' was given!\n')
-    # initialize the nonlinear (outer) solver
-    _sol = dart.Newton(_pbl, linsol, rTol=params['Rel_tol'], aTol=params['Abs_tol'], mIt=params['Max_it'], nthrds=nthrd, vrb=verb)
-
-    # Adjoint (reverse) solver creation
-    if task == 'optimization':
-        if params['LSolver'] == 'PARDISO':
-            alinsol = tbox.Pardiso()
-        elif params['LSolver'] == 'MUMPS':
-            alinsol = tbox.Mumps()
-        else:
-            alinsol = tbox.Gmres()
-            alinsol.setFillFactor(params['G_fill']) if 'G_fill' in params else alinsol.setFillFactor(2)
-            alinsol.setRestart(params['G_restart']) if 'G_restart' in params else alinsol.setRestart(50)
-            alinsol.setTolerance(1e-12)
-        _adj = dart.Adjoint(_sol, _mrf, alinsol)
-    else:
-        _adj = None
-
-    # Return
-    return _dim, _qinf, _msh, _wrtr, _mrf, _pbl, _bnd, _sol, _adj
+#!/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.
+
+
+## Initialize DART
+# Adrien Crovato
+
+def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
+    """Initlialize DART for various API
+
+    Parameters
+    ----------
+    cfg : dict
+        Dictionary of parameters to configure DART
+    scenario : string, optional
+        Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
+    task : string, optional
+        Type of task (available: analysis, optimization) (default: analysis)
+    viscous : bool, optional
+        Whether to configure DART for viscous-inviscid interaction (default: False)
+
+    Returns
+    -------
+    A dictionary containing
+        dim : int
+            Dimension of the problem
+        qinf : float
+            Freestream dynamic pressure
+        msh : tbox.MshData object
+            Mesh data structure
+        wrt : tbox.MshExport object
+            Mesh writer
+        mrf : tbox.MshDeform object
+            Mesh morpher
+        pbl : dart.Probem object
+            Problem definition
+        bnd : Dart.Body object
+            Body of interest
+        blwb : Dart.Blowing object
+            Transpiration B.C. on body
+        blww : Dart.Blowing object
+            Transpiration B.C. on wake
+        sol : dart.Newton object
+            Direct (Newton) solver
+        adj : dart.Adjoint object
+            Adjoint solver
+    """
+    # Imports
+    import dart
+    import tbox
+    import tbox.gmsh as gmsh
+    import math
+
+    # Basic checks and config
+    # dimension
+    if cfg['Dim'] != 2 and cfg['Dim'] != 3:
+        raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n')
+    else:
+        _dim = cfg['Dim']
+    # angles of attack and side slip, and freestream Mach number
+    alpha = cfg.get('AoA', 0.) * math.pi / 180
+    beta = cfg.get('AoS', 0.) * math.pi / 180
+    minf = cfg.get('M_inf', 0.)
+    if _dim == 2:
+        if beta != 0:
+            raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
+        if 'Symmetry' in cfg:
+            raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
+    else:
+        if beta != 0 and 'Symmetry' in cfg:
+            raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
+    # save format
+    if cfg.get('Format', 'vtk') == 'vtk':
+        try:
+            import tboxVtk
+            Writer = tboxVtk.VtkExport
+            print('VTK libraries found! Results will be saved in VTK format.\n')
+        except:
+            Writer = tbox.GmshExport
+            print('VTK libraries not found! Results will be saved in gmsh format.\n')
+    else:
+        Writer = tbox.GmshExport
+        print('Results will be saved in gmsh format.\n')
+    # number of threads and verbosity
+    nthrd = cfg.get('Threads', 1)
+    verb = cfg.get('Verb', 1)
+    # scenario and task type
+    if scenario != 'aerodynamic' and scenario != 'aerostructural':
+        raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
+    if task != 'analysis' and task != 'optimization':
+        raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n')
+    # dynamic pressure
+    if scenario == 'aerostructural':
+        _qinf = cfg['Q_inf']
+    else:
+        _qinf = None
+
+    # Mesh creation
+    pars = cfg.get('Pars', {})
+    _msh = gmsh.MeshLoader(cfg['File']).execute(**pars)
+    # add the wake
+    if _dim == 2:
+        mshCrck = tbox.MshCrack(_msh, _dim)
+        mshCrck.setCrack(cfg['Wake'])
+        mshCrck.addBoundary(cfg['Fluid'])
+        mshCrck.addBoundary(cfg['Farfield'][-1])
+        mshCrck.addBoundary(cfg['Wing'])
+        mshCrck.run()
+    else:
+        for i in range(len(cfg['Wakes'])):
+            mshCrck = tbox.MshCrack(_msh, _dim)
+            mshCrck.setCrack(cfg['Wakes'][i])
+            mshCrck.addBoundary(cfg['Fluid'])
+            mshCrck.addBoundary(cfg['Farfield'][-1])
+            mshCrck.addBoundary(cfg['Wings'][i])
+            if 'Fuselage' in cfg:
+                mshCrck.addBoundary(cfg['Fuselage'])
+            if 'Symmetry' in cfg:
+                mshCrck.addBoundary(cfg['Symmetry'])
+            if 'WakeTips' in cfg:
+                mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
+            mshCrck.run()
+    # save the updated mesh in native (gmsh) format and then set the writer to the requested format
+    tbox.GmshExport(_msh).save(_msh.name)
+    del mshCrck
+    _wrt = Writer(_msh)
+    print('\n')
+
+    # Mesh morpher creation
+    if scenario == 'aerostructural' or task == 'optimization':
+        _mrf = tbox.MshDeform(_msh, tbox.Gmres(2, 1e-6, 50, 1e-8, 500), _dim, nthrds=nthrd, vrb=verb)
+        _mrf.setField(cfg['Fluid'])
+        for tag in cfg['Farfield']:
+            _mrf.addFixed(tag)
+        if _dim == 2:
+            _mrf.addMoving(cfg['Wing'])
+            _mrf.addInternal(cfg['Wake'], cfg['Wake']+'_')
+        else:
+            if 'Fuselage' in cfg:
+                _mrf.addFixed(cfg['Fuselage'])
+            if 'Symmetry' in cfg:
+                _mrf.setSymmetry(cfg['Symmetry'], 1)
+            for i in range(len(cfg['Wings'])):
+                if i == 0:
+                    _mrf.addMoving(cfg['Wings'][i]) # TODO body of interest (FSI/OPT) is always first given body
+                else:
+                    _mrf.addFixed(cfg['Wings'][i])
+                _mrf.addInternal(cfg['Wakes'][i], cfg['Wakes'][i]+'_')
+        _mrf.initialize()
+    else:
+        _mrf = None
+
+    # Problem creation
+    _pbl = dart.Problem(_msh, _dim, alpha, beta, minf, cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'])
+    # add the field
+    _pbl.set(dart.Fluid(_msh, cfg['Fluid'], minf, _dim, alpha, beta))
+    # add the initial condition
+    _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta))
+    # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
+    for i in range (len(cfg['Farfield'])):
+        _pbl.add(dart.Freestream(_msh, cfg['Farfield'][i], _dim, alpha, beta))
+    # add solid boundaries
+    if _dim == 2:
+        bnd = dart.Body(_msh, [cfg['Wing'], cfg['Fluid']])
+        _bnd = bnd
+        _pbl.add(bnd)
+    else:
+        _bnd = None
+        for bd in cfg['Wings']:
+            bnd = dart.Body(_msh, [bd, cfg['Fluid']])
+            if _bnd is None:
+                _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
+            _pbl.add(bnd)
+        if 'Fuselage' in cfg:
+            bnd = dart.Body(_msh, [cfg['Fuselage'], cfg['Fluid']])
+            _pbl.add(bnd)
+    # add Wake/Kutta boundary conditions
+    # TODO refactor Wake "exclusion" for 3D?
+    if _dim == 2:
+        _pbl.add(dart.Wake(_msh, [cfg['Wake'], cfg['Wake']+'_', cfg['Fluid']]))
+        _pbl.add(dart.Kutta(_msh, [cfg['Te'], cfg['Wake']+'_', cfg['Wing'], cfg['Fluid']]))
+    else:
+        for i in range(len(cfg['Wakes'])):
+            if 'WakeExs' in cfg:
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeExs'][i]])) # 3D + fuselage
+            elif 'WakeTips' in cfg:
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeTips'][i]])) # 3D
+            else:
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid']])) # 2.5D
+            _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']]))
+    # add transpiration (blowing) boundary conditions
+    if viscous:
+        _blwb = dart.Blowing(_msh, cfg['Wing'])
+        _blww = dart.Blowing(_msh, cfg['Wake'])
+        _pbl.add(_blwb)
+        _pbl.add(_blww)
+    else:
+        _blwb = None
+        _blww = None
+
+    # Direct (forward) solver creation
+    # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
+    # initialize the linear (inner) solver
+    if cfg['LSolver'] == 'PARDISO':
+        linsol = tbox.Pardiso()
+    elif cfg['LSolver'] == 'MUMPS':
+        linsol = tbox.Mumps()
+    elif cfg['LSolver'] == 'SparseLU':
+        linsol = tbox.SparseLu()
+    elif cfg['LSolver'] == 'GMRES':
+        linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 200)
+    else:
+        raise RuntimeError('Available linear solvers: PARDISO, MUMPS, SparseLU or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
+    # initialize the nonlinear (outer) solver
+    _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)
+
+    # Adjoint (reverse) solver creation
+    if task == 'optimization':
+        _adj = dart.Adjoint(_sol, _mrf)
+    else:
+        _adj = None
+
+    # Return
+    return {
+        'dim': _dim,
+        'qinf': _qinf,
+        'msh': _msh,
+        'wrt': _wrt,
+        'mrf': _mrf,
+        'pbl': _pbl,
+        'bnd': _bnd,
+        'blwb': _blwb,
+        'blww': _blww,
+        'sol': _sol,
+        'adj': _adj
+        }

dart/api/cupydo.py

+# -*- coding: utf-8 -*-
+
+# Copyright 2021 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.
+
+## CUPYDO interface
+# Adrien Crovato
+
+from cupydo.genericSolvers import FluidSolver
+import numpy as np
+import sys
+
+class Dart(FluidSolver):
+    """DART interface for CUPyDO
+    """
+    def __init__(self, p, _nthreads):
+        # load the python module and initialize the solver
+        module = __import__(p['cfdFile'])
+        params = module.getParams()
+        params['Threads'] = _nthreads
+        from .core import init_dart
+        _dart = init_dart(params, scenario='aerostructural', task='analysis')
+        self.qinf, self.msh, self.writer, self.morpher, self.boundary, self.solver = (_dart.get(key) for key in ['qinf', 'msh', 'wrt', 'mrf', 'bnd', 'sol'])
+
+        # count fsi nodes and get their positions
+        self.nNodes = self.boundary.nodes.size()
+        self.nHaloNode = 0
+        self.nPhysicalNodes = self.nNodes - self.nHaloNode
+        self.nodalInitPosX, self.nodalInitPosY, self.nodalInitPosZ = self.getNodalInitialPositions()
+
+        # init save frequency (fsi)
+        self.saveFreq = params.get('SaveFreq', sys.maxsize)
+
+        # generic init
+        FluidSolver.__init__(self, p)
+
+    def run(self, t1, t2):
+        """Run the solver for one steady (time) iteration.
+        """
+        status = self.solver.run()
+        if status > 1:
+            return False
+        self.__setCurrentState()
+        return True
+
+    def __setCurrentState(self):
+        """Compute nodal forces from nodal normalized forces
+        """
+        i = 0
+        for n in self.boundary.nodes:
+            self.nodalLoad_X[i] = self.qinf * self.boundary.nLoads[i][0]
+            self.nodalLoad_Y[i] = self.qinf * self.boundary.nLoads[i][1]
+            self.nodalLoad_Z[i] = self.qinf * self.boundary.nLoads[i][2]
+            i += 1
+
+    def getNodalInitialPositions(self):
+        """Get the initial position of each node
+        """
+        x0 = np.zeros(self.nPhysicalNodes)
+        y0 = np.zeros(self.nPhysicalNodes)
+        z0 = np.zeros(self.nPhysicalNodes)
+        for i in range(self.boundary.nodes.size()):
+            n = self.boundary.nodes[i]
+            x0[i] = n.pos[0]
+            y0[i] = n.pos[1]
+            z0[i] = n.pos[2]
+
+        return (x0, y0, z0)
+
+    def getNodalIndex(self, iVertex):
+        """Get index of each node
+        """
+        no = self.boundary.nodes[iVertex].no
+        return no
+
+    def applyNodalDisplacements(self, dx, dy, dz, dt, haloNodesDisplacements):
+        """Apply displacements coming from solid solver to f/s interface after saving
+        """
+        self.morpher.savePos()
+        for i in range(self.boundary.nodes.size()):
+            self.boundary.nodes[i].pos[0] = self.nodalInitPosX[i] + dx[i]
+            self.boundary.nodes[i].pos[1] = self.nodalInitPosY[i] + dy[i]
+            self.boundary.nodes[i].pos[2] = self.nodalInitPosZ[i] + dz[i]
+
+    def meshUpdate(self, nt):
+        """Deform the mesh using linear elasticity equations
+        """
+        self.morpher.deform()
+
+    def save(self, nt):
+        """Save data on disk at each converged timestep
+        """
+        self.solver.save(self.writer, '_converged')
+        self.writer.save(self.msh.name + '_converged')
+
+    def initRealTimeData(self):
+        """Initialize history file
+        """
+        histFile = open('DartHistory.dat', 'w')
+        histFile.write('{0:>12s}   {1:>12s}   {2:>12s}   {3:>12s}   {4:>12s}\n'.format('Time', 'FSI_Iter', 'C_Lift', 'C_Drag', 'C_Moment'))
+        histFile.close()
+
+    def saveRealTimeData(self, time, nFSIIter):
+        """Save data at each fsi iteration
+        """
+        # history at each iteration
+        histFile = open('DartHistory.dat', 'a')
+        histFile.write('{0:12.6f}   {1:12d}   {2:12.6f}   {3:12.6f}   {4:12.6f}\n'.format(time, nFSIIter, self.solver.Cl, self.solver.Cd, self.solver.Cm))
+        histFile.close()
+        # full solution at user-defined frequency
+        if np.mod(nFSIIter+1, self.saveFreq) == 0:
+            self.solver.save(self.writer, '_{:04d}'.format(int(nFSIIter+1)//int(self.saveFreq)))
+
+    def printRealTimeData(self, time, nFSIIter):
+        """Print data on screen at the end of fsi simulation
+        """
+        print('[DART lift, drag, moment]: {0:6.3f}, {1:6.4f}, {2:6.3f}'.format(self.solver.Cl, self.solver.Cd, self.solver.Cm))
+        print('')
+
+    def exit(self):
+        """Clear memory
+        """
+        del self.boundary
+        del self.solver
+        del self.morpher
+        del self.writer
+        del self.msh
+        
\ No newline at end of file

dart/api/internal/polar.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.
-
-
-## Polar
-# Adrien Crovato
-#
-# Compute aerodynamic load coefficients as a function of angle of attack
-
-import math
-import dart.utils as dU
-import tbox.utils as tU
-import fwk
-from fwk.coloring import ccolors
-from ..core import initDart
-
-class Polar:
-    """Utility to compute the polar of a lifting surface
-
-    Attributes
-    ----------
-    tms : fwk.Timers object
-        dictionary of timers
-    dim : int
-        problem dimensions
-    format : str
-        output format type
-    slice : array of float
-        arrays of y-coordinates to perform slices (for 3D only)
-    tag : int
-        index of physical tag to slice (see gmsh)
-    mach : float
-        freestream Mach number
-    alphas : array of float
-        range of angles of attack to compute
-    Cls : array of float
-        lift coefficients for the different AoAs
-    Cds : array of float
-        drag coefficients for the different AoAs
-    Cms : array of float
-        moment coefficients for the different AoAs
-    """
-    def __init__(self, p):
-        """Setup and configure
-
-        Parameters
-        ----------
-        p : dict
-            dictionary of parameters to configure DART
-        """
-        # basic init
-        self.tms = fwk.Timers()
-        self.tms["total"].start()
-        self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
-        # save some parameters for later use
-        self.format = p['Format']
-        try:
-            self.slice = p['Slice']
-            self.tag = p['TagId']
-        except:
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
-        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
-
-    def run(self):
-        """Compute the polar
-        """
-        # initialize the loop
-        ac = 0 if self.alphas[0] == self.alphas[-1] else 1
-        self.Cls = []
-        self.Cds = []
-        self.Cms = []
-        print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
-        for alpha in self.alphas:
-            # define current angle of attack
-            alpha = alpha*math.pi/180
-            # update problem
-            self.pbl.update(alpha)
-            # run the solver and save the results
-            print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
-            self.tms["solver"].start()
-            self.sol.run()
-            self.tms["solver"].stop()
-            self.sol.save(self.wrtr, ac)
-            # extract Cp
-            if self.dim == 2:
-                Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-                tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
-            elif self.dim == 3 and self.format == 'vtk' and self.slice:
-                dU.writeSlices(self.msh.name, self.slice, self.tag, ac)
-            # extract force coefficients
-            self.Cls.append(self.sol.Cl)
-            self.Cds.append(self.sol.Cd)
-            self.Cms.append(self.sol.Cm)
-            # end loop
-            ac+=1
-
-    def disp(self):
-        """Display the results and draw the polar
-        """
-        # display results
-        print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
-        print("       M    alpha       Cl       Cd       Cm")
-        i = 0
-        while i < len(self.alphas):
-            print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i]))
-            i = i+1
-        print('\n')
-        # display timers
-        self.tms["total"].stop()
-        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
-        print(self.tms)
-        # plot results
-        if self.alphas[0] != self.alphas[-1]:
-            tU.plot(self.alphas, self.Cls, "alpha", "Cl", "")
-            tU.plot(self.alphas, self.Cds, "alpha", "Cd", "")
-            tU.plot(self.alphas, self.Cms, "alpha", "Cm", "")
+#!/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.
+
+
+## Polar
+# Adrien Crovato
+#
+# Compute aerodynamic load coefficients as a function of angle of attack
+
+import math
+import dart.utils as dU
+import tbox.utils as tU
+import fwk
+from fwk.coloring import ccolors
+from ..core import init_dart
+
+class Polar:
+    """Utility to compute the polar of a lifting surface
+
+    Attributes
+    ----------
+    tms : fwk.Timers object
+        dictionary of timers
+    dim : int
+        problem dimensions
+    format : str
+        output format type
+    slice : array of float
+        arrays of y-coordinates to perform slices (for 3D only)
+    tag : int
+        index of physical tag to slice (see gmsh)
+    mach : float
+        freestream Mach number
+    alphas : array of float
+        range of angles of attack to compute
+    Cls : array of float
+        lift coefficients for the different AoAs
+    Cds : array of float
+        drag coefficients for the different AoAs
+    Cms : array of float
+        moment coefficients for the different AoAs
+    """
+    def __init__(self, p):
+        """Setup and configure
+
+        Parameters
+        ----------
+        p : dict
+            dictionary of parameters to configure DART
+        """
+        # basic init
+        self.tms = fwk.Timers()
+        self.tms["total"].start()
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
+        # save some parameters for later use
+        self.format = p.get('Format', 'vtk')
+        self.slice = p.get('Slice')
+        self.tag = p.get('TagId')
+        self.mach = self.pbl.M_inf
+        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
+
+    def run(self):
+        """Compute the polar
+        """
+        # initialize the loop
+        self.Cls = []
+        self.Cds = []
+        self.Cms = []
+        print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
+        for i in range(len(self.alphas)):
+            # define current angle of attack
+            alpha = self.alphas[i]*math.pi/180
+            acs = '_{:04d}'.format(i)
+            # update problem and reset ICs to improve robustness in transonic cases
+            self.pbl.update(alpha)
+            self.sol.reset()
+            # run the solver and save the results
+            print(ccolors.ANSI_BLUE + 'Running the solver for', self.alphas[i], 'degrees', ccolors.ANSI_RESET)
+            self.tms["solver"].start()
+            self.sol.run()
+            self.tms["solver"].stop()
+            self.sol.save(self.wrt, acs)
+            # extract Cp
+            if self.dim == 2:
+                Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
+                tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.4e', ',', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, cp', '')
+            elif self.dim == 3 and self.format == 'vtk' and self.slice:
+                dU.write_slices(self.msh.name, self.slice, self.tag, acs)
+            # extract force coefficients
+            self.Cls.append(self.sol.Cl)
+            self.Cds.append(self.sol.Cd)
+            self.Cms.append(self.sol.Cm)
+
+    def disp(self):
+        """Display the results and draw the polar
+        """
+        # display results
+        print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
+        print("       M    alpha       Cl       Cd       Cm")
+        i = 0
+        while i < len(self.alphas):
+            print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i]))
+            i = i+1
+        print('\n')
+        # display timers
+        self.tms["total"].stop()
+        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
+        print(self.tms)
+        # plot results
+        if self.alphas[0] != self.alphas[-1]:
+            tU.plot(self.alphas, self.Cls, {'xlabel': 'alpha', 'ylabel': 'Cl'})
+            tU.plot(self.alphas, self.Cds, {'xlabel': 'alpha', 'ylabel': 'Cd'})
+            tU.plot(self.alphas, self.Cms, {'xlabel': 'alpha', 'ylabel': 'Cm'})

dart/api/internal/trim.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.
-
-
-## Trim analysis
-# Adrien Crovato
-#
-# Find the angle of attack to match a specified lift coefficient
-
-import math
-import dart.utils as dU
-import tbox.utils as tU
-import fwk
-from fwk.coloring import ccolors
-from ..core import initDart
-
-class Trim:
-    """Utility to perform the trim analysis of a given lifting configuration
-    Find the angle of attack to reach a desired lift coefficient
-
-    Attributes
-    ----------
-    tms : fwk.Timers object
-        dictionary of timers
-    dim : int
-        problem dimensions
-    format : str
-        output format type
-    slice : array of float
-        arrays of y-coordinates to perform slices (for 3D only)
-    tag : int
-        index of physical tag to slice (see gmsh)
-    mach : float
-        freestream Mach number
-    clTgt : float
-        target lift coefficient
-    alpha : float
-        initial angle of attack
-    dCl : float
-        initial (estimated) slope of lift curve
-    """
-    def __init__(self, p):
-        """Setup and configure
-
-        Parameters
-        ----------
-        p : dict
-            dictionary of parameters to configure DART
-        """
-        # basic init
-        self.tms = fwk.Timers()
-        self.tms["total"].start()
-        self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
-        # save some parameters for later use
-        self.format = p['Format']
-        try:
-            self.slice = p['Slice']
-            self.tag = p['TagId']
-        except:
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
-        self.clTgt = p['CL']
-        self.alpha = p['AoA']*math.pi/180
-        self.dCl = p['dCL']
-
-    def run(self):
-        """Perform the trim analysis
-        """
-        # initialize loop
-        it = 0
-        while True:
-            # define current angle of attack
-            print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
-            # update problem
-            self.pbl.update(self.alpha)
-            # run the solver
-            self.tms["solver"].start()
-            success = self.sol.run()
-            if not success:
-                raise RuntimeError('Flow solver diverged!\n')
-            self.tms["solver"].stop()
-            # update slope
-            if it != 0:
-                self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_)
-            # break or store values and update AoA
-            if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl):
-                break
-            else:
-                Cl_ = self.sol.Cl
-                alpha_ = self.alpha
-                dAlpha = (self.clTgt - self.sol.Cl) / self.dCl
-                self.alpha += dAlpha
-                it += 1
-
-        # save results
-        self.sol.save(self.wrtr)
-        # extract Cp
-        if self.dim == 2:
-            Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
-            tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "")
-        elif self.dim == 3 and self.format == 'vtk' and self.slice:
-            dU.writeSlices(self.msh.name, self.slice, self.tag)
-
-    def disp(self):
-        """Display the results
-        """
-        # display results
-        print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
-        print("       M    alpha      dCl       Cl       Cd       Cm")
-        print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm))
-        print('\n')
-        # display timers
-        self.tms["total"].stop()
-        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
-        print(self.tms)
+#!/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.
+
+
+## Trim analysis
+# Adrien Crovato
+#
+# Find the angle of attack to match a specified lift coefficient
+
+import math
+import dart.utils as dU
+import tbox.utils as tU
+import fwk
+from fwk.coloring import ccolors
+from ..core import init_dart
+
+class Trim:
+    """Utility to perform the trim analysis of a given lifting configuration
+    Find the angle of attack to reach a desired lift coefficient
+
+    Attributes
+    ----------
+    tms : fwk.Timers object
+        dictionary of timers
+    dim : int
+        problem dimensions
+    format : str
+        output format type
+    slice : array of float
+        arrays of y-coordinates to perform slices (for 3D only)
+    tag : int
+        index of physical tag to slice (see gmsh)
+    mach : float
+        freestream Mach number
+    clTgt : float
+        target lift coefficient
+    alpha : float
+        current angle of attack
+    dCl : float
+        current slope of lift curve
+    """
+    def __init__(self, p):
+        """Setup and configure
+
+        Parameters
+        ----------
+        p : dict
+            dictionary of parameters to configure DART
+        """
+        # basic init
+        self.tms = fwk.Timers()
+        self.tms["total"].start()
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
+        # save some parameters for later use
+        self.format = p.get('Format', 'vtk')
+        self.slice = p.get('Slice')
+        self.tag = p.get('TagId')
+        self.mach = self.pbl.M_inf
+        self.alpha = self.pbl.alpha
+        self.clTgt = p['CL']
+        self.dCl = p['dCL']
+
+    def run(self):
+        """Perform the trim analysis
+        """
+        # initialize loop
+        it = 0
+        while True:
+            # define current angle of attack
+            print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
+            # update problem and reset ICs to improve robustness in transonic cases
+            self.pbl.update(self.alpha)
+            self.sol.reset()
+            # run the solver
+            self.tms["solver"].start()
+            status = self.sol.run()
+            if status >= 2:
+                raise RuntimeError('Flow solver diverged!\n')
+            self.tms["solver"].stop()
+            # update slope
+            if it != 0:
+                self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_)
+            # break or store values and update AoA
+            if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl):
+                break
+            else:
+                Cl_ = self.sol.Cl
+                alpha_ = self.alpha
+                dAlpha = (self.clTgt - self.sol.Cl) / self.dCl
+                self.alpha += dAlpha
+                it += 1
+
+        # save results
+        self.sol.save(self.wrt)
+        # extract Cp
+        if self.dim == 2:
+            Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
+            tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.4e', ',', 'alpha = '+str(self.alpha*180/math.pi)+' deg\nx, y, cp', '')
+        elif self.dim == 3 and self.format == 'vtk' and self.slice:
+            dU.write_slices(self.msh.name, self.slice, self.tag)
+
+    def disp(self):
+        """Display the results
+        """
+        # display results
+        print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
+        print("       M    alpha      dCl       Cl       Cd       Cm")
+        print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm))
+        print('\n')
+        # display timers
+        self.tms["total"].stop()
+        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
+        print(self.tms)

dart/benchmark/onera.py

@@ -19,13 +19,13 @@
 ## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84 for benchmark
 # Adrien Crovato
 
-import numpy as np
 import dart.default as floD
 import dart.utils as floU
 import tbox.utils as tbxU
 import fwk
 from fwk.testing import *
 from fwk.coloring import ccolors
+import numpy as np
 
 def newton(pbl):
     from fwk.wutils import parseargs
@@ -34,14 +34,10 @@ def newton(pbl):
     # Pardiso and GMRES should give similar perfs
     k = parseargs().k
     try:
-        newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=2)
+        newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=3)
     except:
-        gmres = tbox.Gmres()
-        gmres.setFillFactor(2)
-        gmres.setDropTol(1e-6)
-        gmres.setRestart(50)
-        gmres.setTolerance(1e-5)
-        newton = dart.Newton(pbl, gmres, nthrds=k, vrb=2)
+        gmres = tbox.Gmres(2, 1e-6, 50, 1e-5)
+        newton = dart.Newton(pbl, gmres, nthrds=k, vrb=3)
     return newton
 
 def main():
@@ -66,7 +62,7 @@ def main():
 
     # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
-    pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
+    pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
     tms['pre'].stop()
 
     # solve problem

dart/cases/coyote.py

@@ -30,7 +30,7 @@ def getParam():
     # Arguments
     args = parseargs()
     p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
     # Specific
     scl = 2 # scaling factor for lifting surface mesh size
     wrtems = scl*0.036 # wing root trailing edge mesh size
@@ -44,7 +44,7 @@ def getParam():
     p['Dim'] = 3 # Problem dimension
     p['Format'] = 'vtk' # Save format (vtk or gmsh)
     p['Slice'] = [1.8, 3.5, 5.3] # array of (y) coordinates to perform slice along the span (empty if none)
-    p['TagId'] = 9 # tag number of physical group to be sliced
+    p['TagId'] = 11 # tag number of physical group to be sliced
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
@@ -52,7 +52,8 @@ def getParam():
     p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
     p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
     p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
-    p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
+    p['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.)
+    p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge
     # Freestream
     p['M_inf'] = 0.8 # Freestream Mach number
     p['AoA_begin'] = 1 # Freestream angle of attack (begin)

dart/cases/n0012.py

@@ -26,7 +26,7 @@ def getParam():
     # Arguments
     args = parseargs()
     p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
     # Input/Output
     p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012.geo') # Input file containing the model
     p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model

dart/cases/n0012_3.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.
-
-
-## NACA 0012 wing configuration file for flow polar script
-# Adrien Crovato
-
-def getParam():
-    import os.path
-    from fwk.wutils import parseargs
-    p = {}
-    # Specific
-    span = 1. # span length
-    # Arguments
-    args = parseargs()
-    p['Threads'] = args.k
-    p['Verb'] = args.verb
-    # Input/Output
-    p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012_3.geo') # Input file containing the model
-    p['Pars'] = {'spn' : span, 'lgt' : 3., 'wdt' : 3., 'hgt' : 3., 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} # Parameters for input file model
-    p['Dim'] = 3 # Problem dimension
-    p['Format'] = 'vtk' # Save format (vtk or gmsh)
-    p['Slice'] = [0.01, 0.25, 0.5, 0.75, 0.95] # array of (y) coordinates to perform slice along the span (empty if none)
-    p['TagId'] = 4 # tag number of physical group to be sliced
-    # Markers
-    p['Fluid'] = 'field' # Name of physical group containing the fluid
-    p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
-    p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
-    p['Wings'] = ['wing'] # LIST of names of physical groups containing the lifting surface body
-    p['Wakes'] = ['wake'] # LIST of names of physical group containing the wake
-    p['WakeTips'] = ['wakeTip'] # LIST of names of physical group containing the edge of the wake
-    p['TeTips'] = ['teTip'] # LIST of names of physical group containing the edge of the wake and the trailing edge
-    # Freestream
-    p['M_inf'] = 0.3 # Freestream Mach number
-    p['AoA_begin'] = -3. # Freestream angle of attack (begin)
-    p['AoA_end'] = 3. # Freestream angle of attack (end)
-    p['AoA_step'] = 1. # Freestream angle of attack (step)
-    # Geometry
-    p['S_ref'] = 1.*span # Reference surface length
-    p['c_ref'] = 1. # Reference chord length
-    p['x_ref'] = 0.25 # Reference point for moment computation (x)
-    p['y_ref'] = 0. # Reference point for moment computation (y)
-    p['z_ref'] = 0. # Reference point for moment computation (z)
-    # Numerical
-    p['LSolver'] = 'MUMPS' # Linear solver (PARDISO, GMRES, or MUMPS)
-    p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
-    p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
-    p['Max_it'] = 10 # Solver maximum number of iterations
-    return p
-
-def main():
-    # Script call
-    import dart.api.internal.polar as p
-    polar = p.Polar(getParam())
-    polar.run()
-    polar.disp()
-
-if __name__ == "__main__":
-    main()

dart/cases/n64A410.py

@@ -26,7 +26,7 @@ def getParam():
     # Arguments
     args = parseargs()
     p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
     # Input/Output
     p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n64A410.geo') # Input file containing the model
     p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model