diff --git a/README.md b/README.md
index fcca63f3dfcd228ea3479e5f4cb197eb45721d95..a5a8976e196f8525d9954547ac1456962afa6c22 100644
--- a/README.md
+++ b/README.md
@@ -9,7 +9,7 @@ DART is currently capable of rapidly solving steady transonic flows on arbitrary
* Physical model
- unstructured meshes for 2D and 3D geometries using [gmsh](https://gmsh.info/)
- steady transonic flows
- - viscous-inviscid interation (2D only)
+ - 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/)
- direct (forward) and adjoint (backward) modes
@@ -20,12 +20,12 @@ DART is currently capable of rapidly solving steady transonic flows on arbitrary
- [CUPyDO](https://github.com/ulgltas/CUPyDO)
- [openMDAO](https://openmdao.org/)
-## Build
-Detailed build instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).
+## 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 A., Boman R., Guner H., Terrapon V., Dimitriadis G., Almeida H., Prado A., Breviglieri C., Cabral P., and Silva, G., [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., 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.
diff --git a/dart/CMakeLists.txt b/dart/CMakeLists.txt
index 0ac4587bff42c5e93fce7ce74c10c5ead64ffe07..2d2803f7ecc66b2a8e3cf5a3ea602f8a7f61b5d6 100644
--- a/dart/CMakeLists.txt
+++ b/dart/CMakeLists.txt
@@ -23,7 +23,6 @@ MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
${CMAKE_CURRENT_LIST_DIR}/utils.py
DESTINATION dart)
-INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/interfaces
- ${CMAKE_CURRENT_LIST_DIR}/scripts
+INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
${CMAKE_CURRENT_LIST_DIR}/viscous
DESTINATION dart)
diff --git a/dart/README.md b/dart/README.md
index 62c0342d592267c3e39c94ee677f75c20b43196d..71a9c79ff3dab1a952da7c1cc5059a044104fce1 100644
--- a/dart/README.md
+++ b/dart/README.md
@@ -2,5 +2,6 @@
* [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
-* [cases](dart/cases) : sample configuration files meant to be run with scripts, user-defined case should be derived from these
-* [scripts](dart/scripts): various scripts for running the solver automatically
+* [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
diff --git a/dart/interfaces/__init__.py b/dart/api/__init__.py
similarity index 100%
rename from dart/interfaces/__init__.py
rename to dart/api/__init__.py
diff --git a/dart/scripts/__init__.py b/dart/api/internal/__init__.py
similarity index 100%
rename from dart/scripts/__init__.py
rename to dart/api/internal/__init__.py
diff --git a/dart/scripts/config.py b/dart/api/internal/config.py
similarity index 98%
rename from dart/scripts/config.py
rename to dart/api/internal/config.py
index 64ba2ab4299e56cc954b2159a1ee1790fc3588c8..7c40318b66d7c2b55419632469bb8f8edf8d713b 100644
--- a/dart/scripts/config.py
+++ b/dart/api/internal/config.py
@@ -90,7 +90,7 @@ class Config:
# mesh the airfoil
print(ccolors.ANSI_BLUE + 'Meshing the airfoil...' + ccolors.ANSI_RESET)
- self.msh = gmsh.MeshLoader(p['File'],__file__).execute(**p['Pars'])
+ self.msh = gmsh.MeshLoader(p['File']).execute(**p['Pars'])
if p['Dim'] == 2:
mshCrck = tbox.MshCrack(self.msh, p['Dim'])
mshCrck.setCrack(p['Wake'])
diff --git a/dart/scripts/polar.py b/dart/api/internal/polar.py
similarity index 99%
rename from dart/scripts/polar.py
rename to dart/api/internal/polar.py
index f0f0cf391f7500d9c646a61cd17d9a9597ce47d1..75644e5cde0b313276aecb63c148544b7fc0da47 100644
--- a/dart/scripts/polar.py
+++ b/dart/api/internal/polar.py
@@ -25,7 +25,7 @@ import math
import dart.utils as dU
import tbox.utils as tU
from fwk.coloring import ccolors
-from dart.scripts.config import Config
+from .config import Config
class Polar(Config):
"""Utility to compute the polar of a lifting surface
diff --git a/dart/scripts/trim.py b/dart/api/internal/trim.py
similarity index 99%
rename from dart/scripts/trim.py
rename to dart/api/internal/trim.py
index fc2a15a713e5eb8be8973810304341bfe875b9a3..34b3d0b6e6e6c80f9a510d0fd3e55fdb6c4bc482 100644
--- a/dart/scripts/trim.py
+++ b/dart/api/internal/trim.py
@@ -25,7 +25,7 @@ import math
import dart.utils as dU
import tbox.utils as tU
from fwk.coloring import ccolors
-from dart.scripts.config import Config
+from .config import Config
class Trim(Config):
"""Utility to perform the trim analysis of a given lifting configuration
diff --git a/dart/interfaces/mphys.py b/dart/api/mphys.py
similarity index 97%
rename from dart/interfaces/mphys.py
rename to dart/api/mphys.py
index 8de8b09b753ff25fc10473a4553659d70f4eb121..06f0ca06d1b5598cdd47d6b1ec57f7f0197cc40f 100644
--- a/dart/interfaces/mphys.py
+++ b/dart/api/mphys.py
@@ -1,662 +1,662 @@
-# -*- 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.
-
-## MPHYS interface
-# Adrien Crovato
-#
-# Structure:
-# - AeroBuilder
-# - AeroSurfaceMesh
-# - AeroGroup
-# - Morpher
-# - Solver
-# - Loads
-# - AeroCoefficients
-#
-# Remarks:
-# - for each component that uses the openMDAO matrix-free API, the partial
-# gradients are first computed internally in dart::Adjoint, through
-# om.ImplicitComponent.linearize() or om.ExplicitComponent.compute_partials(),
-# and the matrix-vector product is then computed in dart:: Adjoint through
-# om.ImplicitComponent.apply_linear() or om.ExplicitComponent.compute_jacvec_product().
-# This allows faster gradient evaluations, at the cost of a reasonable
-# increase in memory storage, without interfacing any sparse matrices through
-# python.
-
-import numpy as np
-import openmdao.api as om
-from mphys.builder import Builder
-
-# Surface mesh
-class DartMesh(om.IndepVarComp):
- """Initial surface mesh of moving body
- For compatibility, the node coordinates of connected surfaces are always 3D
-
- Attributes
- ----------
- dim : int
- Problem dimension
- bnd : dart.Body object
- Moving body
- """
- def initialize(self):
- self.options.declare('dim', desc='problem dimension')
- self.options.declare('bnd', desc='moving body', recordable=False)
-
- def setup(self):
- self.bnd = self.options['bnd']
- # Store all the node coordinates and add as output
- x_aero0 = np.zeros(3 * len(self.bnd.nodes)) # surface node coordinates are 3D
- for i in range(len(self.bnd.nodes)):
- for j in range(3):
- x_aero0[3 * i + j] = self.bnd.nodes[i].pos[j]
- self.add_output('x_aero0', val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates'])
-
- def mphys_get_triangulated_surface(self):
- """Triangulate the surface
- Only for 2D surfaces in 3D problems
- Return a list containing a point and two vectors for each surface element
- o p0
- / \
- v0 /___\ v1
- """
- if self.options['dim'] == 2:
- raise RuntimeError('DartMesh - cannot triangulate 1D surface!\n')
- # Create list of point coordinates and vector components
- p0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
- v0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
- v1 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
- for i in range(len(self.bnd.groups[0].tag.elems)):
- e = self.bnd.groups[0].tag.elems[i]
- p0i = e.nodes[0].pos
- v0i = e.nodes[1].pos - e.nodes[0].pos
- v1i = e.nodes[2].pos - e.nodes[0].pos
- for j in range(3):
- p0[i][j] = p0i[j]
- v0[i][j] = v0i[j]
- v1[i][j] = v1i[j]
- return [p0, v0, v1]
-
-# Mesh morphers
-class DartMorpher(om.ImplicitComponent):
- """Volume mesh morpher
- For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
-
- Attributes
- ----------
- dim : int
- Problem dimension, used to set the length of the volume node coordinates vector
- bnd : dart.Body object
- Moving body
- mrf : tbox.MshDeform object
- Mesh morpher
- adj : dart.Adjoint object
- Adjoint solver
- """
- def initialize(self):
- self.options.declare('dim', desc='problem dimension')
- self.options.declare('bnd', desc='moving body', recordable=False)
- self.options.declare('mrf', desc='mesh morpher', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
-
- def setup(self):
- self.dim = self.options['dim']
- self.bnd = self.options['bnd']
- self.mrf = self.options['mrf']
- self.adj = self.options['adj']
- # I/O
- self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
- self.add_output('xv', val=np.zeros(self.dim * len(self.mrf.msh.nodes)), desc='aerodynamic volume node coordinates', tags=['mphys_coupling']) # volume node coordinates can be 2D or 3D
- # Partials
- # self.declare_partials(of=['xv'], wrt=['x_aero'])
-
- def solve_nonlinear(self, inputs, outputs):
- """Deform the volume mesh
- """
- # Update inputs
- self.mrf.savePos()
- for i in range(len(self.bnd.nodes)):
- for j in range(3):
- self.bnd.nodes[i].pos[j] = inputs['x_aero'][3 * i + j]
- # Compute
- self.mrf.deform()
- # Update outputs
- for i in range(len(self.mrf.msh.nodes)):
- for j in range(self.dim):
- outputs['xv'][self.dim * i + j] = self.mrf.msh.nodes[i].pos[j]
-
- def apply_nonlinear(self, inputs, outputs, residuals):
- """Compute the residuals of the mesh coordinates
- """
- raise NotImplementedError('DartMorpher - cannot compute the residuals!\n')
-
- def linearize(self, inputs, outputs, partials):
- """Compute the mesh jacobian matrix
- """
- self.adj.linearizeMesh()
-
- def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
- """Perform the matrix-vector product using the mesh Jacobian
- """
- if mode == 'rev':
- if 'xv' in d_residuals:
- if 'xv' in d_outputs:
- d_outputs['xv'] += self.adj.computeJacVecMesh(d_residuals['xv']) # dX = dRx_dX^T * dRx
- if 'x_aero' in d_inputs:
- for i in range(len(self.bnd.nodes)):
- for j in range(self.dim):
- d_inputs['x_aero'][3 * i + j] -= d_residuals['xv'][self.dim * self.bnd.nodes[i].row + j] # dXs = -dX
- elif mode == 'fwd':
- raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
-
- def solve_linear(self, d_outputs, d_residuals, mode):
- """Solve for the mesh residuals using the mesh Jacobian
- """
- if mode == 'rev':
- d_residuals['xv'] = self.adj.solveMesh(d_outputs['xv']) # dRx = dRx_dX^-T * dX
- elif mode == 'fwd':
- raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
-
-class DartDummyMorpher(om.ExplicitComponent):
- """Dummy volume mesh morpher
- Allow to define and link the surface mesh coordinates to the volume mesh coordinates, in case the mesh do not deform.
- For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
- """
- def initialize(self):
- self.options.declare('dim', desc='problem dimension')
- self.options.declare('msh', desc='mesh data structure', recordable=False)
-
- def setup(self):
- dim = self.options['dim']
- msh = self.options['msh']
- # Set initial values
- xv = np.zeros(dim * len(msh.nodes)) # volume node coordinates can be 2D or 3D
- for i in range(len(msh.nodes)):
- for j in range(dim):
- xv[dim * i + j] = msh.nodes[i].pos[j]
- # I/O
- self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
- self.add_output('xv', val=xv, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
-
- def compute(self, inputs, outputs):
- """DO NOT deform the volume mesh
- """
- pass
-
- def compute_partials(self, inputs, partials):
- """DO NOT compute the partials derivatives
- """
- pass
-
-# Solver
-class DartSolver(om.ImplicitComponent):
- """Aerodynamic solver
-
- Attributes
- ----------
- sol : dart.Newton object
- Direct Newton solver
- adj : dart.Adjoint object
- Adjoint solver
- """
- def initialize(self):
- self.options.declare('sol', desc='direct solver', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
-
- def setup(self):
- self.sol = self.options['sol']
- self.adj = self.options['adj']
- # I/O
- self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
- self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
- self.add_output('phi', val=np.zeros(len(self.sol.pbl.msh.nodes)), desc='dart variables (potential)', tags=['mphys_coupling'])
- # Partials
- # self.declare_partials(of=['phi'], wrt=['aoa', 'xv', 'phi'])
-
- def solve_nonlinear(self, inputs, outputs):
- """Compute the flow variables
- """
- # Update inputs
- # NB: the mesh is already up-to-date mesh already up-to-date because DartMorpher has already been run
- self.sol.pbl.update(inputs['aoa'][0])
- # Compute
- self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
- self.sol.run()
- # Update outputs
- outputs['phi'] = self.sol.phi
-
- def apply_nonlinear(self, inputs, outputs, residuals):
- """Compute the residuals of the flow variables
- """
- raise NotImplementedError('DartSolver - cannot compute the residuals!\n')
-
- def linearize(self, inputs, outputs, partials):
- """Compute the flow jacobian matrix (and the other partials of the flow residuals)
- """
- self.adj.linearizeFlow()
-
- def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
- """Perform the matrix-vector product using the partial gradients of the flow residuals
- """
- if mode == 'rev':
- if 'phi' in d_residuals:
- if 'phi' in d_outputs:
- d_outputs['phi'] += self.adj.computeJacVecFlow(d_residuals['phi']) # dU = dRu_dU^T * dRu
- if 'aoa' in d_inputs:
- d_inputs['aoa'] += self.adj.computeJacVecFlowAoa(d_residuals['phi']) # dA = dRu_dA^T * dRu
- if 'xv' in d_inputs:
- d_inputs['xv'] += self.adj.computeJacVecFlowMesh(d_residuals['phi']) # dX = dRu_dX^T * dRu
- elif mode == 'fwd':
- raise NotImplementedError('DartSolver - forward mode not implemented!\n')
-
- def solve_linear(self, d_outputs, d_residuals, mode):
- """Solve for the flow residuals using the flow Jacobian
- """
- if mode == 'rev':
- d_residuals['phi'] = self.adj.solveFlow(d_outputs['phi']) # dRu = dRu_dU^-T * dU
- elif mode == 'fwd':
- raise NotImplementedError('DartSolver - forward mode not implemented!\n')
-
-# Aerodynamic loads
-class DartLoads(om.ExplicitComponent):
- """Aerodynamic loads on moving body
- For compatibility, the forces of connected surfaces are always 3D
-
- Attributes
- ----------
- qinf : dart.Newton object
- Freestream dynamic pressure
- bnd : dart.Body object
- Moving body
- adj : dart.Adjoint object
- Adjoint solver
- """
- def initialize(self):
- self.options.declare('qinf', desc='freestream dynamic pressure')
- self.options.declare('bnd', desc='moving body', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
-
- def setup(self):
- self.qinf = self.options['qinf']
- self.bnd = self.options['bnd']
- self.adj = self.options['adj']
- # I/O
- self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
- self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
- self.add_output('f_aero', val=np.zeros(3 * len(self.bnd.nodes)), desc='aerodynamic loads', tags=['mphys_coupling'])
- # Partials
- # self.declare_partials(of=['f_aero'], wrt=['xv', 'phi'])
-
- def compute(self, inputs, outputs):
- """Get the forces on moving body
- """
- # NB: inputs already up-to-date because DartSolver has already been run
- for i in range(len(self.bnd.nodes)):
- for j in range(3):
- outputs['f_aero'][3 * i + j] = self.qinf * self.bnd.nLoads[i][j]
-
- def compute_partials(self, inputs, partials):
- """Compute the partials gradients of the loads
- """
- self.adj.linearizeLoads(self.bnd)
-
- def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
- """Perform the matrix-vector product using the partial gradients of the loads
- """
- if mode == 'rev':
- if 'f_aero' in d_outputs:
- if 'xv' in d_inputs:
- d_inputs['xv'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsMesh(d_outputs['f_aero'])) # dX = dL_dX^T * dL
- if 'phi' in d_inputs:
- d_inputs['phi'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsFlow(d_outputs['f_aero'])) # dU = dL_dU^T * dL
- elif mode == 'fwd':
- raise NotImplementedError('DartLoads - forward mode not implemented!\n')
-
-# Aerodynamic coefficients
-class DartCoefficients(om.ExplicitComponent):
- """Aerodynamic load coefficients for full aircraft
- Also write solution to disk
-
- Attributes
- ----------
- iscn : int
- ID of scenario (default: 0)
- sol : dart.Newton object
- Direct Newton sovler
- adj : dart.Adjoint object
- Adjoint solver
- """
- def initialize(self):
- self.options.declare('sol', desc='direct solver', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
- self.options.declare('wrtr', desc='data writer', recordable=False)
- self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
-
- def setup(self):
- self.iscn = 0 # id of scenario
- self.sol = self.options['sol']
- self.adj = self.options['adj']
- # I/O
- self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
- self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
- self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
- self.add_output('cl', val=0., desc='lift coefficient', tags=['mphys_result'])
- self.add_output('cd', val=0., desc='drag coefficient', tags=['mphys_result'])
- # Partials
- self.declare_partials(of=['cl', 'cd'], wrt=['aoa', 'xv', 'phi'])
-
- def mphys_set_iscn(self, iscn):
- """Set the id of the scenario
- Allow to save the mesh and results in different files
- """
- self.iscn = iscn
-
- def compute(self, inputs, outputs):
- """Get the coefficients for the full aircraft
- """
- # NB: inputs already up-to-date because DartSolver has already been run
- # Write to disk
- if self.options['morph'] and str(type(self.options['wrtr'])) == '<class \'tboxw.GmshExport\'>':
- if self.iscn == 0:
- self.options['wrtr'].save(self.sol.pbl.msh.name)
- else:
- self.options['wrtr'].save(self.sol.pbl.msh.name + '_{0:d}'.format(self.iscn))
- self.sol.save(self.options['wrtr'], self.iscn)
- # Update outputs
- outputs['cl'] = self.sol.Cl
- outputs['cd'] = self.sol.Cd
-
- def compute_partials(self, inputs, partials):
- """Compute the partials gradients of the load coefficients
- """
- # Gradients wrt to angle of attack
- dCfAoa = self.adj.computeGradientCoefficientsAoa()
- partials['cl', 'aoa'] = dCfAoa[0]
- partials['cd', 'aoa'] = dCfAoa[1]
- # Gradients wrt to mesh
- dCfMsh = self.adj.computeGradientCoefficientsMesh()
- partials['cl', 'xv'] = dCfMsh[0]
- partials['cd', 'xv'] = dCfMsh[1]
- # Gradients wrt to flow variables
- dCfPhi = self.adj.computeGradientCoefficientsFlow()
- partials['cl', 'phi'] = dCfPhi[0]
- partials['cd', 'phi'] = dCfPhi[1]
-
-# Aerodynamic group
-class DartGroup(om.Group):
- """Aerodynamic group
- Integrate the aerodynamic computations in an aero-structural coupling procedure
-
- Components
- ----------
- - Mesh morpher
- - Direct and adjoint solvers
- - Loads computation
- """
- def initialize(self):
- self.options.declare('qinf', desc='freestream dynamic pressure')
- self.options.declare('bnd', desc='moving body', recordable=False)
- self.options.declare('sol', desc='direct solver', recordable=False)
- self.options.declare('mrf', desc='mesh morpher', recordable=False)
- self.options.declare('adj', desc='adjoint solver', recordable=False)
-
- def setup(self):
- # Components
- if self.options['mrf'] is None:
- self.add_subsystem('morpher', DartDummyMorpher(dim = self.options['sol'].pbl.nDim, msh = self.options['sol'].pbl.msh), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
- else:
- self.add_subsystem('morpher', DartMorpher(dim = self.options['sol'].pbl.nDim, bnd = self.options['bnd'], mrf = self.options['mrf'], adj = self.options['adj']), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
- self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
- if self.options['qinf'] != 0:
- self.add_subsystem('loads', DartLoads(qinf = self.options['qinf'], bnd = self.options['bnd'], adj = self.options['adj']), promotes_inputs=['xv', 'phi'], promotes_outputs=['f_aero'])
-
-# Builder
-class DartBuilder(Builder):
- """Dart builder for MPHYS
-
- Attributes
- ----------
- __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
- __sol : dart.Newton object
- Direct Newton solver
- __adj : dart.Adjoint object
- Adjoint solver
- """
- def __init__(self, params, scenario='aerodynamic', task='analysis'):
- """Instantiate and initialize all the solver components.
- Because we do not use MPI, we do not use the initialize method.
-
- Parameters
- ----------
- params : dict
- Dictonnary of parameters for the mesh, solver, ...
- scenario : string, optional
- Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
- task : string, optional
- Type of task (available: analysis, optimization) (default: analysis)
- """
- # Initialize
- print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
-
- # Imports
- from dartflo import tbox
- from dartflo.tbox import gmsh
- from dartflo 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:
- self.__dim = params['Dim']
- # aoa and aos
- if 'AoA' in params:
- alpha = params['AoA'] * np.pi / 180
- else:
- alpha = 0
- if self.__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'] * np.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
- # 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 "' + scenario + '" was given!\n')
- # dynamic pressure
- if scenario == 'aerostructural':
- self.__qinf = params['Qinf']
- else:
- self.__qinf = 0
-
- # Mesh creation
- self.__msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
- # add the wake
- if self.__dim == 2:
- mshCrck = tbox.MshCrack(self.__msh, self.__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(self.__msh, self.__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(self.__msh).save(self.__msh.name)
- del mshCrck
- self.__wrtr = Writer(self.__msh)
- print('\n')
-
- # Mesh morpher creation
- if scenario == 'aerostructural' or task == 'optimization':
- self.__mrf = tbox.MshDeform(self.__msh, self.__dim)
- self.__mrf.nthreads = nthrd
- self.__mrf.setField(params['Fluid'])
- self.__mrf.addFixed(params['Farfield'])
- if self.__dim == 2:
- self.__mrf.addMoving([params['Wing']])
- self.__mrf.addInternal([params['Wake'], params['Wake']+'_'])
- else:
- if 'Fuselage' in params:
- self.__mrf.addFixed(params['Fuselage'])
- self.__mrf.setSymmetry(params['Symmetry'], 1)
- for i in range(len(params['Wings'])):
- if i == 0:
- self.__mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
- else:
- self.__mrf.addFixed([params['Wings'][i]])
- self.__mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
- self.__mrf.initialize()
- else:
- self.__mrf = None
-
- # Problem creation
- self.__pbl = dart.Problem(self.__msh, self.__dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
- # add the field
- self.__pbl.set(dart.Fluid(self.__msh, params['Fluid'], params['M_inf'], self.__dim, alpha, beta))
- # add the initial condition
- self.__pbl.set(dart.Initial(self.__msh, params['Fluid'], self.__dim, alpha, beta))
- # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
- for i in range (len(params['Farfield'])):
- self.__pbl.add(dart.Freestream(self.__msh, params['Farfield'][i], self.__dim, alpha, beta))
- # add solid boundaries
- if self.__dim == 2:
- bnd = dart.Body(self.__msh, [params['Wing'], params['Fluid']])
- self.__bnd = bnd
- self.__pbl.add(bnd)
- else:
- self.__bnd = None
- for bd in params['Wings']:
- bnd = dart.Body(self.__msh, [bd, params['Fluid']])
- if self.__bnd is None:
- self.__bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
- self.__pbl.add(bnd)
- if 'Fuselage' in params:
- bnd = dart.Body(self.__msh, [params['Fuselage'], params['Fluid']])
- self.__pbl.add(bnd)
- # add Wake/Kutta boundary conditions
- if self.__dim == 2:
- self.__pbl.add(dart.Wake(self.__msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
- self.__pbl.add(dart.Kutta(self.__msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
- else:
- for i in range(len(params['Wakes'])):
- self.__pbl.add(dart.Wake(self.__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
- self.__sol = dart.Newton(self.__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)
- self.__adj = dart.Adjoint(self.__sol, self.__mrf, alinsol)
- else:
- self.__adj = None
-
- # Finalize
- print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
-
- def get_mesh_coordinate_subsystem(self):
- """Return openMDAO component to get the initial surface mesh coordinates
- """
- return DartMesh(dim = self.__dim, bnd = self.__bnd)
-
- def get_coupling_group_subsystem(self):
- """Return openMDAO group containing the morpher, solver and loads
- """
- return DartGroup(qinf = self.__qinf, bnd = self.__bnd, sol = self.__sol, mrf = self.__mrf, adj = self.__adj)
-
- def get_post_coupling_subsystem(self):
- """Return openMDAO component that computes the aero coefficients and writes data to disk
- """
- return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, morph = False if self.__mrf is None else True)
-
- def get_number_of_nodes(self):
- """Return the number of surface nodes
- """
- return len(self.__bnd.nodes)
+# -*- 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.
+
+## MPHYS interface
+# Adrien Crovato
+#
+# Structure:
+# - AeroBuilder
+# - AeroSurfaceMesh
+# - AeroGroup
+# - Morpher
+# - Solver
+# - Loads
+# - AeroCoefficients
+#
+# Remarks:
+# - for each component that uses the openMDAO matrix-free API, the partial
+# gradients are first computed internally in dart::Adjoint, through
+# om.ImplicitComponent.linearize() or om.ExplicitComponent.compute_partials(),
+# and the matrix-vector product is then computed in dart:: Adjoint through
+# om.ImplicitComponent.apply_linear() or om.ExplicitComponent.compute_jacvec_product().
+# This allows faster gradient evaluations, at the cost of a reasonable
+# increase in memory storage, without interfacing any sparse matrices through
+# python.
+
+import numpy as np
+import openmdao.api as om
+from mphys.builder import Builder
+
+# Surface mesh
+class DartMesh(om.IndepVarComp):
+ """Initial surface mesh of moving body
+ For compatibility, the node coordinates of connected surfaces are always 3D
+
+ Attributes
+ ----------
+ dim : int
+ Problem dimension
+ bnd : dart.Body object
+ Moving body
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('bnd', desc='moving body', recordable=False)
+
+ def setup(self):
+ self.bnd = self.options['bnd']
+ # Store all the node coordinates and add as output
+ x_aero0 = np.zeros(3 * len(self.bnd.nodes)) # surface node coordinates are 3D
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ x_aero0[3 * i + j] = self.bnd.nodes[i].pos[j]
+ self.add_output('x_aero0', val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates'])
+
+ def mphys_get_triangulated_surface(self):
+ """Triangulate the surface
+ Only for 2D surfaces in 3D problems
+ Return a list containing a point and two vectors for each surface element
+ o p0
+ / \
+ v0 /___\ v1
+ """
+ if self.options['dim'] == 2:
+ raise RuntimeError('DartMesh - cannot triangulate 1D surface!\n')
+ # Create list of point coordinates and vector components
+ p0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ v0 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ v1 = [[0.]*3 for _ in range(len(self.bnd.groups[0].tag.elems))]
+ for i in range(len(self.bnd.groups[0].tag.elems)):
+ e = self.bnd.groups[0].tag.elems[i]
+ p0i = e.nodes[0].pos
+ v0i = e.nodes[1].pos - e.nodes[0].pos
+ v1i = e.nodes[2].pos - e.nodes[0].pos
+ for j in range(3):
+ p0[i][j] = p0i[j]
+ v0[i][j] = v0i[j]
+ v1[i][j] = v1i[j]
+ return [p0, v0, v1]
+
+# Mesh morphers
+class DartMorpher(om.ImplicitComponent):
+ """Volume mesh morpher
+ For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+
+ Attributes
+ ----------
+ dim : int
+ Problem dimension, used to set the length of the volume node coordinates vector
+ bnd : dart.Body object
+ Moving body
+ mrf : tbox.MshDeform object
+ Mesh morpher
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('bnd', desc='moving body', recordable=False)
+ self.options.declare('mrf', desc='mesh morpher', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.dim = self.options['dim']
+ self.bnd = self.options['bnd']
+ self.mrf = self.options['mrf']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+ self.add_output('xv', val=np.zeros(self.dim * len(self.mrf.msh.nodes)), desc='aerodynamic volume node coordinates', tags=['mphys_coupling']) # volume node coordinates can be 2D or 3D
+ # Partials
+ # self.declare_partials(of=['xv'], wrt=['x_aero'])
+
+ def solve_nonlinear(self, inputs, outputs):
+ """Deform the volume mesh
+ """
+ # Update inputs
+ self.mrf.savePos()
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ self.bnd.nodes[i].pos[j] = inputs['x_aero'][3 * i + j]
+ # Compute
+ self.mrf.deform()
+ # Update outputs
+ for i in range(len(self.mrf.msh.nodes)):
+ for j in range(self.dim):
+ outputs['xv'][self.dim * i + j] = self.mrf.msh.nodes[i].pos[j]
+
+ def apply_nonlinear(self, inputs, outputs, residuals):
+ """Compute the residuals of the mesh coordinates
+ """
+ raise NotImplementedError('DartMorpher - cannot compute the residuals!\n')
+
+ def linearize(self, inputs, outputs, partials):
+ """Compute the mesh jacobian matrix
+ """
+ self.adj.linearizeMesh()
+
+ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
+ """Perform the matrix-vector product using the mesh Jacobian
+ """
+ if mode == 'rev':
+ if 'xv' in d_residuals:
+ if 'xv' in d_outputs:
+ d_outputs['xv'] += self.adj.computeJacVecMesh(d_residuals['xv']) # dX = dRx_dX^T * dRx
+ if 'x_aero' in d_inputs:
+ for i in range(len(self.bnd.nodes)):
+ for j in range(self.dim):
+ d_inputs['x_aero'][3 * i + j] -= d_residuals['xv'][self.dim * self.bnd.nodes[i].row + j] # dXs = -dX
+ elif mode == 'fwd':
+ raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
+
+ def solve_linear(self, d_outputs, d_residuals, mode):
+ """Solve for the mesh residuals using the mesh Jacobian
+ """
+ if mode == 'rev':
+ d_residuals['xv'] = self.adj.solveMesh(d_outputs['xv']) # dRx = dRx_dX^-T * dX
+ elif mode == 'fwd':
+ raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
+
+class DartDummyMorpher(om.ExplicitComponent):
+ """Dummy volume mesh morpher
+ Allow to define and link the surface mesh coordinates to the volume mesh coordinates, in case the mesh do not deform.
+ For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+ """
+ def initialize(self):
+ self.options.declare('dim', desc='problem dimension')
+ self.options.declare('msh', desc='mesh data structure', recordable=False)
+
+ def setup(self):
+ dim = self.options['dim']
+ msh = self.options['msh']
+ # Set initial values
+ xv = np.zeros(dim * len(msh.nodes)) # volume node coordinates can be 2D or 3D
+ for i in range(len(msh.nodes)):
+ for j in range(dim):
+ xv[dim * i + j] = msh.nodes[i].pos[j]
+ # I/O
+ self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+ self.add_output('xv', val=xv, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+
+ def compute(self, inputs, outputs):
+ """DO NOT deform the volume mesh
+ """
+ pass
+
+ def compute_partials(self, inputs, partials):
+ """DO NOT compute the partials derivatives
+ """
+ pass
+
+# Solver
+class DartSolver(om.ImplicitComponent):
+ """Aerodynamic solver
+
+ Attributes
+ ----------
+ sol : dart.Newton object
+ Direct Newton solver
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.sol = self.options['sol']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_output('phi', val=np.zeros(len(self.sol.pbl.msh.nodes)), desc='dart variables (potential)', tags=['mphys_coupling'])
+ # Partials
+ # self.declare_partials(of=['phi'], wrt=['aoa', 'xv', 'phi'])
+
+ def solve_nonlinear(self, inputs, outputs):
+ """Compute the flow variables
+ """
+ # Update inputs
+ # NB: the mesh is already up-to-date mesh already up-to-date because DartMorpher has already been run
+ self.sol.pbl.update(inputs['aoa'][0])
+ # Compute
+ self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
+ self.sol.run()
+ # Update outputs
+ outputs['phi'] = self.sol.phi
+
+ def apply_nonlinear(self, inputs, outputs, residuals):
+ """Compute the residuals of the flow variables
+ """
+ raise NotImplementedError('DartSolver - cannot compute the residuals!\n')
+
+ def linearize(self, inputs, outputs, partials):
+ """Compute the flow jacobian matrix (and the other partials of the flow residuals)
+ """
+ self.adj.linearizeFlow()
+
+ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
+ """Perform the matrix-vector product using the partial gradients of the flow residuals
+ """
+ if mode == 'rev':
+ if 'phi' in d_residuals:
+ if 'phi' in d_outputs:
+ d_outputs['phi'] += self.adj.computeJacVecFlow(d_residuals['phi']) # dU = dRu_dU^T * dRu
+ if 'aoa' in d_inputs:
+ d_inputs['aoa'] += self.adj.computeJacVecFlowAoa(d_residuals['phi']) # dA = dRu_dA^T * dRu
+ if 'xv' in d_inputs:
+ d_inputs['xv'] += self.adj.computeJacVecFlowMesh(d_residuals['phi']) # dX = dRu_dX^T * dRu
+ elif mode == 'fwd':
+ raise NotImplementedError('DartSolver - forward mode not implemented!\n')
+
+ def solve_linear(self, d_outputs, d_residuals, mode):
+ """Solve for the flow residuals using the flow Jacobian
+ """
+ if mode == 'rev':
+ d_residuals['phi'] = self.adj.solveFlow(d_outputs['phi']) # dRu = dRu_dU^-T * dU
+ elif mode == 'fwd':
+ raise NotImplementedError('DartSolver - forward mode not implemented!\n')
+
+# Aerodynamic loads
+class DartLoads(om.ExplicitComponent):
+ """Aerodynamic loads on moving body
+ For compatibility, the forces of connected surfaces are always 3D
+
+ Attributes
+ ----------
+ qinf : dart.Newton object
+ Freestream dynamic pressure
+ bnd : dart.Body object
+ Moving body
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('qinf', desc='freestream dynamic pressure')
+ self.options.declare('bnd', desc='moving body', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ self.qinf = self.options['qinf']
+ self.bnd = self.options['bnd']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
+ self.add_output('f_aero', val=np.zeros(3 * len(self.bnd.nodes)), desc='aerodynamic loads', tags=['mphys_coupling'])
+ # Partials
+ # self.declare_partials(of=['f_aero'], wrt=['xv', 'phi'])
+
+ def compute(self, inputs, outputs):
+ """Get the forces on moving body
+ """
+ # NB: inputs already up-to-date because DartSolver has already been run
+ for i in range(len(self.bnd.nodes)):
+ for j in range(3):
+ outputs['f_aero'][3 * i + j] = self.qinf * self.bnd.nLoads[i][j]
+
+ def compute_partials(self, inputs, partials):
+ """Compute the partials gradients of the loads
+ """
+ self.adj.linearizeLoads(self.bnd)
+
+ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
+ """Perform the matrix-vector product using the partial gradients of the loads
+ """
+ if mode == 'rev':
+ if 'f_aero' in d_outputs:
+ if 'xv' in d_inputs:
+ d_inputs['xv'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsMesh(d_outputs['f_aero'])) # dX = dL_dX^T * dL
+ if 'phi' in d_inputs:
+ d_inputs['phi'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsFlow(d_outputs['f_aero'])) # dU = dL_dU^T * dL
+ elif mode == 'fwd':
+ raise NotImplementedError('DartLoads - forward mode not implemented!\n')
+
+# Aerodynamic coefficients
+class DartCoefficients(om.ExplicitComponent):
+ """Aerodynamic load coefficients for full aircraft
+ Also write solution to disk
+
+ Attributes
+ ----------
+ iscn : int
+ ID of scenario (default: 0)
+ sol : dart.Newton object
+ Direct Newton sovler
+ adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('wrtr', desc='data writer', recordable=False)
+ self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
+
+ def setup(self):
+ self.iscn = 0 # id of scenario
+ self.sol = self.options['sol']
+ self.adj = self.options['adj']
+ # I/O
+ self.add_input('aoa', val=0., units='rad', desc='angle of attack', tags=['mphys_input'])
+ self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
+ self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
+ self.add_output('cl', val=0., desc='lift coefficient', tags=['mphys_result'])
+ self.add_output('cd', val=0., desc='drag coefficient', tags=['mphys_result'])
+ # Partials
+ self.declare_partials(of=['cl', 'cd'], wrt=['aoa', 'xv', 'phi'])
+
+ def mphys_set_iscn(self, iscn):
+ """Set the id of the scenario
+ Allow to save the mesh and results in different files
+ """
+ self.iscn = iscn
+
+ def compute(self, inputs, outputs):
+ """Get the coefficients for the full aircraft
+ """
+ # NB: inputs already up-to-date because DartSolver has already been run
+ # Write to disk
+ if self.options['morph'] and str(type(self.options['wrtr'])) == '<class \'tboxw.GmshExport\'>':
+ if self.iscn == 0:
+ self.options['wrtr'].save(self.sol.pbl.msh.name)
+ else:
+ self.options['wrtr'].save(self.sol.pbl.msh.name + '_{0:d}'.format(self.iscn))
+ self.sol.save(self.options['wrtr'], self.iscn)
+ # Update outputs
+ outputs['cl'] = self.sol.Cl
+ outputs['cd'] = self.sol.Cd
+
+ def compute_partials(self, inputs, partials):
+ """Compute the partials gradients of the load coefficients
+ """
+ # Gradients wrt to angle of attack
+ dCfAoa = self.adj.computeGradientCoefficientsAoa()
+ partials['cl', 'aoa'] = dCfAoa[0]
+ partials['cd', 'aoa'] = dCfAoa[1]
+ # Gradients wrt to mesh
+ dCfMsh = self.adj.computeGradientCoefficientsMesh()
+ partials['cl', 'xv'] = dCfMsh[0]
+ partials['cd', 'xv'] = dCfMsh[1]
+ # Gradients wrt to flow variables
+ dCfPhi = self.adj.computeGradientCoefficientsFlow()
+ partials['cl', 'phi'] = dCfPhi[0]
+ partials['cd', 'phi'] = dCfPhi[1]
+
+# Aerodynamic group
+class DartGroup(om.Group):
+ """Aerodynamic group
+ Integrate the aerodynamic computations in an aero-structural coupling procedure
+
+ Components
+ ----------
+ - Mesh morpher
+ - Direct and adjoint solvers
+ - Loads computation
+ """
+ def initialize(self):
+ self.options.declare('qinf', desc='freestream dynamic pressure')
+ self.options.declare('bnd', desc='moving body', recordable=False)
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('mrf', desc='mesh morpher', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+
+ def setup(self):
+ # Components
+ if self.options['mrf'] is None:
+ self.add_subsystem('morpher', DartDummyMorpher(dim = self.options['sol'].pbl.nDim, msh = self.options['sol'].pbl.msh), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
+ else:
+ self.add_subsystem('morpher', DartMorpher(dim = self.options['sol'].pbl.nDim, bnd = self.options['bnd'], mrf = self.options['mrf'], adj = self.options['adj']), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
+ self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
+ if self.options['qinf'] != 0:
+ self.add_subsystem('loads', DartLoads(qinf = self.options['qinf'], bnd = self.options['bnd'], adj = self.options['adj']), promotes_inputs=['xv', 'phi'], promotes_outputs=['f_aero'])
+
+# Builder
+class DartBuilder(Builder):
+ """Dart builder for MPHYS
+
+ Attributes
+ ----------
+ __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
+ __sol : dart.Newton object
+ Direct Newton solver
+ __adj : dart.Adjoint object
+ Adjoint solver
+ """
+ def __init__(self, params, scenario='aerodynamic', task='analysis'):
+ """Instantiate and initialize all the solver components.
+ Because we do not use MPI, we do not use the initialize method.
+
+ Parameters
+ ----------
+ params : dict
+ Dictonnary of parameters for the mesh, solver, ...
+ scenario : string, optional
+ Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
+ task : string, optional
+ Type of task (available: analysis, optimization) (default: analysis)
+ """
+ # Initialize
+ print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
+
+ # Imports
+ from dartflo import tbox
+ from dartflo.tbox import gmsh
+ from dartflo 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:
+ self.__dim = params['Dim']
+ # aoa and aos
+ if 'AoA' in params:
+ alpha = params['AoA'] * np.pi / 180
+ else:
+ alpha = 0
+ if self.__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'] * np.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
+ # 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 "' + scenario + '" was given!\n')
+ # dynamic pressure
+ if scenario == 'aerostructural':
+ self.__qinf = params['Q_inf']
+ else:
+ self.__qinf = 0
+
+ # Mesh creation
+ self.__msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
+ # add the wake
+ if self.__dim == 2:
+ mshCrck = tbox.MshCrack(self.__msh, self.__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(self.__msh, self.__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(self.__msh).save(self.__msh.name)
+ del mshCrck
+ self.__wrtr = Writer(self.__msh)
+ print('\n')
+
+ # Mesh morpher creation
+ if scenario == 'aerostructural' or task == 'optimization':
+ self.__mrf = tbox.MshDeform(self.__msh, self.__dim)
+ self.__mrf.nthreads = nthrd
+ self.__mrf.setField(params['Fluid'])
+ self.__mrf.addFixed(params['Farfield'])
+ if self.__dim == 2:
+ self.__mrf.addMoving([params['Wing']])
+ self.__mrf.addInternal([params['Wake'], params['Wake']+'_'])
+ else:
+ if 'Fuselage' in params:
+ self.__mrf.addFixed(params['Fuselage'])
+ self.__mrf.setSymmetry(params['Symmetry'], 1)
+ for i in range(len(params['Wings'])):
+ if i == 0:
+ self.__mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
+ else:
+ self.__mrf.addFixed([params['Wings'][i]])
+ self.__mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
+ self.__mrf.initialize()
+ else:
+ self.__mrf = None
+
+ # Problem creation
+ self.__pbl = dart.Problem(self.__msh, self.__dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
+ # add the field
+ self.__pbl.set(dart.Fluid(self.__msh, params['Fluid'], params['M_inf'], self.__dim, alpha, beta))
+ # add the initial condition
+ self.__pbl.set(dart.Initial(self.__msh, params['Fluid'], self.__dim, alpha, beta))
+ # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
+ for i in range (len(params['Farfield'])):
+ self.__pbl.add(dart.Freestream(self.__msh, params['Farfield'][i], self.__dim, alpha, beta))
+ # add solid boundaries
+ if self.__dim == 2:
+ bnd = dart.Body(self.__msh, [params['Wing'], params['Fluid']])
+ self.__bnd = bnd
+ self.__pbl.add(bnd)
+ else:
+ self.__bnd = None
+ for bd in params['Wings']:
+ bnd = dart.Body(self.__msh, [bd, params['Fluid']])
+ if self.__bnd is None:
+ self.__bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
+ self.__pbl.add(bnd)
+ if 'Fuselage' in params:
+ bnd = dart.Body(self.__msh, [params['Fuselage'], params['Fluid']])
+ self.__pbl.add(bnd)
+ # add Wake/Kutta boundary conditions
+ if self.__dim == 2:
+ self.__pbl.add(dart.Wake(self.__msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
+ self.__pbl.add(dart.Kutta(self.__msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
+ else:
+ for i in range(len(params['Wakes'])):
+ self.__pbl.add(dart.Wake(self.__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
+ self.__sol = dart.Newton(self.__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)
+ self.__adj = dart.Adjoint(self.__sol, self.__mrf, alinsol)
+ else:
+ self.__adj = None
+
+ # Finalize
+ print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
+
+ def get_mesh_coordinate_subsystem(self):
+ """Return openMDAO component to get the initial surface mesh coordinates
+ """
+ return DartMesh(dim = self.__dim, bnd = self.__bnd)
+
+ def get_coupling_group_subsystem(self):
+ """Return openMDAO group containing the morpher, solver and loads
+ """
+ return DartGroup(qinf = self.__qinf, bnd = self.__bnd, sol = self.__sol, mrf = self.__mrf, adj = self.__adj)
+
+ def get_post_coupling_subsystem(self):
+ """Return openMDAO component that computes the aero coefficients and writes data to disk
+ """
+ return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, morph = False if self.__mrf is None else True)
+
+ def get_number_of_nodes(self):
+ """Return the number of surface nodes
+ """
+ return len(self.__bnd.nodes)
diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py
index 2d5d214bb5d2ab830ec32cd76e6260e138b7c236..e23edac7696331c98e4ed432539fdc088367e4ca 100644
--- a/dart/cases/coyote.py
+++ b/dart/cases/coyote.py
@@ -24,6 +24,7 @@
# Valid mesh obtained with gmsh 4.4.1
def getParam():
+ import os.path
p = {}
# Specific
scl = 2 # scaling factor for lifting surface mesh size
@@ -32,7 +33,7 @@ def getParam():
fums = 0.1 # fuselage mesh size
fams = 5. # farfield mesh size
# Input/Output
- p['File'] = '../models/coyote.geo' # Input file containing the model
+ p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/coyote.geo') # Input file containing the model
p['Pars'] = {'msWing' : wrtems, 'msTail' : trtems,
'msFus' : fums, 'msFar' : fams} # Parameters for input file model
p['Dim'] = 3 # Problem dimension
@@ -75,7 +76,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.polar as p
+ import dart.api.internal.polar as p
polar = p.Polar(getParam())
polar.run()
polar.disp()
diff --git a/dart/cases/n0012.py b/dart/cases/n0012.py
index 5427c47a9ac1581187ef7da55dd1e59afafdc7fd..f984557d0a4ac95309b850fb30f77014a371cf3f 100644
--- a/dart/cases/n0012.py
+++ b/dart/cases/n0012.py
@@ -20,9 +20,10 @@
# Adrien Crovato
def getParam():
+ import os.path
p = {}
# Input/Output
- p['File'] = '../models/n0012.geo' # Input file containing the model
+ 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
p['Dim'] = 2 # Problem dimension
p['Format'] = 'vtk' # Save format (vtk or gmsh)
@@ -54,7 +55,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.polar as p
+ import dart.api.internal.polar as p
polar = p.Polar(getParam())
polar.run()
polar.disp()
diff --git a/dart/cases/n0012_3.py b/dart/cases/n0012_3.py
index 6039734e9e1d761152ab73e1ed8012600412f6f4..11338607c8935df66e8a9c486a850780820e2b6f 100644
--- a/dart/cases/n0012_3.py
+++ b/dart/cases/n0012_3.py
@@ -20,11 +20,12 @@
# Adrien Crovato
def getParam():
+ import os.path
p = {}
# Specific
span = 1. # span length
# Input/Output
- p['File'] = '../models/n0012_3.geo' # Input file containing the model
+ 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)
@@ -62,7 +63,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.polar as p
+ import dart.api.internal.polar as p
polar = p.Polar(getParam())
polar.run()
polar.disp()
diff --git a/dart/cases/n64A410.py b/dart/cases/n64A410.py
index ba98f07bb25fbc11b2c182c82b29610f727a039f..0ef97fcdde17c9cb70329faa74ece232e62552ac 100644
--- a/dart/cases/n64A410.py
+++ b/dart/cases/n64A410.py
@@ -20,9 +20,10 @@
# Adrien Crovato
def getParam():
+ import os.path
p = {}
# Input/Output
- p['File'] = '../models/n64A410.geo' # Input file containing the model
+ 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
p['Dim'] = 2 # Problem dimension
p['Format'] = 'vtk' # Save format (vtk or gmsh)
@@ -56,7 +57,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.trim as t
+ import dart.api.internal.trim as t
trim = t.Trim(getParam())
trim.run()
trim.disp()
diff --git a/dart/cases/rae2822.py b/dart/cases/rae2822.py
index 67bb16d39064575a98c5570568371bfbc65e15df..cf515702b802afe84263d3dbe38362153024b9a6 100644
--- a/dart/cases/rae2822.py
+++ b/dart/cases/rae2822.py
@@ -20,9 +20,10 @@
# Adrien Crovato
def getParam():
+ import os.path
p = {}
# Input/Output
- p['File'] = '../models/rae2822.geo' # Input file containing the model
+ p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/rae2822.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
p['Dim'] = 2 # Problem dimension
p['Format'] = 'vtk' # Save format (vtk or gmsh)
@@ -59,7 +60,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.polar as p
+ import dart.api.internal.polar as p
polar = p.Polar(getParam())
polar.run()
polar.disp()
diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py
index 8744e6a92cbca86c05b98369ee26c40df2a65130..4255fffd60dff15f1a0e5c7b763343f48a1a821e 100644
--- a/dart/cases/wbht.py
+++ b/dart/cases/wbht.py
@@ -23,6 +23,7 @@
# Gmsh might not generate the wing/tail surfaces properly. Results might not be accurate
def getParam():
+ import os.path
p = {}
# Specific
scl = 2 # scaling factor for lifting surface mesh size
@@ -40,7 +41,7 @@ def getParam():
fucms = 0.01 # fuselage caps mesh size
fams = 5. # farfield mesh size
# Input/Output
- p['File'] = '../models/wbht.geo' # Input file containing the model
+ p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wbht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
'msLeW1' : wtlems, 'msTeW1' : wttems,
'msLeW2' : wwlems, 'msTeW2' : wwtems,
@@ -87,7 +88,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.trim as t
+ import dart.api.internal.trim as t
trim = t.Trim(getParam())
trim.run()
trim.disp()
diff --git a/dart/cases/wht.py b/dart/cases/wht.py
index 3c68e0990d4236a1490929f7fc6e922d1defe3ca..c5f90f507ab3b8dd2615c4a156c5598ef7fcd127 100644
--- a/dart/cases/wht.py
+++ b/dart/cases/wht.py
@@ -20,6 +20,7 @@
# Adrien Crovato
def getParam():
+ import os.path
p = {}
# Specific
scl = 2 # scaling factor for lifting surface mesh size
@@ -35,7 +36,7 @@ def getParam():
tttems = scl*0.0058 # tail tip trailing mesh size
fams = 5. # farfield mesh size
# Input/Output
- p['File'] = '../models/wht.geo' # Input file containing the model
+ p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/wht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
'msLeW1' : wtlems, 'msTeW1' : wttems,
'msLeW2' : wwlems, 'msTeW2' : wwtems,
@@ -78,7 +79,7 @@ def getParam():
def main():
# Script call
- import dart.scripts.trim as t
+ import dart.api.internal.trim as t
trim = t.Trim(getParam())
trim.run()
trim.disp()
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index af6135b7881fcc6b928787fbe1b38703100ab5d8..1d801953e86634863fba0c643f8c2c3c633c975f 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -89,8 +89,8 @@ bool Newton::run()
// Display current freestream conditions
std::cout << std::setprecision(2)
<< "- Mach " << pbl->M_inf << ", "
- << pbl->alpha * 180 / 3.14159 << "° AoA, "
- << pbl->beta * 180 / 3.14159 << "° AoS"
+ << pbl->alpha * 180 / 3.14159 << " deg AoA, "
+ << pbl->beta * 180 / 3.14159 << " deg AoS"
<< std::endl;
// Initialize solver loop
diff --git a/dart/src/wPicard.cpp b/dart/src/wPicard.cpp
index 03981f908dfeb94ba60f267941feb9a249717050..ea218f7d4c02d3dd4e14c3816776fef51755279b 100644
--- a/dart/src/wPicard.cpp
+++ b/dart/src/wPicard.cpp
@@ -83,8 +83,8 @@ bool Picard::run()
// Display current freestream conditions
std::cout << std::setprecision(2)
<< "- Mach " << pbl->M_inf << ", "
- << pbl->alpha * 180 / 3.14159 << "° AoA, "
- << pbl->beta * 180 / 3.14159 << "° AoS"
+ << pbl->alpha * 180 / 3.14159 << " deg AoA, "
+ << pbl->beta * 180 / 3.14159 << " deg AoS"
<< std::endl;
// Initialize solver loop