diff --git a/CMakeLists.txt b/CMakeLists.txt
index c342e270cdbf091f67ca403443ad4792a7009dc5..ec130f661a7d716a7ad8c474bbb159802d10ceb7 100644
--- a/CMakeLists.txt
+++ b/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,17 +56,12 @@ 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)
diff --git a/dart/_src/CMakeLists.txt b/dart/_src/CMakeLists.txt
index c78c46464f4be1f8f7788ca1f565e3e6e7a5e4ba..9cbe7fffcd239bac03021c362fee47f3b2c41b59 100644
--- a/dart/_src/CMakeLists.txt
+++ b/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(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)
-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}
-)
+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})
+INSTALL(TARGETS dartw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/dart/api/core.py b/dart/api/core.py
index 248a6a10a88877f11746837a59529bbc70d447f5..8359e35b75e2f6d53ff24ec03bb9c6cd26a4cee2 100644
--- a/dart/api/core.py
+++ b/dart/api/core.py
@@ -53,11 +53,10 @@ def initDart(cfg, scenario='aerodynamic', task='analysis'):
Adjoint solver
"""
# Imports
- import math
- import fwk.wutils as wu
+ import dart
import tbox
import tbox.gmsh as gmsh
- import dart
+ import math
# Basic checks and config
# dimension
@@ -101,12 +100,11 @@ def initDart(cfg, scenario='aerodynamic', task='analysis'):
nthrd = cfg['Threads']
else:
nthrd = 1
- wu.initMKL(nthrd) # initialize threading layer and number of threads
# verbosity
if 'Verb' in cfg:
verb = cfg['Verb']
else:
- verb = 0
+ 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')
@@ -124,17 +122,21 @@ def initDart(cfg, scenario='aerodynamic', task='analysis'):
if _dim == 2:
mshCrck = tbox.MshCrack(_msh, _dim)
mshCrck.setCrack(cfg['Wake'])
- mshCrck.addBoundaries([cfg['Fluid'], cfg['Farfield'][-1], cfg['Wing']])
+ 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.addBoundaries([cfg['Fluid'], cfg['Farfield'][-1], cfg['Wings'][i]])
+ mshCrck.addBoundary(cfg['Fluid'])
+ mshCrck.addBoundary(cfg['Farfield'][-1])
+ mshCrck.addBoundary(cfg['Wings'][i])
if 'Fuselage' in cfg:
- mshCrck.addBoundaries([cfg['Fuselage']])
+ mshCrck.addBoundary(cfg['Fuselage'])
if 'Symmetry' in cfg:
- mshCrck.addBoundaries([cfg['Symmetry']])
+ mshCrck.addBoundary(cfg['Symmetry'])
if 'WakeTips' in cfg:
mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
mshCrck.run()
@@ -146,22 +148,24 @@ def initDart(cfg, scenario='aerodynamic', task='analysis'):
# Mesh morpher creation
if scenario == 'aerostructural' or task == 'optimization':
- _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd)
+ _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd, vrb=verb)
_mrf.setField(cfg['Fluid'])
- _mrf.addFixed(cfg['Farfield'])
+ for tag in cfg['Farfield']:
+ _mrf.addFixed(tag)
if _dim == 2:
- _mrf.addMoving([cfg['Wing']])
- _mrf.addInternal([cfg['Wake'], cfg['Wake']+'_'])
+ _mrf.addMoving(cfg['Wing'])
+ _mrf.addInternal(cfg['Wake'], cfg['Wake']+'_')
else:
if 'Fuselage' in cfg:
_mrf.addFixed(cfg['Fuselage'])
- _mrf.setSymmetry(cfg['Symmetry'], 1)
+ 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
+ _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.addFixed(cfg['Wings'][i])
+ _mrf.addInternal(cfg['Wakes'][i], cfg['Wakes'][i]+'_')
_mrf.initialize()
else:
_mrf = None
diff --git a/dart/api/internal/polar.py b/dart/api/internal/polar.py
index 69f5968bda90344dcc97352f842cf217c33ceb1b..f4355f3e0f4e35f1e2832cd6263982971a96a6e5 100644
--- a/dart/api/internal/polar.py
+++ b/dart/api/internal/polar.py
@@ -81,34 +81,33 @@ class Polar:
"""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:
+ for i in range(len(self.alphas)):
# define current angle of attack
- alpha = alpha*math.pi/180
- # update problem
+ 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', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
+ 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.wrtr, ac)
+ self.sol.save(self.wrtr, acs)
# 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', '')
+ tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.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)
+ dU.writeSlices(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)
- # end loop
- ac+=1
def disp(self):
"""Display the results and draw the polar
@@ -127,6 +126,6 @@ class Polar:
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", "")
+ tU.plot(self.alphas, self.Cls, 'alpha', 'Cl', '')
+ tU.plot(self.alphas, self.Cds, 'alpha', 'Cd', '')
+ tU.plot(self.alphas, self.Cms, 'alpha', 'Cm', '')
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index 74c0025882f5de8f6be7732433c8500b8a5b109d..b277d9f8d85c2d371f6efeda3187d3017832c70f 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -86,8 +86,9 @@ class Trim:
while True:
# define current angle of attack
print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
- # update problem
+ # 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()
@@ -112,7 +113,7 @@ class Trim:
# extract Cp
if self.dim == 2:
Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
- tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "")
+ tU.write(Cp, f'Cp_{self.msh.name}_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)
diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index e59d6130a610c28eca75414bd025c3edf34aaa54..d23c81e436631b6d30747762b29c13479a8742a6 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -218,7 +218,7 @@ class DartSolver(om.ImplicitComponent):
def initialize(self):
self.options.declare('sol', desc='direct solver', recordable=False)
self.options.declare('adj', desc='adjoint solver', recordable=False)
- self.options.declare('raiseError', default=False, desc='raise AnalysisError when solver not fully converged')
+ self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged')
def setup(self):
self.sol = self.options['sol']
@@ -239,7 +239,7 @@ class DartSolver(om.ImplicitComponent):
# Compute
self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
status = self.sol.run()
- if (status == 1 and self.options['raiseError']) or status == 2: # max number of iterations exceeded or NaN in solution
+ if (status == 1 and self.options['rerr']) or status == 2: # max number of iterations exceeded or NaN in solution
raise om.AnalysisError('DART solver failed to fully converge!\n')
# Update outputs
outputs['phi'] = self.sol.phi
@@ -335,12 +335,9 @@ class DartLoads(om.ExplicitComponent):
# 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
@@ -349,12 +346,8 @@ class DartCoefficients(om.ExplicitComponent):
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('iscn', default=0, desc='ID of the scenario')
- self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
def setup(self):
- self.iscn = self.options['iscn']
self.sol = self.options['sol']
self.adj = self.options['adj']
# I/O
@@ -370,13 +363,6 @@ class DartCoefficients(om.ExplicitComponent):
"""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 self.options['wrtr'].type() == 1:
- 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
@@ -397,9 +383,53 @@ class DartCoefficients(om.ExplicitComponent):
partials['cl', 'phi'] = dCfPhi[0]
partials['cd', 'phi'] = dCfPhi[1]
-# Aerodynamic group
+# Aerodynamic writer
+class DartWriter(om.ExplicitComponent):
+ """Write solution to disk
+
+ Attributes
+ ----------
+ sol : dart.Newton object
+ Direct Newton sovler
+ wrt : dart.MshExport object
+ Mesh and solution data writer
+ sname : string
+ Name of scenario
+ sall : bool
+ Whether to save results at each iteration to different files
+ it : int
+ Optimization iteration count
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('wrt', desc='data writer', recordable=False)
+ self.options.declare('snam', desc='name of the scenario')
+ self.options.declare('sall', desc='whether to save results at each iteration to different files')
+
+ def setup(self):
+ self.sol = self.options['sol']
+ self.wrt = self.options['wrt']
+ self.snam = self.options['snam'] if self.options['snam'] is not None else 'default'
+ self.sall = self.options['sall']
+ self.it = 0 # optimization iteration count
+
+ def compute(self, inputs, outputs):
+ """Write to disk
+ """
+ # Create suffix
+ suffix = '_{:s}'.format(self.snam)
+ if self.sall:
+ suffix += '_{:04d}'.format(self.it)
+ self.it += 1
+ # Write mesh to disk (.msh only)
+ if self.wrt.type() == 1:
+ self.wrt.save(self.sol.pbl.msh.name + suffix)
+ # Write solution to disk
+ self.sol.save(self.wrt, suffix)
+
+# Aerodynamic coupling group
class DartGroup(om.Group):
- """Aerodynamic group
+ """Aerodynamic coupling group
Integrate the aerodynamic computations in an aero-structural coupling procedure
Components
@@ -414,17 +444,39 @@ class DartGroup(om.Group):
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)
- self.options.declare('raiseError', default=False, desc='raise AnalysisError when solver not fully converged')
+ self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged')
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'])
+ 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'], raiseError = self.options['raiseError']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
+ 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'], rerr=self.options['rerr']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
if self.options['qinf'] is not None:
- 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'])
+ 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'])
+
+# Aerodynamic post-coupling group
+class DartPostGroup(om.Group):
+ """Aerodynamic post-coupling group
+ Update the aerodynamic load coefficients and write solution to disk
+
+ Components
+ ----------
+ - Coefficients
+ - Writer
+ """
+ def initialize(self):
+ self.options.declare('sol', desc='direct solver', recordable=False)
+ self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('wrt', desc='data writer', recordable=False)
+ self.options.declare('snam', desc='name of the scenario')
+ self.options.declare('sall', desc='whether to save results at each iteration to different files')
+
+ def setup(self):
+ # Components
+ self.add_subsystem('coeff', DartCoefficients(sol=self.options['sol'], adj=self.options['adj']), promotes_inputs=['aoa', 'xv', 'phi'], promotes_outputs=['cl', 'cd'])
+ self.add_subsystem('writer', DartWriter(sol=self.options['sol'], wrt=self.options['wrt'], snam=self.options['snam'], sall=self.options['sall']))
# Builder
class DartBuilder(Builder):
@@ -436,8 +488,8 @@ class DartBuilder(Builder):
Dimension of the problem
__qinf : float
Freestream dynamic pressure
- __wrtr : tbox.MshExport object
- Mesh writer
+ __wrt : tbox.MshExport object
+ Mesh and solution data writer
__mrf : tbox.MshDeform object
Mesh morpher
__sol : dart.Newton object
@@ -445,7 +497,7 @@ class DartBuilder(Builder):
__adj : dart.Adjoint object
Adjoint solver
"""
- def __init__(self, cfg, scenario='aerodynamic', task='analysis', scenarioId=0, raiseError=False):
+ def __init__(self, cfg, scenario='aerodynamic', task='analysis', saveAll=False, raiseError=False):
"""Instantiate the builder and save the parameters and options
Parameters
@@ -456,8 +508,8 @@ class DartBuilder(Builder):
Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
task : string, optional
Type of task (available: analysis, optimization) (default: analysis)
- scenarioId : int, optional
- ID of the scenario, will be used in the output filename
+ saveAll : bool, optional
+ Whether to write results to different files at each iteration (default: False, results file will be overwritten)
raiseError : bool, optional
Whether to raise an openMDAO.AnalysisError when the solver does not fully converge (default: False)
As of openMDAO V3.9.2, AnalysisError are only handled by pyOptSparse
@@ -465,8 +517,8 @@ class DartBuilder(Builder):
self.__cfg = cfg
self.__scn = scenario
self.__tsk = task
- self.__sid = scenarioId
- self.__raiseError = raiseError
+ self.__sall = saveAll
+ self.__rerr = raiseError
def initialize(self, comm):
"""Instantiate and initialize all the solver components
@@ -478,7 +530,7 @@ class DartBuilder(Builder):
"""
from .core import initDart
def _init():
- self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
+ self.__dim, self.__qinf, _, self.__wrt, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
# If n MPI processes, init DART on rank=0,...,n-1 sequentially to avoid issues with mesh IO
# If mpi4py is not available or only 1 MPI process, init DART as usual
try:
@@ -500,20 +552,20 @@ class DartBuilder(Builder):
except:
_init()
- def get_mesh_coordinate_subsystem(self):
+ def get_mesh_coordinate_subsystem(self, scenario_name=None):
"""Return openMDAO component to get the initial surface mesh coordinates
"""
- return DartMesh(dim = self.__dim, bnd = self.__bnd)
+ return DartMesh(dim=self.__dim, bnd=self.__bnd)
- def get_coupling_group_subsystem(self):
+ def get_coupling_group_subsystem(self, scenario_name=None):
"""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, raiseError = self.__raiseError)
+ return DartGroup(qinf=self.__qinf, bnd=self.__bnd, sol=self.__sol, mrf=self.__mrf, adj=self.__adj, rerr=self.__rerr)
- def get_post_coupling_subsystem(self):
- """Return openMDAO component that computes the aero coefficients and writes data to disk
+ def get_post_coupling_subsystem(self, scenario_name=None):
+ """Return openMDAO group that computes the aero coefficients and writes data to disk
"""
- return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, iscn = self.__sid, morph = False if self.__mrf is None else True)
+ return DartPostGroup(sol=self.__sol, adj=self.__adj, wrt=self.__wrt, snam=scenario_name, sall=self.__sall)
def get_number_of_nodes(self):
"""Return the number of surface nodes
diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py
index b0f9c2157829c582e57c204ac18d53e18246d80d..36f048aa3666e342cda7650de1eddc6e49a00024 100644
--- a/dart/benchmark/onera.py
+++ b/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,10 +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(2, 1e-6, 50, 1e-5)
- newton = dart.Newton(pbl, gmres, nthrds=k, vrb=2)
+ newton = dart.Newton(pbl, gmres, nthrds=k, vrb=3)
return newton
def main():
diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py
index fc493e30324286c0bd5cf268fcb1b797bd2bda72..aa91bf51130009077de14cee8ecdb12e710b0328 100644
--- a/dart/cases/coyote.py
+++ b/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
diff --git a/dart/cases/n0012.py b/dart/cases/n0012.py
index ed9f3dae1a3eca1532ab7c2b236c137d79d8173f..474b4516d82ecf7121462f0e23674b6aadbf0064 100644
--- a/dart/cases/n0012.py
+++ b/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
diff --git a/dart/cases/n64A410.py b/dart/cases/n64A410.py
index d033dcbdbec156fd160b65a4623c30d5c1ca1a53..b6a032343773fbbbe83e4630911d67c57505e12a 100644
--- a/dart/cases/n64A410.py
+++ b/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
diff --git a/dart/cases/rae2822.py b/dart/cases/rae2822.py
index 1b9c6cb4c07e7b9889e62c32c2365f655e930d1b..b18b12ea6cc7469c50a8245e09438f38186a0c62 100644
--- a/dart/cases/rae2822.py
+++ b/dart/cases/rae2822.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/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
diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py
index 73a8f1daaa682aca82263d5d7027baf1f2ecb3d8..0e6a5b00cfaed89cf70613da7d981863a7fe2bee 100644
--- a/dart/cases/wbht.py
+++ b/dart/cases/wbht.py
@@ -44,7 +44,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/wbht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
diff --git a/dart/cases/wht.py b/dart/cases/wht.py
index dd0fc7d10496d2bdafe0eacdc246bbcf7445c00e..1f216c56a6b8d90fc741d41b401c4474106da227 100644
--- a/dart/cases/wht.py
+++ b/dart/cases/wht.py
@@ -39,7 +39,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/wht.geo') # Input file containing the model
p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
diff --git a/dart/default.py b/dart/default.py
index 878ed741783c1c5dbea0eafad99b569a946e5027..c87cc94e64d8ab919264f7914284ce5b8c8eb2aa 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,7 +23,7 @@ from fwk.wutils import parseargs
import dart
import tbox
-def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
+def mesh(dim, file, pars, bnds, wk = 'wake', wktp = None):
"""Initialize mesh and mesh writer
Parameters
@@ -34,7 +34,7 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
file contaning mesh (w/o extention)
pars : dict
parameters for mesh
- bnd : array of str
+ bnds : array of str
names of crack (wake) boundaries physical tag
wk : str (optional)
name of crack (wake) physical tag
@@ -48,7 +48,8 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
# crack the mesh
mshCrck = tbox.MshCrack(msh, dim)
mshCrck.setCrack(wk)
- mshCrck.addBoundaries(bnd)
+ for bnd in bnds:
+ mshCrck.addBoundary(bnd)
if dim == 3 and wktp is not None:
mshCrck.setExcluded(wktp)
mshCrck.run()
@@ -151,7 +152,7 @@ def picard(pbl):
problem formulation
"""
args = parseargs()
- solver = dart.Picard(pbl, tbox.Gmres(1, 1e-6, 30, 1e-8), nthrds=args.k, vrb=args.verb+1)
+ solver = dart.Picard(pbl, tbox.Gmres(1, 1e-6, 30, 1e-8), nthrds=args.k, vrb=args.verb+2)
return solver
def newton(pbl):
@@ -164,7 +165,7 @@ def newton(pbl):
"""
from tbox.solvers import LinearSolver
args = parseargs()
- solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+1)
+ solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+2)
return solver
def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
@@ -190,11 +191,12 @@ def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = '
args = parseargs()
mshDef = tbox.MshDeform(msh, tbox.Gmres(1, 1e-6, 30, 1e-8), dim, nthrds=args.k)
mshDef.setField(fld)
- mshDef.addFixed(fxd)
+ for tag in fxd:
+ mshDef.addFixed(tag)
mshDef.addMoving(mov)
if dim == 3:
mshDef.setSymmetry(sym, 1)
- mshDef.addInternal([wk, wk+'_'])
+ mshDef.addInternal(wk, wk+'_')
mshDef.initialize()
return mshDef
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 40b3b60b40afb7c74b3dbbc366e29bf1535737da..a404e630ee46fe79e5b79edce1f89932ce43240b 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -58,17 +58,23 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
nthreads = sol->nthreads;
verbose = sol->verbose;
- // Linear solvers
+ // Linear solvers (if GMRES, use the same, but with a thighter tolerance)
if (sol->linsol->type() == LSOLTYPE::GMRES_ILUT)
{
- // Use the same GMRES, but with a thighter tolerance
std::shared_ptr<Gmres> gmres = std::make_shared<Gmres>(*dynamic_cast<Gmres *>(sol->linsol.get()));
gmres->setTolerance(1e-8);
alinsol = gmres;
}
else
alinsol = sol->linsol;
- mlinsol = morph->linsol;
+ if (morph->linsol->type() == LSOLTYPE::GMRES_ILUT)
+ {
+ std::shared_ptr<Gmres> gmres = std::make_shared<Gmres>(*dynamic_cast<Gmres *>(morph->linsol.get()));
+ gmres->setTolerance(1e-12);
+ mlinsol = gmres;
+ }
+ else
+ mlinsol = morph->linsol;
// Update element gradients
sol->pbl->initGradElems();
@@ -101,6 +107,14 @@ void Adjoint::run()
// Init
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
+ // Display current freestream conditions
+ if (verbose > 0)
+ std::cout << std::fixed << std::setprecision(2)
+ << "Computing gradients for Mach " << sol->pbl->M_inf << ", "
+ << sol->pbl->alpha * 180 / 3.14159 << "deg AoA, "
+ << sol->pbl->beta * 180 / 3.14159 << "deg AoS"
+ << std::endl;
+
// Compute partial gradients of flow residuals and solve flow adjoint
tms["0-AdjFlo"].start();
this->linearizeFlow(); // dRu/dU, dRu/dA, dRu/dX
@@ -180,7 +194,7 @@ void Adjoint::run()
<< std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lamCdMsh.data(), lamCdMsh.size()).norm() << std::endl;
// Display timers
- if (verbose > 1)
+ if (verbose > 2)
std::cout << "Adjoint solver CPU" << std::endl
<< tms;
std::cout << std::endl;
@@ -415,7 +429,7 @@ void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor
/**
* @brief Build gradients of the load coefficients with respect to the flow variables on a given body
- *
+ *
* Note: this is essentially the same as buildGradientLoadsFlow, but with an additional projection performed inside the loop for performance
*/
void Adjoint::buildGradientCoefficientsFlow(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd)
@@ -661,7 +675,7 @@ void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor
/**
* @brief Build gradients of the loads with respect to mesh coordinates on a given body
- *
+ *
* Note: this is essentially the same as buildGradientLoadsMesh, but with an additional projection performed inside the loop for performance
*/
void Adjoint::buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd)
@@ -716,10 +730,15 @@ void Adjoint::buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowV
/**
* @brief Write the results
*/
-void Adjoint::save(MshExport *mshWriter, int n)
+void Adjoint::save(MshExport *mshWriter, std::string const &suffix)
{
// Write files
- std::cout << "Saving files... " << std::endl;
+ std::cout << "Saving files "
+ << std::setprecision(2)
+ << "(Mach " << sol->pbl->M_inf << ", "
+ << sol->pbl->alpha * 180 / 3.14159 << " deg AoA, "
+ << sol->pbl->beta * 180 / 3.14159 << " deg AoS)"
+ << std::endl;
// setup results
Results results;
results.scalars_at_nodes["lambdaClPhi"] = &lamClFlo;
@@ -727,18 +746,9 @@ void Adjoint::save(MshExport *mshWriter, int n)
results.vectors_at_nodes["gradClMsh"] = &dClMsh;
results.vectors_at_nodes["gradCdMsh"] = &dCdMsh;
// save (whole mesh and bodies)
- if (n > 0)
- {
- mshWriter->save(sol->pbl->msh->name + "_adjoint_" + std::to_string(n), results);
- for (auto bnd : sol->pbl->bodies)
- bnd->save(bnd->groups[0]->tag->name + "adjoint_" + std::to_string(n), results);
- }
- else
- {
- mshWriter->save(sol->pbl->msh->name + "_adjoint", results);
- for (auto bnd : sol->pbl->bodies)
- bnd->save(bnd->groups[0]->tag->name + "_adjoint", results);
- }
+ mshWriter->save(sol->pbl->msh->name + "_adjoint" + suffix, results);
+ for (auto bnd : sol->pbl->bodies)
+ bnd->save(sol->pbl->msh->name + '_' + bnd->groups[0]->tag->name + "_adjoint" + suffix, results);
}
void Adjoint::write(std::ostream &out) const
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index a33efb11100d5775f064e882e91b57d56ea63de0..b2dc8ffd116cf3415e30bde0fef648e4fac4fb54 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -66,7 +66,7 @@ public:
virtual ~Adjoint() { std::cout << "~Adjoint()\n"; }
void run();
- void save(tbox::MshExport *mshWriter, int n = 0);
+ void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
// Partial gradients of the mesh residuals
void linearizeMesh();
diff --git a/dart/src/wF1Ct.h b/dart/src/wF1Ct.h
index 010791b123483047afe138fbb9621992ed61baa5..3041a2fbf7edc24c2185044fc70d9a8ddf4c1e09 100644
--- a/dart/src/wF1Ct.h
+++ b/dart/src/wF1Ct.h
@@ -50,7 +50,7 @@ class DART_API F1CtDrag : public F1Ct
public:
F1CtDrag(int _nDim, double _sRef, double alpha, double _beta = 0);
- void update(double alpha);
+ virtual void update(double alpha) override;
virtual Eigen::Vector3d eval() const override;
virtual Eigen::Vector3d evalGrad() const override;
@@ -69,7 +69,7 @@ class DART_API F1CtSide : public F1Ct
public:
F1CtSide(int _nDim, double _sRef, double alpha, double _beta = 0);
- void update(double alpha);
+ virtual void update(double alpha) override;
virtual Eigen::Vector3d eval() const override;
virtual Eigen::Vector3d evalGrad() const override;
@@ -88,7 +88,7 @@ class DART_API F1CtLift : public F1Ct
public:
F1CtLift(int _nDim, double _sRef, double alpha, double _beta = 0);
- void update(double alpha);
+ virtual void update(double alpha) override;
virtual Eigen::Vector3d eval() const override;
virtual Eigen::Vector3d evalGrad() const override;
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index 23ab5bf6689496c0918338882c5fb92d497ea847..1d53b14ba43af71ade78170de7c4fd284822956f 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -88,14 +88,12 @@ STATUS Newton::run()
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Display current freestream conditions
- if (printOn)
- {
- std::cout << std::setprecision(2)
- << "- Mach " << pbl->M_inf << ", "
- << pbl->alpha * 180 / 3.14159 << " deg AoA, "
- << pbl->beta * 180 / 3.14159 << " deg AoS"
+ if (verbose > 0)
+ std::cout << std::fixed << std::setprecision(2)
+ << "Solving flow for Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << "deg AoA, "
+ << pbl->beta * 180 / 3.14159 << "deg AoS"
<< std::endl;
- }
// Initialize solver loop
nIt = 0; // iteration counter
@@ -121,7 +119,7 @@ STATUS Newton::run()
ls.set(maxLsIt, lsTol, verbose);
// Display residual
- if (printOn)
+ if (verbose > 0)
{
std::cout << std::setw(8) << "N_Iter"
<< std::setw(8) << "L_Iter"
@@ -204,7 +202,7 @@ STATUS Newton::run()
Solver::computeLoad();
// Display residual (at each iteration)
- if (verbose > 0 && printOn)
+ if (verbose > 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
@@ -226,7 +224,7 @@ STATUS Newton::run()
} while (nIt < maxIt);
// Display residual (only last iteration)
- if (verbose == 0 && printOn)
+ if (verbose == 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
@@ -244,29 +242,35 @@ STATUS Newton::run()
// Check the solution
if (relRes < relTol || absRes < absTol)
{
- if (printOn) {std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;}
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::CONVERGED;
}
else if (std::isnan(relRes))
{
- std::cout << ANSI_COLOR_RED << "Newton solver diverged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_RED << "Newton solver diverged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::FAILED;
}
else
{
- if (printOn) {std::cout << ANSI_COLOR_YELLOW << "Newton solver not fully converged!" << ANSI_COLOR_RESET << std::endl;}
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_YELLOW << "Newton solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Newton solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::MAXIT;
}
}
@@ -400,7 +404,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
J.prune(0.);
J.makeCompressed();
- if (verbose > 2)
+ if (verbose > 3)
std::cout << "J (" << J.rows() << "," << J.cols() << ") nnz=" << J.nonZeros() << "\n";
}
@@ -510,7 +514,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
for (auto nod : dBC->nodes)
R(nod->row) = 0.;
- if (verbose > 2)
+ if (verbose > 3)
std::cout << "R (" << R.size() << ")\n";
}
diff --git a/dart/src/wPicard.cpp b/dart/src/wPicard.cpp
index 5d567514923b835699a796af894a5b5b32ace767..15177ebcb14e46fec7581fafc6daf94ae00b2bd7 100644
--- a/dart/src/wPicard.cpp
+++ b/dart/src/wPicard.cpp
@@ -81,11 +81,12 @@ STATUS Picard::run()
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Display current freestream conditions
- std::cout << std::setprecision(2)
- << "- Mach " << pbl->M_inf << ", "
- << pbl->alpha * 180 / 3.14159 << " deg AoA, "
- << pbl->beta * 180 / 3.14159 << " deg AoS"
- << std::endl;
+ if (verbose > 0)
+ std::cout << std::fixed << std::setprecision(2)
+ << "Solving flow for Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << "deg AoA, "
+ << pbl->beta * 180 / 3.14159 << "deg AoS"
+ << std::endl;
// Initialize solver loop
nIt = 0;
@@ -95,11 +96,12 @@ STATUS Picard::run()
Eigen::Map<Eigen::VectorXd> phi_(phi.data(), phi.size()), rPhi_(rPhi.data(), rPhi.size());
Eigen::VectorXd phiOld(phi.size());
- std::cout << std::setw(8) << "N_Iter"
- << std::setw(12) << "Cl"
- << std::setw(12) << "Cd"
- << std::setw(15) << "Rel_Res[phi]"
- << std::setw(15) << "Abs_Res[phi]" << std::endl;
+ if (verbose > 0)
+ std::cout << std::setw(8) << "N_Iter"
+ << std::setw(12) << "Cl"
+ << std::setw(12) << "Cd"
+ << std::setw(15) << "Rel_Res[phi]"
+ << std::setw(15) << "Abs_Res[phi]" << std::endl;
do
{
@@ -122,7 +124,7 @@ STATUS Picard::run()
Solver::computeLoad();
// Display residual (at each iteration)
- if (verbose > 0 || nIt == 0)
+ if (verbose > 1 || (verbose == 1 && nIt == 0))
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
@@ -152,7 +154,7 @@ STATUS Picard::run()
} while (nIt < maxIt);
// Display residual (only last iteration)
- if (verbose == 0)
+ if (verbose == 1)
{
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << nIt
@@ -168,29 +170,35 @@ STATUS Picard::run()
// Check the solution
if (relRes < relTol || absRes < absTol)
{
- std::cout << ANSI_COLOR_GREEN << "Picard solver converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_GREEN << "Picard solver converged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::CONVERGED;
}
else if (std::isnan(relRes))
{
- std::cout << ANSI_COLOR_RED << "Picard sovler diverged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_RED << "Picard sovler diverged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::FAILED;
}
else
{
- std::cout << ANSI_COLOR_YELLOW << "Picard solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 1)
+ if (verbose > 0)
+ std::cout << ANSI_COLOR_YELLOW << "Picard solver not fully converged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 2)
std::cout << "Picard solver CPU" << std::endl
<< tms;
- std::cout << std::endl;
+ if (verbose > 0)
+ std::cout << std::endl;
return STATUS::MAXIT;
}
}
@@ -345,7 +353,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
A.prune(0.);
A.makeCompressed();
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "A (" << A.rows() << "," << A.cols() << ") nnz=" << A.nonZeros() << "\n";
std::cout << "b (" << b.size() << ")\n";
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 68d2b6119cc09b811844dcab4e8a681db8149d2f..32e8c70175ff8009f1450cf65714685ab0d8902c 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -161,10 +161,15 @@ STATUS Solver::run()
/**
* @brief Write the results
*/
-void Solver::save(MshExport *mshWriter, int n)
+void Solver::save(MshExport *mshWriter, std::string const &suffix)
{
// Write files
- std::cout << "Saving files... " << std::endl;
+ std::cout << "Saving files "
+ << std::fixed << std::setprecision(2)
+ << "(Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << "deg AoA, "
+ << pbl->beta * 180 / 3.14159 << "deg AoS)"
+ << std::endl;
// setup results
Results results;
results.scalars_at_nodes["phi"] = φ
@@ -175,19 +180,9 @@ void Solver::save(MshExport *mshWriter, int n)
results.scalars_at_nodes["Mach"] = &M;
results.scalars_at_nodes["Cp"] = &Cp;
// save (whole mesh and bodies)
- if (n > 0)
- {
- mshWriter->save(pbl->msh->name + "_" + std::to_string(n), results);
- for (auto bnd : pbl->bodies)
- bnd->save(bnd->groups[0]->tag->name + "_" + std::to_string(n), results);
- }
- else
- {
- mshWriter->save(pbl->msh->name, results);
- for (auto bnd : pbl->bodies)
- bnd->save(bnd->groups[0]->tag->name, results);
- }
- std::cout << std::endl;
+ mshWriter->save(pbl->msh->name + suffix, results);
+ for (auto bnd : pbl->bodies)
+ bnd->save(pbl->msh->name + '_' + bnd->groups[0]->tag->name + suffix, results);
}
/**
@@ -279,6 +274,6 @@ void Solver::computeFlow()
});
// Check maximum Mach number
- if (*std::max_element(M.begin(), M.end()) >= 1.25)
+ if (*std::max_element(M.begin(), M.end()) >= 1.25 && verbose > 0)
std::cout << ANSI_COLOR_YELLOW << "Max. Mach greater than 1.25!" << ANSI_COLOR_RESET << std::endl;
}
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index aab40ad17beb65aaee05537c73ab5bad693dc710..fdf41515f3278f2056eb7bf977f21dc0ecb18e10 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -78,7 +78,7 @@ public:
void reset();
virtual STATUS run();
- void save(tbox::MshExport *mshWriter, int n = 0);
+ void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
protected:
void computeFlow();
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 5788fe27441057128a222cd330e4d421add55673..df03af2e7fec29ae787f34edacad6ee37b71cc28 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import numpy as np
def main():
# timer
@@ -50,7 +50,7 @@ def main():
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
- morpher = floD.morpher(msh, dim, ['airfoil'])
+ morpher = floD.morpher(msh, dim, 'airfoil')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
diff --git a/dart/tests/bli_Sep.py b/dart/tests/bli_Sep.py
index 363809b55869ba500b58863197d547854a3def0f..cd35876dc54c6ecb9af7e497c4167792e9728dc2 100644
--- a/dart/tests/bli_Sep.py
+++ b/dart/tests/bli_Sep.py
@@ -46,6 +46,7 @@ import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import math
import math
diff --git a/dart/tests/cylinder.py b/dart/tests/cylinder.py
index 07d43117a1d25ebf38f3f0e24cdbb843fc383762..6a10c2b432dc9da509ccae0a388b300429fd1c30 100644
--- a/dart/tests/cylinder.py
+++ b/dart/tests/cylinder.py
@@ -25,11 +25,11 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import math
import dart.default as floD
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import math
def main():
# timer
diff --git a/dart/tests/lift.py b/dart/tests/lift.py
index daa9b6eb2ed31cdfbddb45f0e9b5448b7eb46210..76dea0dbbefff88d58075e3ef701866e3363e2d0 100644
--- a/dart/tests/lift.py
+++ b/dart/tests/lift.py
@@ -25,14 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import math
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
-from matplotlib import pyplot as plt
+import math
def main():
# timer
diff --git a/dart/tests/morpher.py b/dart/tests/morpher.py
index c82f27719227c722d27d4db4d80d0d17d546d4d3..5bca5c11e4ebff7b62457edeb725c770d69267ba 100644
--- a/dart/tests/morpher.py
+++ b/dart/tests/morpher.py
@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import numpy as np
def main():
# timer
@@ -55,7 +55,7 @@ def main():
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
- morpher = floD.morpher(msh, dim, ['airfoil'])
+ morpher = floD.morpher(msh, dim, 'airfoil')
# mesh the reference geometry
pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01, 'angle' : alfa, 'xRot' : xc}
msh_ref, gmshWriter_ref = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
diff --git a/dart/tests/nonlift.py b/dart/tests/nonlift.py
index 223aa336b3ccd995c20e0bfa83e1c344e00ddded..4d6e911a825ddffc0321febeaee4fc4c69b62bbb 100644
--- a/dart/tests/nonlift.py
+++ b/dart/tests/nonlift.py
@@ -26,7 +26,6 @@
# Mesh refinement may have to be performed to obtain physical results.
# The residual might not fully converge if this test is used with gmsh4
-import math
import dart.utils as floU
import dart.default as floD
from matplotlib import pyplot as plt
@@ -34,6 +33,7 @@ import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import math
def main():
# timer
diff --git a/dart/tests/rae_25.py b/dart/tests/rae_25.py
index cc656c11f7211fad86603d24ea9b253615b57698..95ab324502dc9d1614f0822c89e69e3e4e92001a 100644
--- a/dart/tests/rae_25.py
+++ b/dart/tests/rae_25.py
@@ -2,13 +2,13 @@
# -*- 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.
@@ -25,11 +25,11 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import numpy as np
import dart.default as floD
import fwk
from fwk.testing import *
-from fwk.coloring import ccolors
+from fwk.coloring import ccolors
+import numpy as np
def main():
# timer
@@ -49,14 +49,14 @@ def main():
c_ref = 1.0 # reference length
S_ref = spn # reference area
fms = 1.0 # farfield mesh size
- nms = 0.01 # nearfield mesh size
+ nms = 0.01 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms["msh"].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
msh, gmshWriter = floD.mesh(dim, 'models/rae_25.geo', pars, ['field', 'wing', 'symmetry', 'downstream'])
- morpher = floD.morpher(msh, dim, ['wing'])
+ morpher = floD.morpher(msh, dim, 'wing')
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
@@ -122,11 +122,11 @@ def main():
tests.add(CTest('CS', solver.Cs, 0.0000, 1e-3, forceabs=True))
tests.add(CTest('CM', solver.Cm, -0.113, 5e-2)) # 2D -0.117
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 12.7, 5e-2)) # 2D 12.7, FD 12.73 (1e-6)
- tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.059, 5e-3, forceabs=True)) # 2D 0.051, FD 0.0557 (1e-6)
+ tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.055, 1e-2, forceabs=True)) # 2D 0.051, FD 0.0557 (1e-6)
tests.add(CTest('dCl_dMsh', dClX, 22.0, 5e-2))
- tests.add(CTest('dCd_dMsh', dCdX, 1.0, 5e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 1.0, 1e-1))
tests.add(CTest('dCl/dAoA (msh)', dClAoA, 12.9, 5e-2)) # 2D 12.9
- tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.058, 5e-3, forceabs=True)) # 2D 0.051
+ tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.055, 1e-2, forceabs=True)) # 2D 0.051
tests.run()
# eof
diff --git a/dart/tests/rae_3.py b/dart/tests/rae_3.py
index f1f8add3567e2643f5902b94731395bcd24b6f0c..6559ca1b1fa2d0d573213c17211a9c416d6e27b8 100644
--- a/dart/tests/rae_3.py
+++ b/dart/tests/rae_3.py
@@ -2,13 +2,13 @@
# -*- 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.
@@ -25,13 +25,13 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import numpy as np
import dart.default as floD
import dart.utils as floU
from matplotlib import pyplot as plt
import fwk
from fwk.testing import *
-from fwk.coloring import ccolors
+from fwk.coloring import ccolors
+import numpy as np
def main():
# timer
@@ -51,7 +51,7 @@ def main():
c_ref = 1.0 # reference length
S_ref = spn # reference area
fms = 1.0 # farfield mesh size
- nms = 0.02 # nearfield mesh size
+ nms = 0.02 # nearfield mesh size
# parameters for mesh deformation
alfa = 1*np.pi/180
@@ -62,8 +62,8 @@ def main():
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms["msh"].start()
pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms}
- msh, sortElm, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
- morpher = floD.morpher(msh, dim, ['wing'])
+ msh, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ morpher = floD.morpher(msh, dim, 'wing')
tms["msh"].stop()
# set the problem and add fluid, initial/boundary conditions
diff --git a/dart/utils.py b/dart/utils.py
index 9691f77b6a302af2dae89a2e4b80398af2bd8f63..03b98292724d1f3c9ec3fe13b5ddec30908b6045 100644
--- a/dart/utils.py
+++ b/dart/utils.py
@@ -111,7 +111,7 @@ def extract(vElems, vData):
i += 1
return data
-def writeSlices(mshName, ys, wId, n = 0):
+def writeSlices(mshName, ys, wId, sfx=''):
"""Write slice data for each ys coordinate along the span (only works with VTK)
Parameters
@@ -122,8 +122,8 @@ def writeSlices(mshName, ys, wId, n = 0):
y-coordinates of slices
wId : int
index of physical tag to slice (see gmsh)
- n : int (optional, > 0)
- index of mesh/solution file (will be appended to mesh file name)
+ sfx : string (optional)
+ suffix that will be appended to mesh/solution file name
"""
import numpy as np
try:
@@ -134,10 +134,7 @@ def writeSlices(mshName, ys, wId, n = 0):
import tboxVtk.reader as vtkR
import tboxVtk.cutter as vtkC
reader = vtkR.Reader()
- if n == 0:
- reader.open(mshName)
- else:
- reader.open(mshName+'_'+str(n))
+ reader.open(mshName+sfx)
cutter = vtkC.Cutter(reader.reader.GetOutputPort())
for i in range(0, len(ys)):
pts, elems, vals = cutter.extract(cutter.cut(wId, [0., ys[i], 0.], [0., 1., 0.]), 2, ['Cp'])
@@ -147,7 +144,4 @@ def writeSlices(mshName, ys, wId, n = 0):
data = np.hstack((data, vals['Cp']))
tU.sort(elems, data)
data = np.vstack((data, data[0,:]))
- if n == 0:
- tU.write(data, 'slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
- else:
- tU.write(data, 'slice_'+str(n)+'_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
+ tU.write(data, mshName+sfx+'_slice_'+str(i)+'.dat', '%1.5e', ', ', 'x, y, z, x/c, Cp', '')
diff --git a/dart/validation/agard.py b/dart/validation/agard.py
index 86d5b9b767a88ab9107adba31809b85365a4d772..dc18dcb7da37bdbb618c3cf05c12763a1338a6f5 100644
--- a/dart/validation/agard.py
+++ b/dart/validation/agard.py
@@ -19,13 +19,13 @@
## Compute flow around the Agard 445 wing at 1 degrees AOA and Mach 0.80
# Adrien Crovato
-import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tbxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import numpy as np
try:
import tboxVtk
@@ -84,7 +84,7 @@ def main():
tms['post'].start()
if canPost:
floU.writeSlices(msh.name,[0.01, 0.37, 0.75],5)
- data = tbxU.read('slice_1.dat')
+ data = tbxU.read('agard445_slice_1.dat')
tbxU.plot(data[:,3], data[:,4], 'x', 'Cp', 'Pressure distribution at MAC', True)
tms['post'].stop()
diff --git a/dart/validation/onera.py b/dart/validation/onera.py
index 6957094a6f26428f4bc3814c8ca757320c9cc1eb..85f4e671a9b9e8164e463fdefa2fe6e7a69f5d5a 100644
--- a/dart/validation/onera.py
+++ b/dart/validation/onera.py
@@ -19,13 +19,13 @@
## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84
# Adrien Crovato
-import numpy as np
import dart.utils as floU
import dart.default as floD
import tbox.utils as tbxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import numpy as np
try:
import tboxVtk
@@ -82,7 +82,7 @@ def main():
tms['post'].start()
if canPost:
floU.writeSlices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184],5)
- data = tbxU.read('slice_2.dat')
+ data = tbxU.read('oneraM6_slice_2.dat')
tbxU.plot(data[:,3], data[:,4], 'x', 'Cp', 'Pressure distribution at MAC', True)
tms['post'].stop()
diff --git a/ext/amfe b/ext/amfe
index 8a20b0f157021fc286fa14bb9099fba81092da5c..e6eafbe78604a7a27c7ee21f098f4c1e3704b0dc 160000
--- a/ext/amfe
+++ b/ext/amfe
@@ -1 +1 @@
-Subproject commit 8a20b0f157021fc286fa14bb9099fba81092da5c
+Subproject commit e6eafbe78604a7a27c7ee21f098f4c1e3704b0dc