From 641fd0ca10d1c154616204fd2772d62badd4e6f9 Mon Sep 17 00:00:00 2001 From: Paul Dechamps <paul.dechamps@uliege.be> Date: Mon, 29 Aug 2022 11:24:05 +0000 Subject: [PATCH] Added vii to dart --- CMakeLists.txt | 1 + dart/CMakeLists.txt | 1 - dart/api/core.py | 502 ++++---- dart/api/internal/polar.py | 268 ++-- dart/api/internal/trim.py | 268 ++-- dart/api/mphys.py | 1146 ++++++++--------- dart/tests/bli.py | 129 -- dart/viscous/coupler.py | 82 -- dart/viscous/newton.py | 75 -- dart/viscous/solver.py | 489 ------- vii/CMakeLists.txt | 27 + vii/__init__.py | 21 + vii/_src/CMakeLists.txt | 35 + vii/_src/viiw.h | 17 + vii/_src/viiw.i | 57 + {dart/viscous => vii/api}/__init__.py | 0 vii/api/core.py | 110 ++ vii/coupler.py | 97 ++ vii/interfaces/dart/DartInterface.py | 167 +++ vii/interfaces/dart/__init__.py | 1 + .../interfaces/viscous/DAirfoil.py | 24 +- .../interfaces/viscous/DBoundary.py | 0 vii/interfaces/viscous/DGroupBuilder.py | 169 +++ .../interfaces/viscous/DWake.py | 16 +- vii/src/CMakeLists.txt | 27 + vii/src/DBoundaryLayer.cpp | 184 +++ vii/src/DBoundaryLayer.h | 68 + vii/src/DClosures.cpp | 256 ++++ vii/src/DClosures.h | 30 + vii/src/DDiscretization.cpp | 54 + vii/src/DDiscretization.h | 49 + vii/src/DDriver.cpp | 775 +++++++++++ vii/src/DDriver.h | 74 ++ vii/src/DEdge.cpp | 59 + vii/src/DEdge.h | 44 + vii/src/DFluxes.cpp | 242 ++++ vii/src/DFluxes.h | 31 + vii/src/DSolver.cpp | 247 ++++ vii/src/DSolver.h | 51 + vii/src/vii.h | 53 + vii/tests/bli2D.py | 123 ++ vii/utils.py | 143 ++ 42 files changed, 4347 insertions(+), 1865 deletions(-) delete mode 100644 dart/tests/bli.py delete mode 100644 dart/viscous/coupler.py delete mode 100644 dart/viscous/newton.py delete mode 100644 dart/viscous/solver.py create mode 100644 vii/CMakeLists.txt create mode 100644 vii/__init__.py create mode 100644 vii/_src/CMakeLists.txt create mode 100644 vii/_src/viiw.h create mode 100644 vii/_src/viiw.i rename {dart/viscous => vii/api}/__init__.py (100%) create mode 100644 vii/api/core.py create mode 100644 vii/coupler.py create mode 100644 vii/interfaces/dart/DartInterface.py create mode 100644 vii/interfaces/dart/__init__.py rename dart/viscous/airfoil.py => vii/interfaces/viscous/DAirfoil.py (93%) mode change 100644 => 100755 rename dart/viscous/boundary.py => vii/interfaces/viscous/DBoundary.py (100%) mode change 100644 => 100755 create mode 100755 vii/interfaces/viscous/DGroupBuilder.py rename dart/viscous/wake.py => vii/interfaces/viscous/DWake.py (90%) mode change 100644 => 100755 create mode 100644 vii/src/CMakeLists.txt create mode 100644 vii/src/DBoundaryLayer.cpp create mode 100644 vii/src/DBoundaryLayer.h create mode 100644 vii/src/DClosures.cpp create mode 100644 vii/src/DClosures.h create mode 100644 vii/src/DDiscretization.cpp create mode 100644 vii/src/DDiscretization.h create mode 100644 vii/src/DDriver.cpp create mode 100644 vii/src/DDriver.h create mode 100644 vii/src/DEdge.cpp create mode 100644 vii/src/DEdge.h create mode 100644 vii/src/DFluxes.cpp create mode 100644 vii/src/DFluxes.h create mode 100644 vii/src/DSolver.cpp create mode 100644 vii/src/DSolver.h create mode 100644 vii/src/vii.h create mode 100644 vii/tests/bli2D.py create mode 100644 vii/utils.py diff --git a/CMakeLists.txt b/CMakeLists.txt index ec130f6..f7e8db5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -114,6 +114,7 @@ ENDIF() # -- Sub directories ADD_SUBDIRECTORY( ext ) ADD_SUBDIRECTORY( dart ) +ADD_SUBDIRECTORY( vii ) # -- FINAL MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}") diff --git a/dart/CMakeLists.txt b/dart/CMakeLists.txt index 2d2803f..75515b2 100644 --- a/dart/CMakeLists.txt +++ b/dart/CMakeLists.txt @@ -24,5 +24,4 @@ INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py ${CMAKE_CURRENT_LIST_DIR}/utils.py DESTINATION dart) INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api - ${CMAKE_CURRENT_LIST_DIR}/viscous DESTINATION dart) diff --git a/dart/api/core.py b/dart/api/core.py index 8359e35..6518f84 100644 --- a/dart/api/core.py +++ b/dart/api/core.py @@ -1,236 +1,266 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Initialize DART -# Adrien Crovato - -def initDart(cfg, scenario='aerodynamic', task='analysis'): - """Initlialize DART for various API - - Parameters - ---------- - cfg : dict - Dictionary of parameters to configure DART - scenario : string, optional - Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic) - task : string, optional - Type of task (available: analysis, optimization) (default: analysis) - - Returns - ------- - _dim : int - Dimension of the problem - _qinf : float - Freestream dynamic pressure - _msh : tbox.MshData object - Mesh data structure - _wrtr : tbox.MshExport object - Mesh writer - _mrf : tbox.MshDeform object - Mesh morpher - _pbl : dart.Probem object - Problem definition - _bnd : Dart.Body object - Body of interest - _sol : dart.Newton object - Direct (Newton) solver - _adj : dart.Adjoint object - Adjoint solver - """ - # Imports - import dart - import tbox - import tbox.gmsh as gmsh - import math - - # Basic checks and config - # dimension - if cfg['Dim'] != 2 and cfg['Dim'] != 3: - raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n') - else: - _dim = cfg['Dim'] - # aoa and aos - if 'AoA' in cfg: - alpha = cfg['AoA'] * math.pi / 180 - else: - alpha = 0 - if _dim == 2: - if 'AoS' in cfg and cfg['AoS'] != 0: - raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n') - else: - beta = 0 - if 'Symmetry' in cfg: - raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n') - else: - if 'AoS' in cfg: - if cfg['AoS'] != 0 and 'Symmetry' in cfg: - raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n') - beta = cfg['AoS'] * math.pi / 180 - else: - beta = 0 - # save format - if cfg['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 cfg: - nthrd = cfg['Threads'] - else: - nthrd = 1 - # verbosity - if 'Verb' in cfg: - verb = cfg['Verb'] - else: - verb = 1 - # scenario and task type - if scenario != 'aerodynamic' and scenario != 'aerostructural': - raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n') - if task != 'analysis' and task != 'optimization': - raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n') - # dynamic pressure - if scenario == 'aerostructural': - _qinf = cfg['Q_inf'] - else: - _qinf = None - - # Mesh creation - _msh = gmsh.MeshLoader(cfg['File']).execute(**cfg['Pars']) - # add the wake - if _dim == 2: - mshCrck = tbox.MshCrack(_msh, _dim) - mshCrck.setCrack(cfg['Wake']) - mshCrck.addBoundary(cfg['Fluid']) - mshCrck.addBoundary(cfg['Farfield'][-1]) - mshCrck.addBoundary(cfg['Wing']) - mshCrck.run() - else: - for i in range(len(cfg['Wakes'])): - mshCrck = tbox.MshCrack(_msh, _dim) - mshCrck.setCrack(cfg['Wakes'][i]) - mshCrck.addBoundary(cfg['Fluid']) - mshCrck.addBoundary(cfg['Farfield'][-1]) - mshCrck.addBoundary(cfg['Wings'][i]) - if 'Fuselage' in cfg: - mshCrck.addBoundary(cfg['Fuselage']) - if 'Symmetry' in cfg: - mshCrck.addBoundary(cfg['Symmetry']) - if 'WakeTips' in cfg: - mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D) - mshCrck.run() - # save the updated mesh in native (gmsh) format and then set the writer to the requested format - tbox.GmshExport(_msh).save(_msh.name) - del mshCrck - _wrtr = Writer(_msh) - print('\n') - - # Mesh morpher creation - if scenario == 'aerostructural' or task == 'optimization': - _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd, vrb=verb) - _mrf.setField(cfg['Fluid']) - for tag in cfg['Farfield']: - _mrf.addFixed(tag) - if _dim == 2: - _mrf.addMoving(cfg['Wing']) - _mrf.addInternal(cfg['Wake'], cfg['Wake']+'_') - else: - if 'Fuselage' in cfg: - _mrf.addFixed(cfg['Fuselage']) - if 'Symmetry' in cfg: - _mrf.setSymmetry(cfg['Symmetry'], 1) - for i in range(len(cfg['Wings'])): - if i == 0: - _mrf.addMoving(cfg['Wings'][i]) # TODO body of interest (FSI/OPT) is always first given body - else: - _mrf.addFixed(cfg['Wings'][i]) - _mrf.addInternal(cfg['Wakes'][i], cfg['Wakes'][i]+'_') - _mrf.initialize() - else: - _mrf = None - - # Problem creation - _pbl = dart.Problem(_msh, _dim, alpha, beta, cfg['M_inf'], cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref']) - # add the field - _pbl.set(dart.Fluid(_msh, cfg['Fluid'], cfg['M_inf'], _dim, alpha, beta)) - # add the initial condition - _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta)) - # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically) - for i in range (len(cfg['Farfield'])): - _pbl.add(dart.Freestream(_msh, cfg['Farfield'][i], _dim, alpha, beta)) - # add solid boundaries - if _dim == 2: - bnd = dart.Body(_msh, [cfg['Wing'], cfg['Fluid']]) - _bnd = bnd - _pbl.add(bnd) - else: - _bnd = None - for bd in cfg['Wings']: - bnd = dart.Body(_msh, [bd, cfg['Fluid']]) - if _bnd is None: - _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body - _pbl.add(bnd) - if 'Fuselage' in cfg: - bnd = dart.Body(_msh, [cfg['Fuselage'], cfg['Fluid']]) - _pbl.add(bnd) - # add Wake/Kutta boundary conditions - # TODO refactor Wake "exclusion" for 3D? - if _dim == 2: - _pbl.add(dart.Wake(_msh, [cfg['Wake'], cfg['Wake']+'_', cfg['Fluid']])) - _pbl.add(dart.Kutta(_msh, [cfg['Te'], cfg['Wake']+'_', cfg['Wing'], cfg['Fluid']])) - else: - for i in range(len(cfg['Wakes'])): - if 'WakeExs' in cfg: - _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeExs'][i]])) # 3D + fuselage - elif 'WakeTips' in cfg: - _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeTips'][i]])) # 3D - else: - _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid']])) # 2.5D - _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']])) - - # Direct (forward) solver creation - # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones - # initialize the linear (inner) solver - if cfg['LSolver'] == 'PARDISO': - linsol = tbox.Pardiso() - elif cfg['LSolver'] == 'MUMPS': - linsol = tbox.Mumps() - elif cfg['LSolver'] == 'GMRES': - gfil = cfg['G_fill'] if 'G_fill' in cfg else 2 - grst = cfg['G_restart'] if 'G_restart' in cfg else 50 - gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5 - linsol = tbox.Gmres(gfil, 1e-6, grst, gtol) - else: - raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n') - # initialize the nonlinear (outer) solver - _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb) - - # Adjoint (reverse) solver creation - if task == 'optimization': - _adj = dart.Adjoint(_sol, _mrf) - else: - _adj = None - - # Return - return _dim, _qinf, _msh, _wrtr, _mrf, _pbl, _bnd, _sol, _adj +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +## Initialize DART +# Adrien Crovato + +def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False): + """Initlialize DART for various API + + Parameters + ---------- + cfg : dict + Dictionary of parameters to configure DART + scenario : string, optional + Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic) + task : string, optional + Type of task (available: analysis, optimization) (default: analysis) + viscous : bool, optional + Whether to configure DART for viscous-inviscid interaction (default: False) + + Returns + ------- + A dictionary containing + dim : int + Dimension of the problem + qinf : float + Freestream dynamic pressure + msh : tbox.MshData object + Mesh data structure + wrt : tbox.MshExport object + Mesh writer + mrf : tbox.MshDeform object + Mesh morpher + pbl : dart.Probem object + Problem definition + bnd : Dart.Body object + Body of interest + blwb : Dart.Blowing object + Transpiration B.C. on body + blww : Dart.Blowing object + Transpiration B.C. on wake + sol : dart.Newton object + Direct (Newton) solver + adj : dart.Adjoint object + Adjoint solver + """ + # Imports + import dart + import tbox + import tbox.gmsh as gmsh + import math + + # Basic checks and config + # dimension + if cfg['Dim'] != 2 and cfg['Dim'] != 3: + raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n') + else: + _dim = cfg['Dim'] + # aoa and aos + if 'AoA' in cfg: + alpha = cfg['AoA'] * math.pi / 180 + else: + alpha = 0 + if _dim == 2: + if 'AoS' in cfg and cfg['AoS'] != 0: + raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n') + else: + beta = 0 + if 'Symmetry' in cfg: + raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n') + else: + if 'AoS' in cfg: + if cfg['AoS'] != 0 and 'Symmetry' in cfg: + raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n') + beta = cfg['AoS'] * math.pi / 180 + else: + beta = 0 + # save format + if cfg['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 cfg: + nthrd = cfg['Threads'] + else: + nthrd = 1 + # verbosity + if 'Verb' in cfg: + verb = cfg['Verb'] + else: + verb = 1 + # scenario and task type + if scenario != 'aerodynamic' and scenario != 'aerostructural': + raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n') + if task != 'analysis' and task != 'optimization': + raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n') + # dynamic pressure + if scenario == 'aerostructural': + _qinf = cfg['Q_inf'] + else: + _qinf = None + + # Mesh creation + _msh = gmsh.MeshLoader(cfg['File']).execute(**cfg['Pars']) + # add the wake + if _dim == 2: + mshCrck = tbox.MshCrack(_msh, _dim) + mshCrck.setCrack(cfg['Wake']) + mshCrck.addBoundary(cfg['Fluid']) + mshCrck.addBoundary(cfg['Farfield'][-1]) + mshCrck.addBoundary(cfg['Wing']) + mshCrck.run() + else: + for i in range(len(cfg['Wakes'])): + mshCrck = tbox.MshCrack(_msh, _dim) + mshCrck.setCrack(cfg['Wakes'][i]) + mshCrck.addBoundary(cfg['Fluid']) + mshCrck.addBoundary(cfg['Farfield'][-1]) + mshCrck.addBoundary(cfg['Wings'][i]) + if 'Fuselage' in cfg: + mshCrck.addBoundary(cfg['Fuselage']) + if 'Symmetry' in cfg: + mshCrck.addBoundary(cfg['Symmetry']) + if 'WakeTips' in cfg: + mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D) + mshCrck.run() + # save the updated mesh in native (gmsh) format and then set the writer to the requested format + tbox.GmshExport(_msh).save(_msh.name) + del mshCrck + _wrt = Writer(_msh) + print('\n') + + # Mesh morpher creation + if scenario == 'aerostructural' or task == 'optimization': + _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd, vrb=verb) + _mrf.setField(cfg['Fluid']) + for tag in cfg['Farfield']: + _mrf.addFixed(tag) + if _dim == 2: + _mrf.addMoving(cfg['Wing']) + _mrf.addInternal(cfg['Wake'], cfg['Wake']+'_') + else: + if 'Fuselage' in cfg: + _mrf.addFixed(cfg['Fuselage']) + if 'Symmetry' in cfg: + _mrf.setSymmetry(cfg['Symmetry'], 1) + for i in range(len(cfg['Wings'])): + if i == 0: + _mrf.addMoving(cfg['Wings'][i]) # TODO body of interest (FSI/OPT) is always first given body + else: + _mrf.addFixed(cfg['Wings'][i]) + _mrf.addInternal(cfg['Wakes'][i], cfg['Wakes'][i]+'_') + _mrf.initialize() + else: + _mrf = None + + # Problem creation + _pbl = dart.Problem(_msh, _dim, alpha, beta, cfg['M_inf'], cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref']) + # add the field + _pbl.set(dart.Fluid(_msh, cfg['Fluid'], cfg['M_inf'], _dim, alpha, beta)) + # add the initial condition + _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta)) + # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically) + for i in range (len(cfg['Farfield'])): + _pbl.add(dart.Freestream(_msh, cfg['Farfield'][i], _dim, alpha, beta)) + # add solid boundaries + if _dim == 2: + bnd = dart.Body(_msh, [cfg['Wing'], cfg['Fluid']]) + _bnd = bnd + _pbl.add(bnd) + else: + _bnd = None + for bd in cfg['Wings']: + bnd = dart.Body(_msh, [bd, cfg['Fluid']]) + if _bnd is None: + _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body + _pbl.add(bnd) + if 'Fuselage' in cfg: + bnd = dart.Body(_msh, [cfg['Fuselage'], cfg['Fluid']]) + _pbl.add(bnd) + # add Wake/Kutta boundary conditions + # TODO refactor Wake "exclusion" for 3D? + if _dim == 2: + _pbl.add(dart.Wake(_msh, [cfg['Wake'], cfg['Wake']+'_', cfg['Fluid']])) + _pbl.add(dart.Kutta(_msh, [cfg['Te'], cfg['Wake']+'_', cfg['Wing'], cfg['Fluid']])) + else: + for i in range(len(cfg['Wakes'])): + if 'WakeExs' in cfg: + _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeExs'][i]])) # 3D + fuselage + elif 'WakeTips' in cfg: + _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeTips'][i]])) # 3D + else: + _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid']])) # 2.5D + _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']])) + # add transpiration (blowing) boundary conditions + if viscous: + if _dim != 2 or task != 'analysis': + raise RuntimeError('Viscous-inviscid interaction calculations are currently supported only for 2D analysis!\n') + _blwb = dart.Blowing(_msh, cfg['Wing']) + _blww = dart.Blowing(_msh, cfg['Wake']) + _pbl.add(_blwb) + _pbl.add(_blww) + else: + _blwb = None + _blww = None + + # Direct (forward) solver creation + # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones + # initialize the linear (inner) solver + if cfg['LSolver'] == 'PARDISO': + linsol = tbox.Pardiso() + elif cfg['LSolver'] == 'MUMPS': + linsol = tbox.Mumps() + elif cfg['LSolver'] == 'GMRES': + gfil = cfg['G_fill'] if 'G_fill' in cfg else 2 + grst = cfg['G_restart'] if 'G_restart' in cfg else 50 + gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5 + linsol = tbox.Gmres(gfil, 1e-6, grst, gtol) + else: + raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n') + # initialize the nonlinear (outer) solver + _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb) + + # Adjoint (reverse) solver creation + if task == 'optimization': + _adj = dart.Adjoint(_sol, _mrf) + else: + _adj = None + + # Return + return { + 'dim': _dim, + 'qinf': _qinf, + 'msh': _msh, + 'wrt': _wrt, + 'mrf': _mrf, + 'pbl': _pbl, + 'bnd': _bnd, + 'blwb': _blwb, + 'blww': _blww, + 'sol': _sol, + 'adj': _adj + } diff --git a/dart/api/internal/polar.py b/dart/api/internal/polar.py index f4355f3..682dc6e 100644 --- a/dart/api/internal/polar.py +++ b/dart/api/internal/polar.py @@ -1,131 +1,137 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Polar -# Adrien Crovato -# -# Compute aerodynamic load coefficients as a function of angle of attack - -import math -import dart.utils as dU -import tbox.utils as tU -import fwk -from fwk.coloring import ccolors -from ..core import initDart - -class Polar: - """Utility to compute the polar of a lifting surface - - Attributes - ---------- - tms : fwk.Timers object - dictionary of timers - dim : int - problem dimensions - format : str - output format type - slice : array of float - arrays of y-coordinates to perform slices (for 3D only) - tag : int - index of physical tag to slice (see gmsh) - mach : float - freestream Mach number - alphas : array of float - range of angles of attack to compute - Cls : array of float - lift coefficients for the different AoAs - Cds : array of float - drag coefficients for the different AoAs - Cms : array of float - moment coefficients for the different AoAs - """ - def __init__(self, p): - """Setup and configure - - Parameters - ---------- - p : dict - dictionary of parameters to configure DART - """ - # basic init - self.tms = fwk.Timers() - self.tms["total"].start() - self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis') - # save some parameters for later use - self.format = p['Format'] - try: - self.slice = p['Slice'] - self.tag = p['TagId'] - except: - self.slice = None - self.tag = None - self.mach = p['M_inf'] - self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step']) - - def run(self): - """Compute the polar - """ - # initialize the loop - self.Cls = [] - self.Cds = [] - self.Cms = [] - print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET) - for i in range(len(self.alphas)): - # define current angle of attack - alpha = self.alphas[i]*math.pi/180 - acs = '_{:04d}'.format(i) - # update problem and reset ICs to improve robustness in transonic cases - self.pbl.update(alpha) - self.sol.reset() - # run the solver and save the results - print(ccolors.ANSI_BLUE + 'Running the solver for', self.alphas[i], 'degrees', ccolors.ANSI_RESET) - self.tms["solver"].start() - self.sol.run() - self.tms["solver"].stop() - self.sol.save(self.wrtr, acs) - # extract Cp - if self.dim == 2: - Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp) - tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.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, acs) - # extract force coefficients - self.Cls.append(self.sol.Cl) - self.Cds.append(self.sol.Cd) - self.Cms.append(self.sol.Cm) - - def disp(self): - """Display the results and draw the polar - """ - # display results - print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET) - print(" M alpha Cl Cd Cm") - i = 0 - while i < len(self.alphas): - print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i])) - i = i+1 - print('\n') - # display timers - self.tms["total"].stop() - print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET) - print(self.tms) - # plot results - if self.alphas[0] != self.alphas[-1]: - tU.plot(self.alphas, self.Cls, 'alpha', 'Cl', '') - tU.plot(self.alphas, self.Cds, 'alpha', 'Cd', '') - tU.plot(self.alphas, self.Cms, 'alpha', 'Cm', '') +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +## Polar +# Adrien Crovato +# +# Compute aerodynamic load coefficients as a function of angle of attack + +import math +import dart.utils as dU +import tbox.utils as tU +import fwk +from fwk.coloring import ccolors +from ..core import initDart + +class Polar: + """Utility to compute the polar of a lifting surface + + Attributes + ---------- + tms : fwk.Timers object + dictionary of timers + dim : int + problem dimensions + format : str + output format type + slice : array of float + arrays of y-coordinates to perform slices (for 3D only) + tag : int + index of physical tag to slice (see gmsh) + mach : float + freestream Mach number + alphas : array of float + range of angles of attack to compute + Cls : array of float + lift coefficients for the different AoAs + Cds : array of float + drag coefficients for the different AoAs + Cms : array of float + moment coefficients for the different AoAs + """ + def __init__(self, p): + """Setup and configure + + Parameters + ---------- + p : dict + dictionary of parameters to configure DART + """ + # basic init + self.tms = fwk.Timers() + self.tms["total"].start() + _dart = initDart(p, scenario='aerodynamic', task='analysis') + self.dim = _dart['dim'] + self.msh = _dart['msh'] + self.wrt = _dart['wrt'] + self.pbl = _dart['pbl'] + self.bnd = _dart['bnd'] + self.sol = _dart['sol'] + # save some parameters for later use + self.format = p['Format'] + try: + self.slice = p['Slice'] + self.tag = p['TagId'] + except: + self.slice = None + self.tag = None + self.mach = p['M_inf'] + self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step']) + + def run(self): + """Compute the polar + """ + # initialize the loop + self.Cls = [] + self.Cds = [] + self.Cms = [] + print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET) + for i in range(len(self.alphas)): + # define current angle of attack + alpha = self.alphas[i]*math.pi/180 + acs = '_{:04d}'.format(i) + # update problem and reset ICs to improve robustness in transonic cases + self.pbl.update(alpha) + self.sol.reset() + # run the solver and save the results + print(ccolors.ANSI_BLUE + 'Running the solver for', self.alphas[i], 'degrees', ccolors.ANSI_RESET) + self.tms["solver"].start() + self.sol.run() + self.tms["solver"].stop() + self.sol.save(self.wrt, acs) + # extract Cp + if self.dim == 2: + Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp) + tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.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, acs) + # extract force coefficients + self.Cls.append(self.sol.Cl) + self.Cds.append(self.sol.Cd) + self.Cms.append(self.sol.Cm) + + def disp(self): + """Display the results and draw the polar + """ + # display results + print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET) + print(" M alpha Cl Cd Cm") + i = 0 + while i < len(self.alphas): + print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i])) + i = i+1 + print('\n') + # display timers + self.tms["total"].stop() + print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET) + print(self.tms) + # plot results + if self.alphas[0] != self.alphas[-1]: + tU.plot(self.alphas, self.Cls, '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 b277d9f..fd25625 100644 --- a/dart/api/internal/trim.py +++ b/dart/api/internal/trim.py @@ -1,131 +1,137 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Trim analysis -# Adrien Crovato -# -# Find the angle of attack to match a specified lift coefficient - -import math -import dart.utils as dU -import tbox.utils as tU -import fwk -from fwk.coloring import ccolors -from ..core import initDart - -class Trim: - """Utility to perform the trim analysis of a given lifting configuration - Find the angle of attack to reach a desired lift coefficient - - Attributes - ---------- - tms : fwk.Timers object - dictionary of timers - dim : int - problem dimensions - format : str - output format type - slice : array of float - arrays of y-coordinates to perform slices (for 3D only) - tag : int - index of physical tag to slice (see gmsh) - mach : float - freestream Mach number - clTgt : float - target lift coefficient - alpha : float - initial angle of attack - dCl : float - initial (estimated) slope of lift curve - """ - def __init__(self, p): - """Setup and configure - - Parameters - ---------- - p : dict - dictionary of parameters to configure DART - """ - # basic init - self.tms = fwk.Timers() - self.tms["total"].start() - self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis') - # save some parameters for later use - self.format = p['Format'] - try: - self.slice = p['Slice'] - self.tag = p['TagId'] - except: - self.slice = None - self.tag = None - self.mach = p['M_inf'] - self.clTgt = p['CL'] - self.alpha = p['AoA']*math.pi/180 - self.dCl = p['dCL'] - - def run(self): - """Perform the trim analysis - """ - # initialize loop - it = 0 - while True: - # define current angle of attack - print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET) - # update problem and reset ICs to improve robustness in transonic cases - self.pbl.update(self.alpha) - self.sol.reset() - # run the solver - self.tms["solver"].start() - status = self.sol.run() - if status >= 2: - raise RuntimeError('Flow solver diverged!\n') - self.tms["solver"].stop() - # update slope - if it != 0: - self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_) - # break or store values and update AoA - if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl): - break - else: - Cl_ = self.sol.Cl - alpha_ = self.alpha - dAlpha = (self.clTgt - self.sol.Cl) / self.dCl - self.alpha += dAlpha - it += 1 - - # save results - self.sol.save(self.wrtr) - # extract Cp - if self.dim == 2: - Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp) - tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - elif self.dim == 3 and self.format == 'vtk' and self.slice: - dU.writeSlices(self.msh.name, self.slice, self.tag) - - def disp(self): - """Display the results - """ - # display results - print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET) - print(" M alpha dCl Cl Cd Cm") - print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm)) - print('\n') - # display timers - self.tms["total"].stop() - print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET) - print(self.tms) +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +## Trim analysis +# Adrien Crovato +# +# Find the angle of attack to match a specified lift coefficient + +import math +import dart.utils as dU +import tbox.utils as tU +import fwk +from fwk.coloring import ccolors +from ..core import initDart + +class Trim: + """Utility to perform the trim analysis of a given lifting configuration + Find the angle of attack to reach a desired lift coefficient + + Attributes + ---------- + tms : fwk.Timers object + dictionary of timers + dim : int + problem dimensions + format : str + output format type + slice : array of float + arrays of y-coordinates to perform slices (for 3D only) + tag : int + index of physical tag to slice (see gmsh) + mach : float + freestream Mach number + clTgt : float + target lift coefficient + alpha : float + initial angle of attack + dCl : float + initial (estimated) slope of lift curve + """ + def __init__(self, p): + """Setup and configure + + Parameters + ---------- + p : dict + dictionary of parameters to configure DART + """ + # basic init + self.tms = fwk.Timers() + self.tms["total"].start() + _dart = initDart(p, scenario='aerodynamic', task='analysis') + self.dim = _dart['dim'] + self.msh = _dart['msh'] + self.wrt = _dart['wrt'] + self.pbl = _dart['pbl'] + self.bnd = _dart['bnd'] + self.sol = _dart['sol'] + # save some parameters for later use + self.format = p['Format'] + try: + self.slice = p['Slice'] + self.tag = p['TagId'] + except: + self.slice = None + self.tag = None + self.mach = p['M_inf'] + self.clTgt = p['CL'] + self.alpha = p['AoA']*math.pi/180 + self.dCl = p['dCL'] + + def run(self): + """Perform the trim analysis + """ + # initialize loop + it = 0 + while True: + # define current angle of attack + print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET) + # update problem and reset ICs to improve robustness in transonic cases + self.pbl.update(self.alpha) + self.sol.reset() + # run the solver + self.tms["solver"].start() + status = self.sol.run() + if status >= 2: + raise RuntimeError('Flow solver diverged!\n') + self.tms["solver"].stop() + # update slope + if it != 0: + self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_) + # break or store values and update AoA + if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl): + break + else: + Cl_ = self.sol.Cl + alpha_ = self.alpha + dAlpha = (self.clTgt - self.sol.Cl) / self.dCl + self.alpha += dAlpha + it += 1 + + # save results + self.sol.save(self.wrt) + # extract Cp + if self.dim == 2: + Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp) + tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') + elif self.dim == 3 and self.format == 'vtk' and self.slice: + dU.writeSlices(self.msh.name, self.slice, self.tag) + + def disp(self): + """Display the results + """ + # display results + print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET) + print(" M alpha dCl Cl Cd Cm") + print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm)) + print('\n') + # display timers + self.tms["total"].stop() + print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET) + print(self.tms) diff --git a/dart/api/mphys.py b/dart/api/mphys.py index d23c81e..d61d932 100644 --- a/dart/api/mphys.py +++ b/dart/api/mphys.py @@ -1,573 +1,573 @@ -# -*- 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 - For compatibility, the output is distributed, though it should always be serial - - 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', distributed=True, val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates']) - - def 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 - For compatibility, the input is distributed, though it should always be serial - - 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', distributed=True, 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 - For compatibility, the input is distributed, though it should always be serial - """ - 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', distributed=True, 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) - self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged') - - 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) - status = self.sol.run() - 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 - - 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 - For compatibility, the output is distributed, though it should always be serial - - 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', distributed=True, 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 - - Attributes - ---------- - 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) - - 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_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 compute(self, inputs, outputs): - """Get the coefficients for the full aircraft - """ - # NB: inputs already up-to-date because DartSolver has already been run - # 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 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 coupling 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) - 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']) - 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'], 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']) - -# 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): - """Dart builder for MPHYS - - Attributes - ---------- - __dim : int - Dimension of the problem - __qinf : float - Freestream dynamic pressure - __wrt : tbox.MshExport object - Mesh and solution data writer - __mrf : tbox.MshDeform object - Mesh morpher - __sol : dart.Newton object - Direct Newton solver - __adj : dart.Adjoint object - Adjoint solver - """ - def __init__(self, cfg, scenario='aerodynamic', task='analysis', saveAll=False, raiseError=False): - """Instantiate the builder and save the parameters and options - - Parameters - ---------- - cfg : dict - Dictonnary of parameters to configure DART - scenario : string, optional - Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic) - task : string, optional - Type of task (available: analysis, optimization) (default: analysis) - 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 - """ - self.__cfg = cfg - self.__scn = scenario - self.__tsk = task - self.__sall = saveAll - self.__rerr = raiseError - - def initialize(self, comm): - """Instantiate and initialize all the solver components - - Parameters - ---------- - comm: mpi4py.MPI.Comm object - MPI communicator - """ - from .core import initDart - def _init(): - 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: - from mpi4py import MPI - # If several processes are provided by OpenMDAO, DART is requested to run with MPI, which is not supported - if comm.size > 1: - raise RuntimeError('DartBuiler.initialize - DART cannot be run using MPI, except for MultiPoint analysis!\n') - # Else, we are running a MultiPoint analysis and need to ID the actual process - true_comm = MPI.COMM_WORLD - if true_comm.size > 1: - if true_comm.rank != 0: - true_comm.recv(source=true_comm.rank-1) - _init() - if true_comm.rank != true_comm.size-1: - true_comm.send([], dest=true_comm.rank+1) - true_comm.barrier() - else: - _init() - except: - _init() - - 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) - - 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, rerr=self.__rerr) - - def get_post_coupling_subsystem(self, scenario_name=None): - """Return openMDAO group that computes the aero coefficients and writes data to disk - """ - 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 - """ - 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 + For compatibility, the output is distributed, though it should always be serial + + 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', distributed=True, val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates']) + + def 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 + For compatibility, the input is distributed, though it should always be serial + + 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', distributed=True, 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 + For compatibility, the input is distributed, though it should always be serial + """ + 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', distributed=True, 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) + self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged') + + 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) + status = self.sol.run() + 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 + + 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 + For compatibility, the output is distributed, though it should always be serial + + 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', distributed=True, 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 + + Attributes + ---------- + 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) + + 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_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 compute(self, inputs, outputs): + """Get the coefficients for the full aircraft + """ + # NB: inputs already up-to-date because DartSolver has already been run + # 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 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 coupling 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) + 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']) + 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'], 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']) + +# 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): + """Dart builder for MPHYS + + Attributes + ---------- + __dim : int + Dimension of the problem + __qinf : float + Freestream dynamic pressure + __wrt : tbox.MshExport object + Mesh and solution data writer + __mrf : tbox.MshDeform object + Mesh morpher + __sol : dart.Newton object + Direct Newton solver + __adj : dart.Adjoint object + Adjoint solver + """ + def __init__(self, cfg, scenario='aerodynamic', task='analysis', saveAll=False, raiseError=False): + """Instantiate the builder and save the parameters and options + + Parameters + ---------- + cfg : dict + Dictonnary of parameters to configure DART + scenario : string, optional + Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic) + task : string, optional + Type of task (available: analysis, optimization) (default: analysis) + 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 + """ + self.__cfg = cfg + self.__scn = scenario + self.__tsk = task + self.__sall = saveAll + self.__rerr = raiseError + + def initialize(self, comm): + """Instantiate and initialize all the solver components + + Parameters + ---------- + comm: mpi4py.MPI.Comm object + MPI communicator + """ + from .core import initDart + def _init(): + self.__dart = 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: + from mpi4py import MPI + # If several processes are provided by OpenMDAO, DART is requested to run with MPI, which is not supported + if comm.size > 1: + raise RuntimeError('DartBuiler.initialize - DART cannot be run using MPI, except for MultiPoint analysis!\n') + # Else, we are running a MultiPoint analysis and need to ID the actual process + true_comm = MPI.COMM_WORLD + if true_comm.size > 1: + if true_comm.rank != 0: + true_comm.recv(source=true_comm.rank-1) + _init() + if true_comm.rank != true_comm.size-1: + true_comm.send([], dest=true_comm.rank+1) + true_comm.barrier() + else: + _init() + except: + _init() + + def get_mesh_coordinate_subsystem(self, scenario_name=None): + """Return openMDAO component to get the initial surface mesh coordinates + """ + return DartMesh(dim=self.__dart['dim'], bnd=self.__dart['bnd']) + + def get_coupling_group_subsystem(self, scenario_name=None): + """Return openMDAO group containing the morpher, solver and loads + """ + return DartGroup(qinf=self.__dart['qinf'], bnd=self.__dart['bnd'], sol=self.__dart['sol'], mrf=self.__dart['mrf'], adj=self.__dart['adj'], rerr=self.__rerr) + + def get_post_coupling_subsystem(self, scenario_name=None): + """Return openMDAO group that computes the aero coefficients and writes data to disk + """ + return DartPostGroup(sol=self.__dart['sol'], adj=self.__dart['adj'], wrt=self.__dart['wrt'], snam=scenario_name, sall=self.__sall) + + def get_number_of_nodes(self): + """Return the number of surface nodes + """ + return len(self.__dart['bnd'].nodes) diff --git a/dart/tests/bli.py b/dart/tests/bli.py deleted file mode 100644 index 9fbf80b..0000000 --- a/dart/tests/bli.py +++ /dev/null @@ -1,129 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -import dart.viscous.solver as floVS -import dart.viscous.coupler as floVC -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors -import math - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 5*math.pi/180 - M_inf = 0.5 - c_ref = 1 - dim = 2 - tol = 2e-3 # tolerance for the VII (usually one drag count) #TODO tolerance has been decreased with new Kutta, hopefully will increase with Paul's dev - - # mesh the geometry - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - 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']) - tms['msh'].stop() - - # set the problem and add fluid, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - isolver = floD.newton(pbl) - vsolver = floVS.Solver(Re) - coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter) - coupler.run() - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - # visualize solution and plot results - floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58 - tests.add(CTest('Cd', vsolver.Cd, 0.0062, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.03, forceabs=True)) # Xfoil 0.056 - tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.03, forceabs=True)) - elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69 - tests.add(CTest('Cd', vsolver.Cd, 0.0067, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.03, forceabs=True)) # Xfoil 0.038 - tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.03, forceabs=True)) - elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39 - tests.add(CTest('Cd', vsolver.Cd, 0.011, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.03, forceabs=True)) # Xfoil 0.008 - tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.03, forceabs=True)) - else: - raise Exception('Test not defined for this flow') - - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/viscous/coupler.py b/dart/viscous/coupler.py deleted file mode 100644 index 6c339cf..0000000 --- a/dart/viscous/coupler.py +++ /dev/null @@ -1,82 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Viscous-inviscid coupler (quasi-simultaneous coupling) -# Amaury Bilocq - -import numpy as np -from fwk.coloring import ccolors -from dart.viscous.airfoil import Airfoil -from dart.viscous.wake import Wake - -class Coupler: - def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer): - self.isolver =_isolver # inviscid solver - self.vsolver = _vsolver # viscous solver - self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects - self.tol = _tol # tolerance of the VII - self.writer = _writer - - def run(self): - ''' Perform coupling - ''' - # initialize loop - it = 0 - converged = 0 # temp - CdOld = self.vsolver.Cd - while converged == 0: - print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET) - # run inviscid solver - self.isolver.run() - print('--- Viscous solver parameters ---') - print('Tolerance (drag count):', self.tol * 1e4) - print('--- Viscous problem definition ---') - print('Reynolds number:', self.vsolver.Re) - print('') - for n in range(0, len(self.group)): - print('Computing for', self.group[n].name, '...', end = ' ') - for i in range(0, len(self.group[n].boundary.nodes)): - self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0] - self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1] - self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2] - self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row] - self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row] - # run viscous solver - self.vsolver.run(self.group[n]) - if self.group[n].name == 'airfoil': - self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1] - self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0 - self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0 - # get blowing velocity from viscous solver and update inviscid solver - for i in range(0, self.group[n].nE): - self.group[n].boundary.setU(i, self.group[n].u[i]) - print('done!') - print(' Iter Cd xtr_top xtr_bot Error') - print('{0:8d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd)) - print('') - # Converged or not - if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol: - converged = 1 - else: - CdOld = self.vsolver.Cd - it += 1 - self.vsolver.it += 1 - # save results - self.isolver.save(self.writer) - self.vsolver.writeFile() - print('\n') diff --git a/dart/viscous/newton.py b/dart/viscous/newton.py deleted file mode 100644 index 3715fa8..0000000 --- a/dart/viscous/newton.py +++ /dev/null @@ -1,75 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Newton raphson method coupled with linear solver -# Amaury Bilocq - -import numpy as np -from fwk.coloring import ccolors - -def newton(f, x, maxIt, tol=1.0e-8): - '''Newton procedure to solve the non linear set of boundary layer equations. - Boundary layer equations are parabolic in x -> downstream marching is used - An implicit marching model is applied (Crank Nicolson) -> matrix inversion is performed - ''' - # Numerical Jacobian - def jaco(f, x): - eps = 1.0e-8 - n = len(x) - jac = np.zeros((n,n)) - xp = x.copy() - xm = x.copy() - for j in range(n): - # perturb solution - xTmp = x[j] - xp[j] += eps - xm[j] -= eps - # compute jacobian - jac[:,j] = (f(xp) - f(xm)) / (2 * eps) - # reset - xp[j] = xTmp - xm[j] = xTmp - return jac - # Backtrack line search - def ls(f, x, dx): - a = 1.0 - ires = np.linalg.norm(f(x)) - while a > 1/64: - res = np.linalg.norm(f(x + a * dx)) - if res < ires: - break - else: - a /= 2 - return a - # Solve - for i in range(maxIt): - # solve - jac = jaco(f,x) - dx = np.linalg.solve(jac, -f(x)) - # line search - a = ls(f, x, dx) - x += a * dx - fx = f(x) - err = np.linalg.norm(fx) - # safeguard - x[0] = max(x[0],0.) - x[1] = max(x[1],1.0005) - if i != 0 and err < tol: - return x - # Raise error but return solution and error - raise RuntimeError(ccolors.ANSI_YELLOW + 'Newton solver - too many iterations!\n' + ccolors.ANSI_RESET, x, err) diff --git a/dart/viscous/solver.py b/dart/viscous/solver.py deleted file mode 100644 index 9b96228..0000000 --- a/dart/viscous/solver.py +++ /dev/null @@ -1,489 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Integral boundary layer equations viscous solver -# Amaury Bilocq -# -# TODO list -# - Improve user-input -# -* let the user define the transition location -# -* let the user define the turbulence intensity (which will affect the amplification factor) -# - Improve the physics (robustness) -# -* wake should start turbulent -# -* use empirical correlations to fix Cteq at transition -# -- split the panel where the transition happens instead of averaging -# -* understand why Ct becomes negative in turbulent solver and prevent it -# - Improve the numerics (robustness) -# -* test the convergence by monitoring the change in the interaction velocity (vt) between 2 VVI iterations -# -- implement a more sophisticated line search (3 points quadratic?) -# -* try to use underrelaxation -# -- compute the jacobian analytically -# -* use an inverse coupling procedure when laminar separation occurs -# - Increase range of validity -# -- split the inviscid and viscous mesh -# -- check transonic predictions -# -- implement quasi-2D method for 3D flows - -from dart.viscous.boundary import Boundary -from dart.viscous.wake import Wake -import dart.viscous.newton as Newton -import numpy as np -from fwk.coloring import ccolors - -class Solver: - def __init__(self, _Re): - '''Initialize the viscous solver - ''' - self.Re = _Re # Reynolds number - self.gamma = 1.4 # heat capacity of air - self.Cd = 0 # drag coefficient - self.it = 0 # viscous inviscid iteration for interaction law - self.xtr = [1, 1] # set transition at TE reference coordinates - - def dictionary(self): - '''Create dictionary with the coordinates and the boundary layer parameters - ''' - wData = {} - wData['x'] = self.data[:,0] - wData['y'] = self.data[:,1] - wData['z'] = self.data[:,2] - wData['x/c'] = self.data[:,0]/max(self.data[:,0]) - wData['Cd'] = self.blScalars[0] - wData['Cdf'] = self.blScalars[1] - wData['Cdp'] = self.blScalars[2] - wData['delta*'] = self.blVectors[0] - wData['theta'] = self.blVectors[1] - wData['H'] = self.blVectors[2] - wData['Hk'] = self.blVectors[3] - wData['H*'] = self.blVectors[4] - wData['H**'] = self.blVectors[5] - wData['Cf'] = self.blVectors[6] - wData['Cdissip'] = self.blVectors[7] - wData['Ctau'] = self.blVectors[8] - wData['delta'] = self.blVectors[9] - wData['xtrTop'] = self.xtr[0] - wData['xtrBot'] = self.xtr[1] - return wData - - def writeFile(self): - '''Write the results in airfoil_viscous.dat - ''' - wData = self.dictionary() - toW = ['delta*', 'H', 'Cf'] # list of variable to be written (should be a user input) - # Write - print('writing file: airfoil_viscous.dat...', end = '') - f = open('airfoil_viscous.dat', 'w+') - f.write('$Aerodynamic coefficients\n') - f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n') - f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'], wData['Cdf'], wData['xtrTop'], wData['xtrBot'])) - f.write('$Boundary layer variables\n') - f.write(' x y z x/c') - for s in toW: - f.write(' {0:>15s}'.format(s)) - f.write('\n') - for i in range(len(self.data[:,0])): - f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i])) - for s in toW: - f.write(' {0:15.6f}'.format(wData[s][i])) - f.write('\n') - f.close() - print('done!') - - def __getBLcoordinates(self, data): - '''Transform the reference coordinates into airfoil coordinates - ''' - nN = len(data[:,0]) - nE = nN-1 #nbr of element along the surface - vt = np.zeros(nN) - xx = np.zeros(nN) # position along the surface of the airfoil - dx = np.zeros(nE) # distance along two nodes - dv = np.zeros(nE) # speed difference according to element - alpha = np.zeros(nE) # angle of the element for the rotation matrix - for i in range(0,nE): - alpha[i] = np.arctan2((data[i+1,1]-data[i,1]),(data[i+1,0]-data[i,0])) - dx[i] = np.sqrt((data[i+1,0]-data[i,0])**2+(data[i+1,1]-data[i,1])**2) - xx[i+1] = dx[i]+xx[i] - vt[i] = (data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # tangential speed at node 1 element i - vt[i+1] = (data[i+1,3]*np.cos(alpha[i]) + data[i+1,4]*np.sin(alpha[i])) # tangential speed at node 2 element i - dv[i] = (vt[i+1] - vt[i])/dx[i] #velocity gradient according to element i. central scheme with half point - return xx, dv, vt, nN, alpha - - def __amplRoutine(self, Hk, Ret, theta): - '''Compute the amplitude of the TS waves envelop for the transition - ''' - Dgr = 0.08 - Hk1 = 3.5 - Hk2 = 4 - Hmi = (1/(Hk-1)) - logRet_crit = 2.492*(1/(Hk-1))**0.43 + 0.7*(np.tanh(14*Hmi-9.24)+1.0) # value at which the amplification starts to grow - if Ret <=0: - Ret = 1 - if np.log10(Ret) < logRet_crit - Dgr : - Ax = 0 - else: - Rnorm = (np.log10(Ret) - (logRet_crit-Dgr))/(2*Dgr) - if Rnorm <=0: - Rfac = 0 - if Rnorm >= 1: - Rfac = 1 - else: - Rfac = 3*Rnorm**2 - 2*Rnorm**3 - # normal envelope amplification rate - m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi) - arg = 3.87*Hmi-2.52 - dndRet = 0.028*(Hk-1)-0.0345*np.exp(-(arg**2)) - Ax = (m*dndRet/theta)*Rfac - # set correction for separated profiles - if Hk > Hk1: - Hnorm = (Hk - Hk1)/(Hk2-Hk1) - if Hnorm >=1: - Hfac = 1 - else: - Hfac = 3*Hnorm**2 - 2*Hnorm**3 - Ax1 = Ax - Gro = 0.3+0.35*np.exp(-0.15*(Hk-5)) - Tnr = np.tanh(1.2*(np.log10(Ret)-Gro)) - Ax2 = (0.086*Tnr - 0.25/(Hk-1)**1.5)/theta - if Ax2 < 0: - Ax2 = 0 - Ax = Hfac*Ax2 + (1-Hfac)*Ax1 - if Ax < 0: - Ax = 0 - return Ax - - def __laminarClosure(self, theta, H, Ret,Me, name): - ''' Laminar closure derived from the Falkner-Skan one-parameter profile family - Nishida and Drela (1996) - ''' - Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor - if name == "airfoil": - Hk = max(Hk, 1.02) - Hk = min(Hk,6) - else: - Hk = max(Hk, 1.00005) - Hk = min(Hk,6) - Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter - delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta) - Hks = Hk - 4.35 - if Hk <= 4.35: - dens = Hk + 1 - Hstar = 0.0111*Hks**2/dens - 0.0278*Hks**3/dens + 1.528 -0.0002*(Hks*Hk)**2 - elif Hk > 4.35: - Hstar = 1.528 + 0.015*Hks**2/Hk - Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction - if Hk < 5.5: - tmp = (5.5-Hk)**3 / (Hk+1) - Cf = (-0.07 + 0.0727*tmp)/Ret - elif Hk >= 5.5: - tmp = 1 - 1/(Hk-4.5) - Cf = (-0.07 + 0.015*tmp**2) /Ret - if Hk < 4: - Cd = ( 0.00205*(4.0-Hk)**5.5 + 0.207 ) * (Hstar/(2*Ret)) - elif Hk >= 4: - HkCd = (Hk-4)**2 - denCd = 1+0.02*HkCd - Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret)) - if name == 'wake': - Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)* (Hstar/(2*Ret)) - Cf = 0 - return Cd, Cf, Hstar, Hstar2, Hk, delta - - def __turbulentClosure(self, theta, H, Ct, Ret, Me, name): - ''' Turbulent closure derived from Nishida and Drela (1996) - ''' - # eliminate absurd transients - Ct = min(Ct, 0.30) - Ct = max(Ct, 0.0000001) - Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) - if name == 'wake': - Hk = max(Hk, 1.00005) - Hk = min(Hk,10) - else: - Hk = max(Hk, 1.00005) - Hk = min(Hk,10) - Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2 - gam = self.gamma - 1 - Fc = np.sqrt(1+0.5*gam*Me**2) - if Ret <= 400: - H0 = 4 - elif Ret > 400: - H0 = 3 + (400/Ret) - if Ret > 200: - Ret = Ret - elif Ret <= 200: - Ret = 200 - if Hk <= H0: - Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret - elif Hk > H0: - Hstar = (Hk-H0)**2*(0.007*np.log(Ret)/(Hk-H0+4/np.log(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret - Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction - logRt = np.log(Ret/Fc) - logRt = np.max((logRt,3)) - arg = np.max((-1.33*Hk, -20)) - Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity - delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta) - Ctcon = 0.5/(6.7**2*0.75) - if name == 'wake': - if Us > 0.99995: - Us = 0.99995 - n = 2 # two wake halves - Cf = 0 # no friction inside the wake - Hkc = Hk - 1 - Cdw = 0 # wall contribution - Cdd = (0.995-Us)*Ct**2 # turbulent outer layer contribution - Cdl = 0.15*(0.995-Us)**2/Ret # laminar stress contribution to outer layer - else: - if Us > 0.95: - Us = 0.98 - n = 1 - Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1)) - Hkc = max(Hk-1-18/Ret, 0.01) - # dissipation coefficient - Hmin = 1 + 2.1/np.log(Ret) - Fl = (Hk-1)/(Hmin-1) - Dfac = 0.5+0.5*np.tanh(Fl) - Cdw = 0.5*(Cf*Us) * Dfac - Cdd = (0.995-Us)*Ct**2 - Cdl = 0.15*(0.995-Us)**2/Ret - Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5 - Cd = n*(Cdw+ Cdd + Cdl) - Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2) - return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta - - - def __boundaryLayer(self, data, group, savef=False): - '''Discretize and solve the boundary layer equations for laminar/turbulent flows - ''' - xx, dv, vtOld, nN, alpha = self.__getBLcoordinates(data) - Me = data[:,5] - rhoe = data[:,6] - deltaStarOld = data[:,7] - if self.it == 0: - xxOld = xx - else: - xxOld = data[:,8] - #Initialize variables - blwVel = np.zeros(nN-1) # Blowing velocity - vt = np.zeros(nN) # Boundary layer velocity - theta = np.zeros(nN) # Momentum thickness - H = np.zeros(nN) # Shape factor - deltaStar = np.zeros(nN) # Displacement thickness - Hstar = np.zeros(nN) # Kinetic energy shape parameter - Hstar2 = np.zeros(nN) # Density shape parameter - Hk = np.zeros(nN) # Kinematic shape factor - Ret = np.zeros(nN) #Re based on momentum thickness - Cd = np.zeros(nN) # Dissipation coefficient - Cf = np.zeros(nN) # Skin friction - n = np.zeros(nN) # Logarithm of the maximum amplification ratio - Ax = np.zeros(nN) # Amplification factor - Ct = np.zeros(nN) # Shear stress coefficient - Cteq = np.zeros(nN) # Equilibrium shear stress coefficient - delta = np.zeros(nN) # BL thickness - # Boundary conditions - if group.name == 'airfoil': - theta[0], H[0], n[0], Ct[0] = group.initialConditions(self.Re, dv[0]) - elif group.name == 'wake': - theta[0] = group.T0 - H[0] = group.H0 - n[0] = group.n0 - Ct[0] = group.Ct0 - vt[0] = vtOld[0] - Ret[0] = vt[0]*theta[0]*self.Re - if Ret[0] < 1: - Ret[0] = 1 - vt[0] = Ret[0]/(theta[0]*self.Re) - deltaStar[0] = theta[0]*H[0] - # Define transition location ( N = 9) and compute stagnation point - ntr = 9 - Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name) - if n[0] < ntr: - turb = 0 # we begin with laminar flow - else : - turb = 1 # typically for the wake - xtr = 0 # if only turbulent flow - itTurb = 0 - # Solve the boundary layer equations - for i in range(1,nN): - # x[0] = theta; x[1] = H; x[2] = n for laminar/Ct for turbulent; x[3] = vt from interaction law - def f_lam(x): - f = np.zeros(len(x)) - Ret[i] = np.maximum(x[3]*x[0]*self.Re, 1) - Cd[i], Cf[i], Hstar[i], Hstar2[i], Hk[i], delta[i] = self.__laminarClosure(x[0], x[1], Ret[i],Me[i], group.name) - Ta = upw*x[0]+(1-upw)*theta[i-1] - Ha = upw*x[1]+(1-upw)*H[i-1] - uea = upw*x[3]+(1-upw)*vt[i-1] - Reta = upw*Ret[i] + (1-upw)*Ret[i-1] - Mea = upw*Me[i]+(1-upw)*Me[i-1] - Cda, Cfa, Hstara, Hstar2a, Hka, deltaa = self.__laminarClosure(Ta , Ha, Reta, Mea, group.name) - dx = xx[i] - xx[i-1] - Axa = self.__amplRoutine(Hka, Reta, Ta) - dTheta = x[0]-theta[i-1] - due = x[3]-vt[i-1] - dH = x[1]-H[i-1] - dn = x[2] - n[i-1] - dHstar = Hstar[i]-Hstar[i-1] - Hstara = upw*Hstar[i]+(1-upw)*Hstar[i-1] - f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2 - f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2 - f[2] = dn - dx*Axa - f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1]) - return f - def f_turb(x): - f = np.zeros(len(x)) - if group.name == 'wake': - dissipRatio = 0.9 # wall/wake dissipation length ratio - else: - dissipRatio = 1 - Ret[i] = x[3]*x[0]*self.Re - Ta = upw*x[0]+(1-upw)*theta[i-1] - xxa = upw*xx[i]+(1-upw)*xx[i-1] - Ha = upw*x[1]+(1-upw)*H[i-1] - Cta = upw*x[2]+(1-upw)*Ct[i-1] - uea = upw*x[3]+(1-upw)*vt[i-1] - Reta = np.maximum(upw*(x[3]*x[0]*self.Re)+(1-upw)*(vt[i-1]*theta[i-1]*self.Re), 1) - Mea = upw*Me[i]+(1-upw)*Me[i-1] - Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name) - Cda, Cfa, Hstara, Hstar2a, Hka,Cteqa, deltaa = self.__turbulentClosure(Ta, Ha,Cta, Reta,Mea, group.name) - dx = xx[i]-xx[i-1] - dTheta = x[0]-theta[i-1] - due = x[3]-vt[i-1] - dH = x[1]-H[i-1] - dHstar = Hstar[i]-Hstar[i-1] - dCt = x[2]-Ct[i-1] # here it is Ct^1/2 which is represented - Ta = upw*x[0]+(1-upw)*theta[i-1] - Cta = upw*x[2]+(1-upw)*Ct[i-1] - f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2 - f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2 - f[2] = (2*deltaa/Cta)*(dCt/dx) - 5.6*(Cteqa-Cta*dissipRatio)-2*deltaa*(4/(3*Ha*Ta)*(Cfa/2-((Hka-1)/(6.7*Hka*dissipRatio))**2)-1/uea * due/dx) - f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1]) - return f - # Laminar solver - if turb == 0: - if n[i-1] > 8 or i < 5: - upw = 1 # backward Euler - else: - upw = 0.5 # Trapezoidal scheme - x = np.array([theta[i-1],H[i-1], n[i-1], vt[i-1]]) - try: - sol = Newton.newton(f_lam, x, 200) - except RuntimeError as e: - sol = e.args[1] - print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Laminar solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET) - # Determine if the amplification waves start to grow or not - logRet_crit = 2.492*(1/(Hk[i]-1))**0.43 + 0.7*(np.tanh(14*(1/(Hk[i]-1))-9.24)+1.0) - if np.log10(Ret[i]) <= logRet_crit - 0.08 : - n[i] = 0 - else: - n[i] = sol[2] - xtr = 1 - # Turbulent solver - elif turb == 1: - if i <2 and group.name =='wake': - upw = 1 # begin of the wake - elif itTurb <2 and group.name == 'airfoil': - upw = 1 - itTurb += 1 - else: - upw = 0.5 # trapezoidal scheme - x = np.array([theta[i-1], H[i-1], Ct[i-1], vt[i-1]]) - try: - sol = Newton.newton(f_turb, x, 200) - except RuntimeError as e: - sol = e.args[1] - print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Turbulent solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET) - Ct[i] = sol[2] - theta[i] = sol[0] - H[i] = sol[1] - vt[i] = sol[3] - # Transition solver - if n[i] >= ntr and turb ==0: - turb = 1 - upw = 1 - xtr = (ntr-n[i-1])*((data[i,0]-data[i-1,0])/(n[i]-n[i-1]))+data[i-1,0] - avTurb = (data[i,0]-xtr)/(data[i,0]-data[i-1,0]) # percentage of the subsection corresponding to a turbulent flow - avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow - # Averaging both solutions - Cteq[i-1]= self.__turbulentClosure(theta[i-1], H[i-1], 0.03, Ret[i-1], Me[i-1], group.name)[5] - Ct[i-1] = avLam*Cteq[i-1] * 0.7 - x_turb = np.array([theta[i-1], 1.515, Ct[i-1] , vt[i-1]]) - try: - sol_turb = Newton.newton(f_turb, x_turb, 200) - except RuntimeError as e: - sol = e.args[1] - print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Transition solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET) - theta[i] = avLam*sol[0]+avTurb*sol_turb[0] - H[i] = avLam*sol[1]+avTurb*sol_turb[1] - if group.name == 'airfoil': - Ct[i] = max(avTurb*sol_turb[2],0.03) - else: - Ct[i] = max(avTurb*sol_turb[2],0.05) - vt[i] = avLam*sol[3]+avTurb*sol_turb[3] - Cd[i-1], Cf[i-1], Hstar[i-1], Hstar2[i-1], Hk[i-1], delta[i-1] = self.__laminarClosure(theta[i-1], H[i-1], Ret[i-1],Me[i-1], group.name) - Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name) - deltaStar[i] = H[i]*theta[i] - blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1])) - # Normalize the friction and dissipation coefficients by the freestream velocity - Cf = vt**2*Cf - Cd = vt**2*Cd - # Compute the various drag coefficients (total, friction, profile) - Cdtot = 2*theta[-1]*(vt[-1]**((H[-1]+5)/2)) - Cdf = np.trapz(Cf, data[:,0]) - Cdp = Cdtot-Cdf - # Outputs - blScalars = [Cdtot, Cdf, Cdp] - blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta] - return blwVel, xtr, xx, blScalars, blVectors - - - def run(self, group): - ''' Run the viscous solver to solve the BL equations and give the outputs to the coupler - ''' - group.connectListNodes, group.connectListElems, data = group.connectList() - # Compute BLE for the airfoil (upper section(data[0]) and lower section (data[1])) - if len(data) == 2: - ublwVel, xtrTop, uxx, ublScalars, ublVectors = self.__boundaryLayer(data[0], group, True) - lblwVel, xtrBot, lxx, lblScalars, lblVectors = self.__boundaryLayer(data[1], group) - # Put the upper and lower values together - self.data = np.concatenate((data[0], data[1])) - self.blScalars = np.add(ublScalars, lblScalars) - self.blVectors = np.hstack((ublVectors, lblVectors)) - blwVel = np.concatenate((ublwVel, lblwVel)) - self.xtr = [xtrTop, xtrBot] - # Boundary conditions for the wake - group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]] - group.HEnd = [ublVectors[2][-1], lblVectors[2][-1]] - group.CtEnd = [ublVectors[8][-1], lblVectors[8][-1]] - # Delete the stagnation point doublon for variables in the VII loop - lblVectors[0] = np.delete(lblVectors[0],0,0) #deltastar - lxx = np.delete(lxx,0,0) # airfoil coordinates - group.deltaStar = np.concatenate((ublVectors[0], lblVectors[0])) - group.xx = np.concatenate((uxx, lxx)) - # Compute BLE for the wake - else: - blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group) - group.deltaStar = blVectors[0] - # The drag is measured at the end of the wake - self.Cd = blScalars[0] - self.Cdp = self.Cd-self.blScalars[1] # Cdf comes from the airfoil as there is no friction in the wake - self.blScalars[0] = self.Cd - self.blScalars[2] = self.Cdp - self.Cdf = self.blScalars[1] - # Sort the following results in reference frame - group.deltaStar = group.deltaStar[group.connectListNodes.argsort()] - group.xx = group.xx[group.connectListNodes.argsort()] - group.u = blwVel[group.connectListElems.argsort()] diff --git a/vii/CMakeLists.txt b/vii/CMakeLists.txt new file mode 100644 index 0000000..5bec1c8 --- /dev/null +++ b/vii/CMakeLists.txt @@ -0,0 +1,27 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Add source dir +ADD_SUBDIRECTORY( src ) +ADD_SUBDIRECTORY( _src ) + +# Add test dir +MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests) + +# Add to install +INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py + ${CMAKE_CURRENT_LIST_DIR}/utils.py + DESTINATION vii) +INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/pyVII + DESTINATION vii) diff --git a/vii/__init__.py b/vii/__init__.py new file mode 100644 index 0000000..43be8da --- /dev/null +++ b/vii/__init__.py @@ -0,0 +1,21 @@ +# -*- coding: utf-8 -*- +# dart MODULE initialization file + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +import fwk +import tbox +from viiw import * diff --git a/vii/_src/CMakeLists.txt b/vii/_src/CMakeLists.txt new file mode 100644 index 0000000..c49de25 --- /dev/null +++ b/vii/_src/CMakeLists.txt @@ -0,0 +1,35 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# CMake input file of the SWIG wrapper around "viiw.so" + +INCLUDE(${SWIG_USE_FILE}) + +FILE(GLOB SRCS *.h *.cpp *.inl *.swg) +FILE(GLOB ISRCS *.i) + +SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON) +SET(CMAKE_SWIG_FLAGS "-interface" "_viiw") # avoids "import _module_d" with MSVC/Debug +SWIG_ADD_LIBRARY(viiw LANGUAGE python SOURCES ${ISRCS} ${SRCS}) +SET_PROPERTY(TARGET viiw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON) +MACRO_DebugPostfix(viiw) + +TARGET_INCLUDE_DIRECTORIES(viiw PRIVATE ${PROJECT_SOURCE_DIR}/vii/_src + ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src + ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src + ${PYTHON_INCLUDE_PATH}) +TARGET_LINK_LIBRARIES(viiw PRIVATE vii tbox fwk ${PYTHON_LIBRARIES}) + +INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/viiw.py DESTINATION ${CMAKE_INSTALL_PREFIX}) +INSTALL(TARGETS viiw DESTINATION ${CMAKE_INSTALL_PREFIX}) diff --git a/vii/_src/viiw.h b/vii/_src/viiw.h new file mode 100644 index 0000000..ec45463 --- /dev/null +++ b/vii/_src/viiw.h @@ -0,0 +1,17 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "DDriver.h" \ No newline at end of file diff --git a/vii/_src/viiw.i b/vii/_src/viiw.i new file mode 100644 index 0000000..cc8f9da --- /dev/null +++ b/vii/_src/viiw.i @@ -0,0 +1,57 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +// SWIG input file of the 'dart' module + +%feature("autodoc","1"); + +%module(docstring= +"'viiw' module: Viscous inviscid interaction for dart 2021-2022 +(c) ULiege - A&M", +directors="0", +threads="1" +) viiw +%{ + +#include "vii.h" + +#include "fwkw.h" +#include "tboxw.h" +#include "viiw.h" + +%} + + +%include "fwkw.swg" + +// ----------- MODULES UTILISES ------------ +%import "tboxw.i" + +// ----------- VII CLASSES ---------------- +%include "vii.h" + +%shared_ptr(vii::Driver); + +%immutable vii::Driver::Re; // read-only variable (no setter) +%immutable vii::Driver::Cdt; +%immutable vii::Driver::Cdt_sec; +%immutable vii::Driver::Cdf; +%immutable vii::Driver::Cdf_sec; +%immutable vii::Driver::Cdp; +%immutable vii::Driver::Cdp_sec; +%immutable vii::Driver::CFL0; +%immutable vii::Driver::verbose; +%include "DDriver.h" \ No newline at end of file diff --git a/dart/viscous/__init__.py b/vii/api/__init__.py similarity index 100% rename from dart/viscous/__init__.py rename to vii/api/__init__.py diff --git a/vii/api/core.py b/vii/api/core.py new file mode 100644 index 0000000..19ccbb8 --- /dev/null +++ b/vii/api/core.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2022 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +## Initialize VII computation +# Paul Dechamps + +def initVII(cfg, icfg, iSolverName='DART'): + """ + Inputs + ------ + - cfg: Dict with options + - Re (float): Flow Reynolds number. + - Minf (float): Freestream Mach number (only used for time step computation). + - CFL0 (float): Initial CFL number. + - Verb (int): Verbosity level of the viscous solver. + + - couplIter (int): Maximum number of coupling iterations. + - couplTol (float): Tolerance to end computation (drag count between 2 iterations). + - iterPrint (int): Number of iterations between which the coupler outputs its state. + - resetInv (bool): Resets the inviscid solver at each iteration (True), reuses previous + state as initial condition (False). + Note + All inputs have default values except Re. + - icg: Dict with options. Inviscid solver configuration. + - iSolverName (string): Name of the inviscid solver. + + Outputs + ------- + dict with keys + - coupler: Coupler python object to run the viscous-inviscid algorithm. + - solversAPI: Interface between the solver objects to ensure proper communication. + - vSolver: vii::Driver object to run the viscous computation. + - iSolverObjects: Dict containing. + - All inviscid solver dependencies (depend on iSolverName). + + @todo: - Correctly handle 3D computation parameters. + """ + # Imports + from vii.coupler import Coupler + import vii + + if iSolverName == 'DART': + from vii.interfaces.dart.DartInterface import DartInterface as interface + from dart.api.core import initDart + iSolverObjects = initDart(icfg, scenario=icfg['scenario'], task=icfg['task'], viscous = 1) + + # Check viscous solver parameters. + if 'Re' in cfg and cfg['Re'] > 0: + _Re = cfg['Re'] + else: + raise RuntimeError('Missing or invalid Reynolds number') + if 'Minf' in cfg and cfg['Minf'] > 0: + _Minf = cfg['Minf'] + else: + _Minf = 0.1 + if 'CFL0' in cfg and cfg['CFL0'] > 0: + _CFL0 = cfg['CFL0'] + else: + _CFL0 = 1 + if 'Verb' in cfg and 0<= cfg['Verb'] <= 3: + _verbose = cfg['Verb'] + else: + _verbose = 1 + _span = 0 + _nSections = 1 + + # Check coupler parameters. + if 'couplIter' in cfg and cfg['couplIter'] > 0: + __couplIter = cfg['couplIter'] + else: + __couplIter = 150 + if 'couplTol' in cfg: + __couplTol = cfg['couplTol'] + else: + __couplTol = 1e-4 + if 'iterPrint' in cfg: + __iterPrint = cfg['iterPrint'] + else: + __iterPrint = 1 + if 'resetInv' in cfg: + __resetInv = cfg['resetInv'] + else: + __resetInv = False + + # Viscous solver object. + vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _span, verbose =_verbose) + # Solvers interface. + solversAPI = interface(iSolverObjects['sol'], vSolver, [iSolverObjects['blwb'], iSolverObjects['blww']]) + # Coupler + coupler = Coupler(solversAPI, vSolver, _maxCouplIter = __couplIter, _couplTol = __couplTol, _iterPrint = __iterPrint, _resetInv = __resetInv) + + return {'coupler': coupler, + 'solversAPI': solversAPI, + 'vSolver': vSolver, + 'iSolverObjects': iSolverObjects} diff --git a/vii/coupler.py b/vii/coupler.py new file mode 100644 index 0000000..7c6ead6 --- /dev/null +++ b/vii/coupler.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# VII Coupler +# Paul Dechamps + +import fwk +from fwk.coloring import ccolors +import math +import numpy as np + +class Coupler: + def __init__(self, iSolverAPI, vSolver, _maxCouplIter=150, _couplTol=1e-4, _iterPrint=1, _resetInv=False): + self.iSolverAPI = iSolverAPI + self.vSolver = vSolver + self.maxCouplIter = _maxCouplIter + self.couplTol = _couplTol + self.resetInv = _resetInv + self.iterPrint = _iterPrint if self.iSolverAPI.getVerbose() == 0 and self.vSolver.verbose == 0 else 1 + self.tms = fwk.Timers() + + def run(self): + # Aerodynamic coefficients. + aeroCoeffs = np.zeros((0, 2)) + + # Convergence parameters. + couplIter = 0 + cdPrev = 0.0 + + while couplIter < self.maxCouplIter: + # Run inviscid solver. + self.tms['inviscid'].start() + if self.resetInv: + self.iSolverAPI.reset() + self.iSolverAPI.run() + self.tms['inviscid'].stop() + + # Write inviscid Cp distribution. + if couplIter == 0: + self.iSolverAPI.writeCp() + + # Impose inviscid boundary in the viscous solver. + self.tms['processing'].start() + self.iSolverAPI.imposeInvBC(couplIter) + self.tms['processing'].stop() + + # Run viscous solver. + self.tms['viscous'].start() + exitCode = self.vSolver.run(couplIter) + self.tms['viscous'].stop() + + # Impose blowing boundary condition in the inviscid solver. + self.tms['processing'].start() + self.iSolverAPI.imposeBlwVel() + self.tms['processing'].stop() + + aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolverAPI.getCl(), self.vSolver.Cdt])) + + # Check convergence status. + error = abs((self.vSolver.Cdt - cdPrev) / self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1 + if error <= self.couplTol: + print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET) + return aeroCoeffs + cdPrev = self.vSolver.Cdt + + if couplIter == 0: + print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'.format(self.iSolverAPI.getAoA()*180/math.pi, self.iSolverAPI.getMinf(), self.vSolver.getRe()/1e6)) + print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error')) + print(' ----------------------------------------------------------') + if couplIter % self.iterPrint == 0: + print(' {:>4.0f}| {:>10.4f} {:>10.4f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error))) + if self.iSolverAPI.getVerbose() != 0 or self.vSolver.verbose != 0: + print('') + couplIter += 1 + self.tms['processing'].stop() + else: + print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET) + print('') + print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error')) + print(' ----------------------------------------------------------') + print(ccolors.ANSI_RED, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET) + return aeroCoeffs \ No newline at end of file diff --git a/vii/interfaces/dart/DartInterface.py b/vii/interfaces/dart/DartInterface.py new file mode 100644 index 0000000..217cc44 --- /dev/null +++ b/vii/interfaces/dart/DartInterface.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +# Dart interface +# Paul Dechamps + +import numpy as np + +import dart.utils as floU +import tbox.utils as tboxU +from ..viscous.DAirfoil import Airfoil +from ..viscous.DWake import Wake +from ..viscous import DGroupBuilder + +class DartInterface: + def __init__(self, dartSolver, vSolver, blw): + self.solver = dartSolver + self.vSolver = vSolver + self.createProblem(blw[0], blw[1]) + self.nElmUpperPrev = 0 + + def createProblem(self, _boundaryAirfoil, _boundaryWake): + """Create data structure for information transfer. + + Args + ---- + - _boundaryAirfoil: dart::Blowing data structure for the airfoil. + - _boundaryWake: dart::Blowing data structure for the wake. + """ + self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] + self.blwVel = [None, None, None] + self.deltaStar = [None, None, None] + self.xx = [None, None, None] + + def run(self): + return self.solver.run() + + def writeCp(self): + """ Write Cp distribution around the airfoil on file. + """ + Cp = floU.extract(self.group[0].boundary.tag.elems, self.solver.Cp) + tboxU.write(Cp, 'Cp_Inviscid.dat', '%1.5e',', ', 'x, y, z, Cp', '') + + def save(self, writer): + self.solver.save(writer) + + def getAoA(self): + return self.solver.pbl.alpha + + def getMinf(self): + return self.solver.pbl.M_inf + + def getCl(self): + return self.solver.Cl + + def getCm(self): + return self.solver.Cm + + def getVerbose(self): + return self.solver.verbose + + def reset(self): + self.solver.reset() + + def imposeInvBC(self, couplIter): + """ Extract the inviscid paramaters at the edge of the boundary layer + and use it as a boundary condition in the viscous solver. + + Params + ------ + - couplIter (int): Coupling iteration. + """ + # Retreive parameters. + for n in range(0, len(self.group)): + for i in range(0, len(self.group[n].boundary.nodes)): + self.group[n].v[i,0] = self.solver.U[self.group[n].boundary.nodes[i].row][0] + self.group[n].v[i,1] = self.solver.U[self.group[n].boundary.nodes[i].row][1] + self.group[n].v[i,2] = self.solver.U[self.group[n].boundary.nodes[i].row][2] + self.group[n].Me[i] = self.solver.M[self.group[n].boundary.nodes[i].row] + self.group[n].rhoe[i] = self.solver.rho[self.group[n].boundary.nodes[i].row] + + dataIn = [None, None] + for n in range(0, len(self.group)): + self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList(couplIter) + self.fMeshDict, _, _ = DGroupBuilder.mergeVectors(dataIn, 0, 1) + fMeshDict = self.fMeshDict.copy() + + if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter == 0): + # Initialize mesh + self.vSolver.setGlobMesh(0, 0, fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['yGlobCoord'][:fMeshDict['bNodes'][1]], fMeshDict['zGlobCoord'][:fMeshDict['bNodes'][1]]) + self.vSolver.setGlobMesh(0, 1, fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setGlobMesh(0, 2, fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][2]:], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][2]:]) + + # Initialize parameters at the edge of the boundary layer. + self.vSolver.setxxExt(0, 0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) + self.vSolver.setxxExt(0, 1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setxxExt(0, 2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:]) + self.vSolver.setDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]]) + self.vSolver.setDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:]) + else: + # Parameters at the edge of the boundary layer only. + self.vSolver.setxxExt(0, 0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]]) + self.vSolver.setxxExt(0, 1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setxxExt(0, 2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:]) + self.vSolver.setDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]]) + self.vSolver.setDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:]) + + # Inviscid boundary condition. + self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) + self.vSolver.setVelocities(0, 0, fMeshDict['vx'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]]) + self.vSolver.setVelocities(0, 1, fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setVelocities(0, 2, fMeshDict['vx'][fMeshDict['bNodes'][2]:], fMeshDict['vy'][fMeshDict['bNodes'][2]:], fMeshDict['vz'][fMeshDict['bNodes'][2]:]) + self.vSolver.setMe(0, 0, fMeshDict['Me'][:fMeshDict['bNodes'][1]]) + self.vSolver.setMe(0, 1, fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setMe(0, 2, fMeshDict['Me'][fMeshDict['bNodes'][2]:]) + self.vSolver.setRhoe(0, 0, fMeshDict['rho'][:fMeshDict['bNodes'][1]]) + self.vSolver.setRhoe(0, 1, fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]]) + self.vSolver.setRhoe(0, 2, fMeshDict['rho'][fMeshDict['bNodes'][2]:]) + + def imposeBlwVel(self): + """ Extract the solution of the viscous calculation (blowing velocity) + and use it as a boundary condition in the inviscid solver. + """ + # Retreive and set blowing velocities. + for iRegion in range(3): + self.blwVel[iRegion] = np.asarray(self.vSolver.extractBlowingVel(0, iRegion)) + self.deltaStar[iRegion] = np.asarray(self.vSolver.extractDeltaStar(0, iRegion)) + self.xx[iRegion] = np.asarray(self.vSolver.extractxx(0, iRegion)) + blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1])) + + # Remove stagnation point doublon. + deltaStarAF = np.concatenate((self.deltaStar[0], self.deltaStar[1][1:])) + xxAF = np.concatenate((self.xx[0], self.xx[1][1:])) + + # Blowing velocity, displacement thickness and mesh (local coordinates) distributions in the wake. + blwVelWK = self.blwVel[2] + deltaStarWK = self.deltaStar[2] + xxWK = self.xx[2] + # Update data structures for information transfer. + self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()] + self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()] + self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()] + self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()] + self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()] + self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()] + + # Impose blowing velocities. + for n in range(0, len(self.group)): + for i in range(0, self.group[n].nE): + self.group[n].boundary.setU(i, self.group[n].u[i]) \ No newline at end of file diff --git a/vii/interfaces/dart/__init__.py b/vii/interfaces/dart/__init__.py new file mode 100644 index 0000000..7c68785 --- /dev/null +++ b/vii/interfaces/dart/__init__.py @@ -0,0 +1 @@ +# -*- coding: utf-8 -*- \ No newline at end of file diff --git a/dart/viscous/airfoil.py b/vii/interfaces/viscous/DAirfoil.py old mode 100644 new mode 100755 similarity index 93% rename from dart/viscous/airfoil.py rename to vii/interfaces/viscous/DAirfoil.py index 2bd069d..79d9494 --- a/dart/viscous/airfoil.py +++ b/vii/interfaces/viscous/DAirfoil.py @@ -19,7 +19,8 @@ ## Airfoil around which the boundary layer is computed # Amaury Bilocq -from dart.viscous.boundary import Boundary +from .DBoundary import Boundary + import numpy as np class Airfoil(Boundary): @@ -29,12 +30,12 @@ class Airfoil(Boundary): self.T0 = 0 # initial condition for the momentum thickness self.H0 = 0 self.n0 = 0 - self.Ct0 = 0 + self.ct0 = 0 self.TEnd = [0,0] self.HEnd = [0,0] - self.CtEnd = [0,0] + self.ctEnd = [0,0] self.nEnd = [0,0] - + def initialConditions(self, Re, dv): if dv > 0: self.T0 = np.sqrt(0.075/(Re*dv)) @@ -42,12 +43,13 @@ class Airfoil(Boundary): self.T0 = 1e-6 self.H0 = 2.23 # One parameter family Falkner Skan self.n0 = 0 - self.Ct0 = 0 - return self.T0, self.H0, self.n0, self.Ct0 + self.ct0 = 0 + return self.T0, self.H0, self.n0, self.ct0 + - def connectList(self): - ''' Sort the value read by the viscous solver/ Create list of connectivity - ''' + def connectList(self, couplIter): + """ Sort the value read by the viscous solver/ Create list of connectivity + """ N1 = np.zeros(self.nN, dtype=int) # Node number connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems @@ -72,7 +74,9 @@ class Airfoil(Boundary): eData[i,1] = self.boundary.tag.elems[i].nodes[0].no eData[i,2] = self.boundary.tag.elems[i].nodes[1].no # Find the stagnation point - idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2)) + if couplIter == 0: + self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2)) + idxStag = self.idxStag globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes # Find TE nodes (assuming suction side element is above pressure side element in the y direction) idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes diff --git a/dart/viscous/boundary.py b/vii/interfaces/viscous/DBoundary.py old mode 100644 new mode 100755 similarity index 100% rename from dart/viscous/boundary.py rename to vii/interfaces/viscous/DBoundary.py diff --git a/vii/interfaces/viscous/DGroupBuilder.py b/vii/interfaces/viscous/DGroupBuilder.py new file mode 100755 index 0000000..d4fcbca --- /dev/null +++ b/vii/interfaces/viscous/DGroupBuilder.py @@ -0,0 +1,169 @@ +import numpy as np + +def mergeVectors(dataIn, mapMesh, nFactor): + """Merges 3 groups upper/lower sides and wake. + """ + for k in range(len(dataIn)): + if len(dataIn[k]) == 2: # Airfoil case. + xx_up, dv_up, vtExt_up, _, alpha_up = getBLcoordinates(dataIn[k][0]) + xx_lw, dv_lw, vtExt_lw, _, alpha_lw = getBLcoordinates(dataIn[k][1]) + else: # Wake case + xx_wk, dv_wk, vtExt_wk, _, alpha_wk = getBLcoordinates(dataIn[k]) + + if mapMesh: + # Save initial mesh. + cMesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][1:, 0], dataIn[1][:, 0])) + cMeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][1:, 0])] + cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk)) + cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk)) + cxxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][1:, 8], dataIn[1][:, 8])) + + # Create fine mesh. + fMeshUp = createFineMesh(dataIn[0][0][:, 0], nFactor) + fMeshLw = createFineMesh(dataIn[0][1][:, 0], nFactor) + fMeshWk = createFineMesh(dataIn[1][:, 0], nFactor) + fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk)) + fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])] + + fxx_up = interpolateToFineMesh(xx_up, fMeshUp, nFactor) + fxx_lw = interpolateToFineMesh(xx_lw, fMeshLw, nFactor) + fxx_wk = interpolateToFineMesh(xx_wk, fMeshWk, nFactor) + fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk)) + + fvtExt_up = interpolateToFineMesh(vtExt_up, fMeshUp, nFactor) + fvtExt_lw = interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor) + fvtExt_wk = interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor) + fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk)) + + fMe_up = interpolateToFineMesh(dataIn[0][0][:, 5], fMeshUp, nFactor) + fMe_lw = interpolateToFineMesh(dataIn[0][1][:, 5], fMeshLw, nFactor) + fMe_wk = interpolateToFineMesh(dataIn[1][:, 5], fMeshWk, nFactor) + fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk)) + + frho_up = interpolateToFineMesh(dataIn[0][0][:, 6], fMeshUp, nFactor) + frho_lw = interpolateToFineMesh(dataIn[0][1][:, 6], fMeshLw, nFactor) + frho_wk = interpolateToFineMesh(dataIn[1][:, 6], fMeshWk, nFactor) + frho = np.concatenate((frho_up, frho_lw[1:], frho_wk)) + + fdStarExt_up = interpolateToFineMesh(dataIn[0][0][:, 7], fMeshUp, nFactor) + fdStarExt_lw = interpolateToFineMesh(dataIn[0][1][:, 7], fMeshLw, nFactor) + fdStarExt_wk = interpolateToFineMesh(dataIn[1][:, 7], fMeshWk, nFactor) + + fxxExt_up = interpolateToFineMesh(dataIn[0][0][:, 8], fMeshUp, nFactor) + fxxExt_lw = interpolateToFineMesh(dataIn[0][1][:, 8], fMeshLw, nFactor) + fxxExt_wk = interpolateToFineMesh(dataIn[1][:, 8], fMeshWk, nFactor) + + fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk)) + fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk)) + fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk)) + + # Create dictionnaries for the fine and coarse meshes. + fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx, + 'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv} + + cMe = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][1:, 5], dataIn[1][:, 5])) + cRho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][1:, 6], dataIn[1][:, 6])) + cdStarExt = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][1:, 7], dataIn[1][:, 7])) + cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx, + 'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv} + + dataUpper = np.zeros((len(fMeshUp), 0)) + dataLower = np.zeros((len(fMeshLw), 0)) + dataWake = np.zeros((len(fMeshWk), 0)) + for iParam in range(len(dataIn[0][0][0, :])): + interpParamUp = interpolateToFineMesh(dataIn[0][0][:, iParam], fMeshUp, nFactor) + dataUpper = np.column_stack((dataUpper, interpParamUp)) + interpParamLw = interpolateToFineMesh(dataIn[0][1][:, iParam], fMeshLw, nFactor) + dataLower = np.column_stack((dataLower, interpParamLw)) + interpParamWk = interpolateToFineMesh(dataIn[1][:, iParam], fMeshWk, nFactor) + dataWake = np.column_stack((dataWake, interpParamWk)) + # Remove stagnation point doublon. + dataLower = np.delete(dataLower, 0, axis=0) + else: + Mesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][:, 0], dataIn[1][:, 0])) + Meshy = np.concatenate((dataIn[0][0][:, 1], dataIn[0][1][:, 1], dataIn[1][:, 1])) + Meshz = np.concatenate((dataIn[0][0][:, 2], dataIn[0][1][:, 2], dataIn[1][:, 2])) + MeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]) + len(dataIn[1][:, 0])] + + # Concatenate all parameters (upper side, lower side, wake) + xx = np.concatenate((xx_up, xx_lw, xx_wk)) + vtExt = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk)) + vx = np.concatenate((dataIn[0][0][:, 3], dataIn[0][1][:, 3], dataIn[1][:, 3])) + vy = np.concatenate((dataIn[0][0][:, 4], dataIn[0][1][:, 4], dataIn[1][:, 4])) + vz = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5])) + alpha = np.concatenate((alpha_up, alpha_lw, alpha_wk)) + dv = np.concatenate((dv_up, dv_lw, dv_wk)) + Me = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5])) + rho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][:, 6], dataIn[1][:, 6])) + dStarext = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][:, 7], dataIn[1][:, 7])) + xxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][:, 8], dataIn[1][:, 8])) + + # Create dictionnaries for the fine and coarse meshes which are equivalent if meshMap is disabled. + fMeshDict = {'globCoord': Mesh, 'yGlobCoord': Meshy, 'zGlobCoord': Meshz, 'bNodes': MeshbNodes, 'locCoord': xx, + 'vtExt': vtExt, 'vx': vx, 'vy': vy, 'vz': vz, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv} + cMeshDict = fMeshDict.copy() + dataUpper = dataIn[0][0] + dataLower = dataIn[0][1][1:, :] + dataWake = dataIn[1] + + data = np.concatenate((dataUpper, dataLower, dataWake), axis=0) + + return fMeshDict, cMeshDict, data + +def getBLcoordinates(data): + """Transform the reference coordinates into airfoil coordinates. + """ + nN = len(data[:, 0]) + nE = nN-1 # Nbr of element along the surface. + vt = np.zeros(nN) # Outer flow velocity. + xx = np.zeros(nN) # Position along the surface of the airfoil. + dx = np.zeros(nE) # Distance along two nodes. + dv = np.zeros(nE) # Speed difference according to element. + alpha = np.zeros(nE) # Angle of the element for the rotation matrix. + + for i in range(0, nE): + alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0])) + dx[i] = np.sqrt((data[i+1, 0] - data[i, 0]) **2 + (data[i+1, 1]-data[i, 1])**2) + xx[i+1] = dx[i] + xx[i] + # Tangential speed at node 1 element i + vt[i] = (data[i, 3] * np.cos(alpha[i]) + data[i, 4] * np.sin(alpha[i])) + # Tangential speed at node 2 element i + vt[i+1] = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i])) + # Velocity gradient according to element i. + dv[i] = (vt[i+1] - vt[i]) / dx[i] + return xx, dv, vt, nN, alpha + +def interpolateToFineMesh(cVector, fMesh, nFactor): + """ Interpolate variables of cVector on the fine mesh fMesh. + + Params + ------ + - cVector (numpy.ndarray): Distribution to be interpolated. + - fMesh (numpy.ndarray): Local tangential coordinate (xi) of the fine mesh. + - nFactor (int): Number of times each cell was split when creating the fine mesh. + """ + fVector = np.zeros(len(fMesh)-1) + for cPoint in range(len(cVector) - 1): + dv = cVector[cPoint + 1] - cVector[cPoint] + for fPoint in range(nFactor): + fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + \ + fPoint * dv / nFactor + fVector = np.append(fVector, cVector[-1]) + return fVector + +def createFineMesh(cMesh, nFactor): + """ Create fine mesh distribution in the local frame of reference. + + Params + ------ + - cMesh (numpy.ndarray): Initial mesh (used for the inviscid calculations. + - nFactor (int): Number of times each cell will be split to create the fine mesh. + """ + fMesh = [] + for iPoint in range(len(cMesh) - 1): + dx = cMesh[iPoint + 1] - cMesh[iPoint] + for iRefinedPoint in range(nFactor): + fMesh = np.append( + fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor) + fMesh = np.append(fMesh, cMesh[-1]) + return fMesh diff --git a/dart/viscous/wake.py b/vii/interfaces/viscous/DWake.py old mode 100644 new mode 100755 similarity index 90% rename from dart/viscous/wake.py rename to vii/interfaces/viscous/DWake.py index 91e1f32..9e770a9 --- a/dart/viscous/wake.py +++ b/vii/interfaces/viscous/DWake.py @@ -17,10 +17,9 @@ ## Wake behind airfoil (around which the boundary layer is computed) -# Amaury Bilocq +# Amaury Bilocq & Paul Dechamps -from dart.viscous.boundary import Boundary -from dart.viscous.airfoil import Airfoil +from .DBoundary import Boundary import numpy as np class Wake(Boundary): @@ -30,11 +29,11 @@ class Wake(Boundary): self.T0 = 0 # initial condition for the momentum thickness self.H0 = 0 self.n0 = 9 # wake is always turbulent - self.Ct0 = 0 + self.ct0 = 0 - def connectList(self): - ''' Sort the value read by the viscous solver/ Create list of connectivity - ''' + def connectList(self, couplIter): + """ Sort the value read by the viscous solver/ Create list of connectivity + """ N1 = np.zeros(self.nN, dtype=int) # Node number connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems @@ -59,7 +58,8 @@ class Wake(Boundary): eData[i,1] = self.boundary.tag.elems[i].nodes[0].no eData[i,2] = self.boundary.tag.elems[i].nodes[1].no connectListNodes = data[:,1].argsort() - connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0] + # connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0] + connectListElems[0] = 0 for i in range(1, len(eData[:,0])): connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0] data[:,0] = data[connectListNodes,0] diff --git a/vii/src/CMakeLists.txt b/vii/src/CMakeLists.txt new file mode 100644 index 0000000..1e17fcd --- /dev/null +++ b/vii/src/CMakeLists.txt @@ -0,0 +1,27 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# CMake input file of dart.so + +FILE(GLOB SRCS *.h *.cpp *.inl *.hpp) + +ADD_LIBRARY(vii SHARED ${SRCS}) +MACRO_DebugPostfix(vii) +TARGET_INCLUDE_DIRECTORIES(vii PUBLIC ${PROJECT_SOURCE_DIR}/vii/src) + +TARGET_LINK_LIBRARIES(vii tbox) + +INSTALL(TARGETS vii DESTINATION ${CMAKE_INSTALL_PREFIX}) + +SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)") diff --git a/vii/src/DBoundaryLayer.cpp b/vii/src/DBoundaryLayer.cpp new file mode 100644 index 0000000..6301a13 --- /dev/null +++ b/vii/src/DBoundaryLayer.cpp @@ -0,0 +1,184 @@ +#include "DBoundaryLayer.h" +#include "DClosures.h" +#include "DDiscretization.h" +#include "DEdge.h" +#include <cmath> +#include <iomanip> + +using namespace vii; + +/** + * @brief Construct a new BoundaryLayer::BoundaryLayer object. + * + * @param _xtrF Forced transition location (default=-1; free transition). + */ +BoundaryLayer::BoundaryLayer(double _xtrF) +{ + if (_xtrF < 0. && _xtrF != -1.) + throw std::runtime_error("Fatal error: Invalid transition location.\n"); + + xtrF = _xtrF; + xtr = 1.0; + closSolver = new Closures(); + mesh = new Discretization(); + blEdge = new Edge(); + transMarker = 0; +} + +BoundaryLayer::~BoundaryLayer() +{ + delete closSolver; + delete mesh; + delete blEdge; + // std::cout << "~vii::BoundaryLayer()\n"; +} + +/** + * @brief Set the mesh and resize all boundary layer quantities accordingly. + * + * @param x Position x in the global frame of reference. + * @param y Position y in the global frame of reference. + * @param z Position z in the global frame of reference. + */ +void BoundaryLayer::setMesh(std::vector<double> x, + std::vector<double> y, + std::vector<double> z) +{ + size_t nMarkers = x.size(); + mesh->setGlob(x, y, z); + + // Resize and reset all variables if needed. + if (mesh->getDiffnElm() != 0) + { + u.resize(nVar * nMarkers, 0.); + Ret.resize(nMarkers, 0.); + cd.resize(nMarkers, 0.); + cf.resize(nMarkers, 0.); + Hstar.resize(nMarkers, 0.); + Hstar2.resize(nMarkers, 0.); + Hk.resize(nMarkers, 0.); + ctEq.resize(nMarkers, 0.); + us.resize(nMarkers, 0.); + deltaStar.resize(nMarkers, 0.); + delta.resize(nMarkers, 0.); + regime.resize(nMarkers, 0); + } +} + +/** + * @brief Set the boundary condition at the stagnation point. + * + * @param Re Freestream Reynolds number. + */ +void BoundaryLayer::setStagBC(double Re) +{ + double dv0 = (blEdge->getVt(1) - blEdge->getVt(0)) / (mesh->getLoc(1) - mesh->getLoc(0)); + if (dv0 > 0) + u[0] = sqrt(0.075 / (Re * dv0)); // Theta + else + u[0] = 1e-6; // Theta + u[1] = 2.23; // H + u[2] = 0.; // N + u[3] = blEdge->getVt(0); // ue + u[4] = 0.; // Ctau + + Ret[0] = u[3] * u[0] * Re; + if (Ret[0] < 1) + { + Ret[0] = 1.; + u[3] = Ret[0] / (u[0] * Re); + } + deltaStar[0] = u[0] * u[1]; + + std::vector<double> lamParam(8, 0.); + closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name); + + Hstar[0] = lamParam[0]; + Hstar2[0] = lamParam[1]; + Hk[0] = lamParam[2]; + cf[0] = lamParam[3]; + cd[0] = lamParam[4]; + delta[0] = lamParam[5]; + ctEq[0] = lamParam[6]; + us[0] = lamParam[7]; +} + +/** + * @brief Set the boundary condition at the first wake point. + * + * @param Re Freestream Reynolds number + * @param UpperCond Variables at the last point on the upper side. + * @param LowerCond Variables at the last point on the lower side. + */ +void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond) +{ + if (name != "wake") + std::cout << "Warning: Imposing wake boundary condition on " << name << std::endl; + u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low). + u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / u[0]; // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low). + u[2] = 9.; // Amplification factor. + u[3] = blEdge->getVt(0); // Inviscid velocity. + u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). + + Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness. + + if (Ret[0] < 1) + { + Ret[0] = 1; + u[3] = Ret[0] / (u[0] * Re); + } + deltaStar[0] = u[0] * u[1]; + + // Laminar closures. + std::vector<double> lamParam(8, 0.); + closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name); + Hstar[0] = lamParam[0]; + Hstar2[0] = lamParam[1]; + Hk[0] = lamParam[2]; + cf[0] = lamParam[3]; + cd[0] = lamParam[4]; + delta[0] = lamParam[5]; + ctEq[0] = lamParam[6]; + us[0] = lamParam[7]; + + for (size_t k = 0; k < mesh->getnMarkers(); ++k) + regime[k] = 1; +} + +/** + * @brief Compute the blowing velocity in each station of the region. + * + */ +void BoundaryLayer::computeBlwVel() +{ + deltaStar[0] = u[0] * u[1]; + blowingVelocity.resize(mesh->getnMarkers() - 1, 0.); + for (size_t iPoint = 1; iPoint < mesh->getnMarkers(); ++iPoint) + { + deltaStar[iPoint] = u[iPoint * nVar + 0] * u[iPoint * nVar + 1]; + blowingVelocity[iPoint - 1] = (blEdge->getRhoe(iPoint) * blEdge->getVt(iPoint) * deltaStar[iPoint] - blEdge->getRhoe(iPoint - 1) * blEdge->getVt(iPoint - 1) * deltaStar[iPoint - 1]) / (blEdge->getRhoe(iPoint) * (mesh->getLoc(iPoint) - mesh->getLoc(iPoint - 1))); + } +} + +/** + * @brief Print solution at a given station. + * + * @param iPoint Marker id. + * @note Function used for debugging. + */ +void BoundaryLayer::printSolution(size_t iPoint) const +{ + std::cout << "Pt " << iPoint << "Reg " << regime[iPoint] << std::endl; + std::cout << "x " << mesh->getx(iPoint) << "xx " << mesh->getLoc(iPoint) << "xxExt " << mesh->getExt(iPoint) << std::endl; + std::cout << std::scientific << std::setprecision(15); + std::cout << std::setw(10) << "T " << u[iPoint * nVar + 0] << "\n" + << std::setw(10) << "H " << u[iPoint * nVar + 1] << "\n" + << std::setw(10) << "N " << u[iPoint * nVar + 2] << "\n" + << std::setw(10) << "ue " << u[iPoint * nVar + 3] << "\n" + << std::setw(10) << "Ct " << u[iPoint * nVar + 4] << "\n" + << std::setw(10) << "dStar (BL) " << deltaStar[iPoint] << "\n" + << std::setw(10) << "dStar (old) " << blEdge->getDeltaStar(iPoint) << "\n" + << std::setw(10) << "vt " << blEdge->getVt(iPoint) << "\n" + << std::setw(10) << "Me " << blEdge->getMe(iPoint) << "\n" + << std::setw(10) << "rhoe " << blEdge->getRhoe(iPoint) << std::endl; +} \ No newline at end of file diff --git a/vii/src/DBoundaryLayer.h b/vii/src/DBoundaryLayer.h new file mode 100644 index 0000000..f6556ac --- /dev/null +++ b/vii/src/DBoundaryLayer.h @@ -0,0 +1,68 @@ +#ifndef DBoundaryLayer_H +#define DBoundaryLayer_H + +#include "vii.h" +#include "DClosures.h" +#include "DDiscretization.h" +#include "DEdge.h" + +namespace vii +{ + +/** + * @brief Boundary layer region upper/lower side or wake. + * @author Paul Dechamps + */ +class VII_API BoundaryLayer +{ + +private: + unsigned int const nVar = 5; ///< Number of variables of the partial differential problem. + double nCrit = 9.0; ///< Critical amplification factor. + +public: + Discretization *mesh; ///< 1D surface boundary layer mesh. + Edge *blEdge; ///< Quantites at the inviscid boundary (edge of the boundary layer). + Closures *closSolver; ///< Closure relations class. + std::string name; ///< Name of the region. + + // Boundary layer variables. + std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct). + std::vector<double> Ret; ///< Reynolds number based on the momentum thickness (theta). + std::vector<double> cd; ///< Local dissipation coefficient. + std::vector<double> cf; ///< Local friction coefficient. + std::vector<double> Hstar; ///< Kinetic energy shape parameter (thetaStar/theta). + std::vector<double> Hstar2; ///< Density shape parameter (deltaStarStar/theta). + std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)). + std::vector<double> ctEq; ///< Equilibrium shear stress coefficient (turbulent BL). + std::vector<double> us; ///< Equivalent normalized wall slip velocity. + std::vector<double> delta; ///< Boundary layer thickness. + std::vector<double> deltaStar; ///< Dispacement thickness (int(1-rho*u/rhoe*ue d_eta)). + + // Transition related variables. + std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime. + std::vector<double> blowingVelocity; ///< Blowing velocity. + size_t transMarker; ///< Marker id of the transition location. + double xtrF; ///< Forced transition location in the global frame of reference. + double xtr; ///< Transition location in the global frame of reference. + + BoundaryLayer(double _xtrF = -1.0); + ~BoundaryLayer(); + + // Boundary conditions. + void setStagBC(double Re); + void setWakeBC(double Re, std::vector<double> upperConditions, std::vector<double> lowerConditions); + void setMesh(std::vector<double> x, + std::vector<double> y, + std::vector<double> z); + + // Getters. + size_t getnVar() const { return nVar; }; + double getnCrit() const { return nCrit; }; + + // Others + void printSolution(size_t iPoint) const; + void computeBlwVel(); +}; +} // namespace vii +#endif // DBoundaryLayer_H diff --git a/vii/src/DClosures.cpp b/vii/src/DClosures.cpp new file mode 100644 index 0000000..badb021 --- /dev/null +++ b/vii/src/DClosures.cpp @@ -0,0 +1,256 @@ +#include "DClosures.h" +#include <cmath> +#include <iostream> + +using namespace vii; + +/** + * @brief Laminar closure relations at a given station. + * + * @param closureVars Vector where the computed variables will be stored. + * @param theta Momentum thickness. + * @param H Shape factor of the boundary layer. + * @param Ret Reynolds number based on the momentum thickness. + * @param Me Local Mach number. + * @param name Name of the region. + */ +void Closures::laminarClosures(std::vector<double> &closureVars, double theta, + double H, double Ret, const double Me, + const std::string name) +{ + H = std::max(H, 1.00005); + + // Kinematic shape factor H*. + double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + // Density shape parameter. + double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me; + + // Boundary layer thickness. + double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta); + + double Hstar = 0.; + double Hks = Hk - 4.35; + if (Hk <= 4.35) + Hstar = 0.0111 * Hks * Hks / (Hk + 1) - 0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 - 0.0002 * (Hks * Hk) * (Hks * Hk); + else + Hstar = 1.528 + 0.015 * Hks * Hks / Hk; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); + + // Friction coefficient. + double tmp = 0.; + double cf = 0.; + if (Hk < 5.5) + { + tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0); + cf = (-0.07 + 0.0727 * tmp) / Ret; + } + else if (Hk >= 5.5) + { + tmp = 1.0 - 1.0 / (Hk - 4.5); + cf = (-0.07 + 0.015 * tmp * tmp) / Ret; + } + + // Dissipation coefficient. + double Cd = 0.; + if (Hk < 4) + Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret)); + else + { + double HkCd = (Hk - 4) * (Hk - 4); + double denCd = 1 + 0.02 * HkCd; + Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret)); + } + + // Wake relations. + if (name == "wake") + { + Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2 * Ret)); + cf = 0.; + } + + double us = 0.; + double ctEq = 0.; + + closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; +} + +/** + * @brief Turbulent closure relations at a given station. + * + * @param closureVars Vector where the computed variables will be stored. + * @param theta Momentum thickness. + * @param H Shape factor of the boundary layer. + * @param Ct Shear stress coefficient. + * @param Ret Reynolds number based on the momentum thickness. + * @param Me Local Mach number. + * @param name Name of the region. + */ +void Closures::turbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name) +{ + H = std::max(H, 1.00005); + Ct = std::min(Ct, 0.3); + // Ct = std::max(std::min(Ct, 0.30), 0.0000001); + + double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me; + double gamma = 1.4 - 1.; + double Fc = std::sqrt(1 + 0.5 * gamma * Me * Me); + + double H0 = 0.; + if (Ret <= 400) + H0 = 4.0; + else + H0 = 3 + (400 / Ret); + if (Ret <= 200) + Ret = 200; + + double Hstar = 0.; + if (Hk <= H0) + Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret; + else + Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); + + double logRt = std::log(Ret / Fc); + logRt = std::max(logRt, 3.0); + double arg = std::max(-1.33 * Hk, -20.0); + + // Equivalent normalized wall slip velocity. + double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H)); + + // Boundary layer thickness. + double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta); + + double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); + double Hkc = 0.; + double Cdw = 0.; + double Cdd = 0.; + double Cdl = 0.; + + double Hmin = 0.; + double Fl = 0.; + double Dfac = 0.; + + double cf = 0.; + double Cd = 0.; + + double n = 1.0; + + double ctEq = 0.; + + if (name == "wake") + { + if (us > 0.99995) + us = 0.99995; + n = 2.0; // two wake halves + cf = 0.0; // no friction inside the wake + Hkc = Hk - 1; + Cdw = 0.0; // Wall contribution. + Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution. + Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; // Laminar stress contribution to outer layer. + ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5. + } + else + { + if (us > 0.95) + us = 0.98; + n = 1.0; + cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4 - (Hk / 0.875)) - 1)); + Hkc = std::max(Hk - 1 - 18 / Ret, 0.01); + // Dissipation coefficient. + Hmin = 1 + 2.1 / std::log(Ret); + Fl = (Hk - 1) / (Hmin - 1); + Dfac = 0.5 + 0.5 * std::tanh(Fl); + Cdw = 0.5 * (cf * us) * Dfac; + Cdd = (0.995 - us) * Ct * Ct; + Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; + ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + } + if (n * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0) + std::cout << "Negative sqrt encountered " << std::endl; + + // Dissipation coefficient. + Cd = n * (Cdw + Cdd + Cdl); + closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us}; +} + +/** + * @brief Turbulent closure relations at a given station. + * + * @param closureVars Computed equilibrium shear stress coefficient. + * @param theta Momentum thickness. + * @param H Shape factor of the boundary layer. + * @param Ct Shear stress coefficient. + * @param Ret Reynolds number based on the momentum thickness. + * @param Me Local Mach number. + * @param name Name of the region. + */ +void Closures::turbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name) +{ + H = std::max(H, 1.00005); + Ct = std::min(Ct, 0.3); + + double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); + if (name == "wake") + Hk = std::max(Hk, 1.00005); + else + Hk = std::max(Hk, 1.05000); + + double H0 = 0.; + if (Ret <= 400) + H0 = 4.0; + else + H0 = 3 + (400 / Ret); + if (Ret <= 200) + Ret = 200; + + double Hstar = 0.; + if (Hk <= H0) + Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret; + else + Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret; + + // Whitfield's minor additional compressibility correction. + Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); + + // Equivalent normalized wall slip velocity. + double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H)); + + double Ctcon = 0.5 / (6.7 * 6.7 * 0.75); + + double Hkc = 0.; + double ctEq = 0.; + + if (name == "wake") + { + if (us > 0.99995) + us = 0.99995; + Hkc = Hk - 1; + ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + } + else + { + if (us > 0.95) + us = 0.98; + Hkc = std::max(Hk - 1 - 18 / Ret, 0.01); + ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5 + } + if (Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0) + std::cout << "Negative sqrt encountered " << std::endl; + + closureVars = ctEq; +} \ No newline at end of file diff --git a/vii/src/DClosures.h b/vii/src/DClosures.h new file mode 100644 index 0000000..b08744f --- /dev/null +++ b/vii/src/DClosures.h @@ -0,0 +1,30 @@ +#ifndef DCLOSURES_H +#define DCLOSURES_H + +#include "vii.h" +#include <vector> +#include <string> + +namespace vii +{ + +/** + * @brief Boundary layer closure relations. + * @author Paul Dechamps + */ +class VII_API Closures +{ + +public: + static void laminarClosures(std::vector<double> &closureVars, double theta, + double H, double Ret, const double Me, + const std::string name); + static void turbulentClosures(std::vector<double> &closureVars, double theta, + double H, double Ct, double Ret, double Me, + const std::string name); + static void turbulentClosures(double &closureVars, double theta, + double H, double Ct, double Ret, const double Me, + const std::string name); +}; +} // namespace vii +#endif // DClosures_H \ No newline at end of file diff --git a/vii/src/DDiscretization.cpp b/vii/src/DDiscretization.cpp new file mode 100644 index 0000000..29b7033 --- /dev/null +++ b/vii/src/DDiscretization.cpp @@ -0,0 +1,54 @@ +#include "DDiscretization.h" +#include <cmath> + +using namespace vii; + +Discretization::Discretization() +{ + nMarkers = 0; + nElm = 0; + prevnMarkers = 0; +} + +Discretization::~Discretization() {} + +void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z) +{ + if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size()) + throw std::runtime_error("vii::Discretization Wrong mesh sizes.\n"); + + nMarkers = _x.size(); + x.resize(nMarkers, 0.); + y.resize(nMarkers, 0.); + z.resize(nMarkers, 0.); + nElm = nMarkers - 1; + x = _x; + y = _y; + z = _z; + + computeBLcoord(); +} + +/** + * @brief Compute nodes coordinates in the local frame of reference. + */ +void Discretization::computeBLcoord() +{ + loc.resize(nMarkers, 0.); + alpha.resize(nElm, 0.); + dx.resize(nElm, 0.); + + for (size_t iPoint = 0; iPoint < nElm; ++iPoint) + { + alpha[iPoint] = std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]); + // dx = sqrt(delta x ^2 + delta y ^2). + dx[iPoint] = std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint])); + loc[iPoint + 1] = loc[iPoint] + dx[iPoint]; + } +} + +void Discretization::setExt(std::vector<double> extM) +{ + ext.resize(extM.size(), 0.); + ext = extM; +}; diff --git a/vii/src/DDiscretization.h b/vii/src/DDiscretization.h new file mode 100644 index 0000000..f88480c --- /dev/null +++ b/vii/src/DDiscretization.h @@ -0,0 +1,49 @@ +#ifndef DDiscretization_H +#define DDiscretization_H + +#include "vii.h" +#include "DBoundaryLayer.h" + +namespace vii +{ + +/** + * @brief Boundary layer mesh. + * @author Paul Dechamps + */ +class VII_API Discretization +{ +private: + std::vector<double> x; ///< x coordinate in the global frame of reference. + std::vector<double> y; ///< y coordinate in the global frame of reference. + std::vector<double> z; ///< z coordinate in the global frame of reference. + std::vector<double> loc; ///< xi coordinate in the local frame of reference. + std::vector<double> ext; ///< xi coordinate in the local frame of reference, at the edge of the boundary layer. + std::vector<double> dx; ///< Cell size. + std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference. + size_t nMarkers; ///< Number of points in the domain. + size_t nElm; ///< Number of cells in the domain. + +public: + size_t prevnMarkers; ///< Number of points in the domain at the previous iteration. + + Discretization(); + ~Discretization(); + + // Getters. + size_t getnMarkers() const { return nMarkers; }; + double getLoc(size_t iPoint) const { return loc[iPoint]; }; + double getx(size_t iPoint) const { return x[iPoint]; }; + double gety(size_t iPoint) const { return y[iPoint]; }; + double getz(size_t iPoint) const { return z[iPoint]; }; + double getExt(size_t iPoint) const { return ext[iPoint]; }; + double getAlpha(size_t iElm) const { return alpha[iElm]; }; + size_t getDiffnElm() const { return nMarkers - prevnMarkers; }; + + // Setters. + void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z); + void setExt(std::vector<double> extM); + void computeBLcoord(); +}; +} // namespace vii +#endif // WDiscretization_H \ No newline at end of file diff --git a/vii/src/DDriver.cpp b/vii/src/DDriver.cpp new file mode 100644 index 0000000..a5d8b1a --- /dev/null +++ b/vii/src/DDriver.cpp @@ -0,0 +1,775 @@ +#include "DDriver.h" +#include "DBoundaryLayer.h" +#include "DSolver.h" +#include <iomanip> + +#define ANSI_COLOR_RED "\x1b[1;31m" +#define ANSI_COLOR_GREEN "\x1b[1;32m" +#define ANSI_COLOR_YELLOW "\x1b[1;33m" +#define ANSI_COLOR_RESET "\x1b[0m" + +using namespace vii; + +/** + * @brief Driver of the viscous solver. + * + * @param _Re Freestream Reynolds number. + * @param _Minf Freestream Mach number. + * @param _CFL0 Initial CFL number for pseudo-time integration. + * @param nSections Number of sections in the computation (1 for 2D cases). + * @param _Span Span of the considered wing (only used for drag computation). + * @param _verbose Verbosity level of the solver 0<=verbose<=1. + */ +Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _Span, unsigned int _verbose) +{ + // Freestream parameters. + Re = _Re; + CFL0 = _CFL0; + verbose = _verbose; + + // Initialzing boundary layer + sections.resize(nSections); + convergenceStatus.resize(nSections); + for (size_t iSec = 0; iSec < nSections; ++iSec) + { + convergenceStatus[iSec].resize(3); + for (size_t i = 0; i < 3; ++i) + { + BoundaryLayer *region = new BoundaryLayer(); + sections[iSec].push_back(region); + } + sections[iSec][0]->name = "upper"; + sections[iSec][1]->name = "lower"; + sections[iSec][2]->name = "wake"; + } + + Cdt_sec.resize(sections.size(), 0.); + Cdf_sec.resize(sections.size(), 0.); + Cdp_sec.resize(sections.size(), 0.); + Cdt = 0.; + Cdf = 0.; + Cdp = 0.; + span = _Span; + + // Pre/post processing flags. + stagPtMvmt.resize(sections.size(), false); + + // Time solver initialization. + tSolver = new Solver(_CFL0, _Minf, Re); + + // Numerical parameters. + std::cout << "--- Viscous solver setup ---\n" + << "Reynolds number: " << Re << " \n" + << "Freestream Mach number: " << _Minf << " \n" + << "Initial CFL number: " << CFL0 << " \n" + << std::endl; + std::cout << "Boundary layer initialized." << std::endl; +} + +/** + * @brief Driver and subsequent dependencies destruction. + */ +Driver::~Driver() +{ + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + for (size_t i = 0; i < sections[iSec].size(); ++i) + delete sections[iSec][i]; + delete tSolver; + std::cout << "~vii::Driver()\n"; +} // end ~Driver + +/** + * @brief Runs viscous operations. Temporal solver parametrisation is first adapted, + * Solutions are initialized than time marched towards steady state successively for each points. + * + * @param couplIter Current coupling iteration. + * @return int Solver exit code. + */ +int Driver::run(unsigned int const couplIter) +{ + bool lockTrans = false; // Flagging activation of transition routines. + int pointExitCode; // Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter). + + double maxMach = 0.; + + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + { + tms["0-CheckIn"].start(); + if (verbose > 1) + std::cout << "Sec " << iSec << " "; + + // Check for stagnation point movement. + if (couplIter != 0 && (sections[iSec][0]->mesh->getDiffnElm())) + { + stagPtMvmt[iSec] = true; + if (verbose > 2) + std::cout << "Stagnation point moved by " << sections[iSec][0]->mesh->getDiffnElm() << " elements." << std::endl; + } + else + stagPtMvmt[iSec] = false; + + // Set boundary conditions. + sections[iSec][0]->xtr = 1.; // Upper side initial xtr. + sections[iSec][1]->xtr = 1.; // Lower side initial xtr. + sections[iSec][2]->xtr = 0.; // Wake initial xtr (fully turbulent). + sections[iSec][0]->setStagBC(Re); + sections[iSec][1]->setStagBC(Re); + tms["0-CheckIn"].stop(); + + // Loop over the regions (upper, lower and wake). + for (size_t iRegion = 0; iRegion < sections[iSec].size(); ++iRegion) + { + convergenceStatus[iSec][iRegion].resize(0); + + // Maximum Mach number in the region. + if (sections[iSec][iRegion]->blEdge->getMaxMach() > maxMach) + maxMach = sections[iSec][iRegion]->blEdge->getMaxMach(); + + if (verbose > 1) + std::cout << sections[iSec][iRegion]->name << " "; + + // Check if safe mode is required. + if (couplIter == 0 || sections[iSec][iRegion]->mesh->getDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0) + { + tSolver->setCFL0(1); + tSolver->setitFrozenJac(1); + tSolver->setinitSln(true); + + if (sections[iSec][iRegion]->mesh->getDiffnElm()) + { + std::vector<double> xxExtCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.); + for (size_t iPoint = 0; iPoint < xxExtCurr.size(); ++iPoint) + xxExtCurr[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint); + sections[iSec][iRegion]->mesh->setExt(xxExtCurr); + + std::vector<double> DeltaStarZeros(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0); + sections[iSec][iRegion]->blEdge->setDeltaStar(DeltaStarZeros); + + std::vector<double> ueOld(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0); + for (size_t iPoint = 0; iPoint < ueOld.size(); ++iPoint) + ueOld[iPoint] = sections[iSec][iRegion]->blEdge->getVt(iPoint); + sections[iSec][iRegion]->blEdge->setVtOld(ueOld); + } + + // Keyword annoncing iteration safety level. + if (verbose > 1) + std::cout << "restart. "; + } + else + { + tSolver->setCFL0(CFL0); + tSolver->setitFrozenJac(5); + tSolver->setinitSln(false); + + // Keyword annoncing iteration safety level. + if (verbose > 1) + std::cout << "update. "; + } + + // Save number of points in each regions. + sections[iSec][iRegion]->mesh->prevnMarkers = sections[iSec][iRegion]->mesh->getnMarkers(); + + // Reset regime. + lockTrans = false; + if (iRegion < 2) // Airfoil + for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++) + sections[iSec][iRegion]->regime[k] = 0; + else // Wake + { + for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++) + sections[iSec][iRegion]->regime[k] = 1; + + // Impose wake boundary condition. + std::vector<double> upperConditions(sections[iSec][iRegion]->getnVar(), 0.); + std::vector<double> lowerConditions(sections[iSec][iRegion]->getnVar(), 0.); + for (size_t k = 0; k < upperConditions.size(); ++k) + { + upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->mesh->getnMarkers() - 1) * sections[iSec][0]->getnVar() + k]; + lowerConditions[k] = sections[iSec][0]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k]; + } + sections[iSec][iRegion]->setWakeBC(Re, upperConditions, lowerConditions); + lockTrans = true; + } + + // Loop over points + for (size_t iPoint = 1; iPoint < sections[iSec][iRegion]->mesh->getnMarkers(); ++iPoint) + { + // Initialize solution at point + tSolver->initialCondition(iPoint, sections[iSec][iRegion]); + // Solve equations + tms["1-Solver"].start(); + pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]); + // sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.; + tms["1-Solver"].stop(); + + if (pointExitCode != 0) + convergenceStatus[iSec][iRegion].push_back(iPoint); + + // Transition + tms["2-Transition"].start(); + if (!lockTrans) + { + // Check if perturbation waves are growing or not. + checkWaves(iPoint, sections[iSec][iRegion]); + + // Amplification factor is compared to critical amplification factor. + if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit()) + { + averageTransition(iPoint, sections[iSec][iRegion], 0); + if (verbose > 1) + { + std::cout << std::fixed; + std::cout << std::setprecision(2); + std::cout << sections[iSec][iRegion]->xtr + << " (" << sections[iSec][iRegion]->transMarker << ") "; + } + lockTrans = true; + } + } + tms["2-Transition"].stop(); + } + tms["3-Blowing"].start(); + sections[iSec][iRegion]->computeBlwVel(); + tms["3-Blowing"].stop(); + } + if (verbose > 1) + std::cout << std::endl; + } + + for (size_t i = 0; i < sections.size(); ++i) + computeSectionalDrag(sections[i], i); + computeTotalDrag(); + + // Output information. + solverExitCode = outputStatus(maxMach); + + return solverExitCode; // exit code +} + +/** + * @brief Check whether or not instabilities are growing in the laminar boundary layer. + * + * @param iPoint Marker id. + * @param bl BoundaryLayer region. + */ +void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl) +{ + + double logRet_crit = 2.492 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0); + + if (bl->Ret[iPoint] > 0) // Avoid log of negative number. + { + if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08) + bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow. + } + else + bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow. +} + +/** + * @brief Computes turbulent solution @ the transition point and averages solutions. + * + * @param iPoint Marker id. + * @param bl BoundaryLayer region. + * @param forced Flag for forced transition. + */ +void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced) +{ + // Averages solution on transition marker. + size_t nVar = bl->getnVar(); + + // Save laminar solution. + std::vector<double> lamSol(nVar, 0.0); + for (size_t k = 0; k < lamSol.size(); ++k) + lamSol[k] = bl->u[iPoint * nVar + k]; + + // Set up turbulent portion boundary condition. + bl->transMarker = iPoint; // Save transition marker. + + // Loop over remaining points in the region. + for (size_t k = iPoint; k < bl->mesh->getnMarkers(); ++k) + bl->regime[k] = 1; // Set Turbulent regime. + + // Compute transition location. + bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1); + + // Percentage of the subsection corresponding to a turbulent flow. + double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)); + double avLam = 1. - avTurb; + + // Impose boundary condition. + double Cteq_trans; + bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name); + bl->ctEq[iPoint - 1] = Cteq_trans; + bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1]; + + // Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above). + // These initial guess values do not influence the converged solution. + bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point). + bl->u[iPoint * nVar + 1] = 1.515; // Because we expect a significant drop of the shape factor. + + // Solve point in turbulent configuration. + int exitCode = tSolver->integration(iPoint, bl); + + if (exitCode != 0 && verbose > 1) + std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl; + + // Average both solutions (Now solution @ iPoint is the turbulent solution). + bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta. + bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H. + bl->u[iPoint * nVar + 2] = 9.; // N. + bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue. + bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct. + + // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name); + bl->Hstar[iPoint - 1] = lamParam[0]; + bl->Hstar2[iPoint - 1] = lamParam[1]; + bl->Hk[iPoint - 1] = lamParam[2]; + bl->cf[iPoint - 1] = lamParam[3]; + bl->cd[iPoint - 1] = lamParam[4]; + bl->delta[iPoint - 1] = lamParam[5]; + bl->ctEq[iPoint - 1] = lamParam[6]; + bl->us[iPoint - 1] = lamParam[7]; + + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures(turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; +} + +/** + * @brief Friction, pressure and total drag computation qt each station of the configuration. + * + * @tparam Cdt = theta * Ue^((H+5)/2). + * @tparam Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower. + * @tparam Cdp = Cdtot - Cdf. + * + * @param bl BoundaryLayer region. + * @param i Considered section. + */ +void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i) +{ + unsigned int nVar = bl[0]->getnVar(); + size_t lastWkPt = (bl[2]->mesh->getnMarkers() - 1) * nVar; + + // Normalize friction and dissipation coefficients. + std::vector<double> frictioncoeff_upper(bl[0]->mesh->getnMarkers(), 0.0); + for (size_t k = 0; k < frictioncoeff_upper.size(); ++k) + frictioncoeff_upper[k] = bl[0]->blEdge->getVt(k) * bl[0]->blEdge->getVt(k) * bl[0]->cf[k]; + + std::vector<double> frictioncoeff_lower(bl[1]->mesh->getnMarkers(), 0.0); + for (size_t k = 0; k < frictioncoeff_lower.size(); ++k) + frictioncoeff_lower[k] = bl[1]->blEdge->getVt(k) * bl[1]->blEdge->getVt(k) * bl[1]->cf[k]; + + // Compute friction drag coefficient. + double Cdf_upper = 0.0; + for (size_t k = 1; k < bl[0]->mesh->getnMarkers(); ++k) + Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->getx(k) - bl[0]->mesh->getx(k - 1)) / 2; + double Cdf_lower = 0.0; + for (size_t k = 1; k < bl[1]->mesh->getnMarkers(); ++k) + Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->getx(k) - bl[1]->mesh->getx(k - 1)) / 2; + + // Friction drag coefficient. + Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake + // Total drag coefficient. + Cdt_sec[i] = 2 * bl[2]->u[lastWkPt + 0] * pow(bl[2]->u[lastWkPt + 3], (bl[2]->u[lastWkPt + 1] + 5) / 2); + // Pressure drag coefficient. + Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i]; +} + +/** + * @brief Compute total drag coefficient on the airfoil/wing. + * Performs the sectional drag integration for 3D computations. + * + */ +void Driver::computeTotalDrag() +{ + Cdt = 0.; + Cdf = 0.; + Cdp = 0.; + + if (sections.size() == 1 || span == 0.) + { + Cdt = Cdt_sec[0]; + Cdf = Cdf_sec[0]; + Cdp = Cdp_sec[0]; + } + else + { + Cdt += (sections[0][0]->mesh->getz(0) - 0) * Cdt_sec[0]; + Cdp += (sections[0][0]->mesh->getz(0) - 0) * Cdp_sec[0]; + Cdf += (sections[0][0]->mesh->getz(0) - 0) * Cdf_sec[0]; + double dz = 0.; + for (size_t k = 1; k < sections.size(); ++k) + { + dz = (sections[k][0]->mesh->getz(0) - sections[k - 1][0]->mesh->getz(0)); + Cdt += dz * Cdt_sec[k]; + Cdp += dz * Cdp_sec[k]; + Cdf += dz * Cdf_sec[k]; + } + Cdt /= span; + Cdp /= span; + Cdf /= span; + } +} + +/** + * @brief Outputs solver status depending on the verbosity level (verbose). + * + * @param maxMach Maximum Mach number identified on the configuration. + * @return int Solver exit code (0 converged, !=0 not converged). + */ +int Driver::outputStatus(double maxMach) const +{ + int solverExitCode = 0; + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) + if (convergenceStatus[iSec][iReg].size() != 0) + { + solverExitCode = -1; + break; + } + + if (verbose > 2 && solverExitCode != 0) + { + // Output unconverged points. + std::cout << "Unconverged points" << std::endl; + std::cout << std::setw(10) << "Section" + << std::setw(12) << "Region" + << std::setw(16) << "Markers" << std::endl; + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + { + unsigned int count = 0; + for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg) + { + std::cout << std::setw(10) << iSec; + std::cout << std::setw(12) << iReg << " "; + if (convergenceStatus[iSec][iReg].size() != 0) + for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); ++iPoint) + { + std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" << sections[iSec][iReg]->mesh->getx(convergenceStatus[iSec][iReg][iPoint]) << "), "; + ++count; + } + std::cout << "\n"; + } + std::cout << "total: " << count << "\n"; + } + std::cout << "" << std::endl; + } + + if (verbose > 0) + { + // Compute mean transition locations. + std::vector<double> meanTransitions(2, 0.); + for (size_t iSec = 0; iSec < sections.size(); ++iSec) + for (size_t iReg = 0; iReg < 2; ++iReg) + meanTransitions[iReg] += sections[iSec][iReg]->xtr; + for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg) + meanTransitions[iReg] /= sections.size(); + + std::cout << std::fixed << std::setprecision(2) + << "Viscous solver for Re " << Re / 1e6 << "e6, Mach max " + << maxMach + << std::endl; + + std::cout << std::setw(10) << "xTr Upper" + << std::setw(12) << "xTr Lower" + << std::setw(12) << "Cd" << std::endl; + std::cout << std::fixed << std::setprecision(2); + std::cout << std::setw(10) << meanTransitions[0] + << std::setw(12) << meanTransitions[1]; + std::cout << std::fixed << std::setprecision(5); + std::cout << std::setw(12) << Cdt << "\n"; + + if (verbose > 0) + { + if (solverExitCode == 0) + std::cout << ANSI_COLOR_GREEN << "Viscous solver converged." << ANSI_COLOR_RESET << std::endl; + else + std::cout << ANSI_COLOR_YELLOW << "Viscous solver not converged." << ANSI_COLOR_RESET << std::endl; + } + + if (verbose > 2) + { + std::cout << "Viscous solver CPU" << std::endl + << tms; + } + } + + if (verbose > 0) + std::cout << "" << std::endl; + return solverExitCode; +} + +/** + * @brief Set nodes coordinates in the global frame of reference. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _x Vector of x coordinates of the nodes in the region. + * @param _y Vector of y coordinates of the nodes in the region. + * @param _z Vector of z coordinates of the nodes in the region. + */ +void Driver::setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z) +{ + sections[iSec][iRegion]->setMesh(_x, _y, _z); +} + +/** + * @brief Set the velocity at the nodes in the global frame of reference. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param vx Vector of velocities projected on the x axis of the global frame of reference. + * @param vy Vector of velocities projected on the y axis of the global frame of reference. + * @param vz Vector of velocities projected on the z axis of the global frame of reference. + */ +void Driver::setVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz) +{ + if (vx.size() != vy.size() || vy.size() != vz.size() || vx.size() != vz.size()) + throw std::runtime_error("dart::Driver Wrong velocity vector sizes."); + if (vx.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + throw std::runtime_error("dart::Driver Velocity vector is not coherent with the mesh."); + + std::vector<double> vt(vx.size(), 0.0); + + for (size_t iPoint = 0; iPoint < vx.size() - 1; ++iPoint) + { + vt[iPoint] = vx[iPoint] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint)); + vt[iPoint + 1] = vx[iPoint + 1] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint + 1] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint)); + } + + sections[iSec][iRegion]->blEdge->setVt(vt); +} + +/** + * @brief Set the Mach number at the nodes. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _Me Vector of Mach numbers at the nodes. + */ +void Driver::setMe(size_t iSec, size_t iRegion, std::vector<double> _Me) +{ + if (_Me.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + throw std::runtime_error("dart::Driver Mach number vector is not coherent with the mesh."); + sections[iSec][iRegion]->blEdge->setMe(_Me); +} + +/** + * @brief Set the density at the nodes. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _Rhoe Vector of density at the nodes. + */ +void Driver::setRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe) +{ + if (_Rhoe.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + throw std::runtime_error("dart::Driver Density vector is not coherent with the mesh."); + sections[iSec][iRegion]->blEdge->setRhoe(_Rhoe); +} + +/** + * @brief Set the nodes coordinates in the local frame of reference at the edge of the boundary layer. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _xxExt Vector of coordinates of the nodes. + */ +void Driver::setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt) +{ + if (_xxExt.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + throw std::runtime_error("dart::Driver External mesh vector is not coherent with the mesh."); + sections[iSec][iRegion]->mesh->setExt(_xxExt); +} + +/** + * @brief Set the displacement thickness (previous iteration). + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _DeltaStarExt Vector of displacement thicknesses at the nodes. + */ +void Driver::setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt) +{ + if (_DeltaStarExt.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + { + std::cout << "size DeltaStar: " << _DeltaStarExt.size() << ". nMarkers: " << sections[iSec][iRegion]->mesh->getnMarkers() << std::endl; + throw std::runtime_error("dart::Driver External delta star vector is not coherent with the mesh."); + } + sections[iSec][iRegion]->blEdge->setDeltaStar(_DeltaStarExt); +} + +/** + * @brief Returns x coordinates of the nodes in the local frame of reference in the considered region. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @return std::vector<double> + */ +std::vector<double> Driver::extractxCoord(size_t iSec, size_t iRegion) const +{ + std::vector<double> xCoord(sections[iSec][iRegion]->mesh->getnMarkers(), 0.); + for (size_t iPoint = 0; iPoint < xCoord.size(); ++iPoint) + xCoord[iPoint] = sections[iSec][iRegion]->mesh->getx(iPoint); + return xCoord; +} + +/** + * @brief Returns blowing velocities (defined on the cells' cg) in the considered region. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @return std::vector<double> + */ +std::vector<double> Driver::extractBlowingVel(size_t iSec, size_t iRegion) const +{ + return sections[iSec][iRegion]->blowingVelocity; +} +std::vector<double> Driver::extractxx(size_t iSec, size_t iRegion) const +{ + std::vector<double> xx(sections[iSec][iRegion]->mesh->getnMarkers(), 0.); + for (size_t iPoint = 0; iPoint < xx.size(); ++iPoint) + xx[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint); + return xx; +} + +/** + * @brief Returns the displacement thickness in the considered region. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @return std::vector<double> + */ +std::vector<double> Driver::extractDeltaStar(size_t iSec, size_t iRegion) const +{ + return sections[iSec][iRegion]->deltaStar; +} + +/** + * @brief Returns the velocity of the interaction law in the considered region. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @return std::vector<double> + */ +std::vector<double> Driver::extractUe(size_t iSec, size_t iRegion) const +{ + std::vector<double> ueCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0); + for (size_t i = 0; i < ueCurr.size(); ++i) + ueCurr[i] = sections[iSec][iRegion]->u[i * sections[iSec][iRegion]->getnVar() + 3]; + return ueCurr; +} + +/** + * @brief Set the velocity at the edge of the BL at the previous iteration in the considered region. + * + * @param iSec id of the considered section. + * @param iRegion id of the considered region. + * @param _ue Velocity vector. + * @return std::vector<double> + */ +void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue) +{ + if (_ue.size() != sections[iSec][iRegion]->mesh->getnMarkers()) + throw std::runtime_error("dart::Driver External velocity vector is not coherent with the mesh."); + sections[iSec][iRegion]->blEdge->setVtOld(_ue); +} + +/** + * @brief Returns the integral boundary layer solution in a section (sorted from upper Te to lower Te than wake). + * + * @param iSec id of the considered section. + * @return std::vector<std::vector<double>> + */ +std::vector<std::vector<double>> Driver::getSolution(size_t iSec) +{ + size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers(); + size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers() - 1; // No stagnation point doublon. + size_t nVar = sections[iSec][0]->getnVar(); + + std::vector<std::vector<double>> Sln(16); + for (size_t iVector = 0; iVector < Sln.size(); ++iVector) + Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.); + // Upper side. + size_t revPoint = nMarkersUp - 1; + for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint) + { + Sln[0][iPoint] = sections[iSec][0]->u[revPoint * nVar + 0]; + Sln[1][iPoint] = sections[iSec][0]->u[revPoint * nVar + 1]; + Sln[2][iPoint] = sections[iSec][0]->u[revPoint * nVar + 2]; + Sln[3][iPoint] = sections[iSec][0]->u[revPoint * nVar + 3]; + Sln[4][iPoint] = sections[iSec][0]->u[revPoint * nVar + 4]; + Sln[5][iPoint] = sections[iSec][0]->deltaStar[revPoint]; + + Sln[6][iPoint] = sections[iSec][0]->Ret[revPoint]; + Sln[7][iPoint] = sections[iSec][0]->cd[revPoint]; + Sln[8][iPoint] = sections[iSec][0]->cf[revPoint]; + Sln[9][iPoint] = sections[iSec][0]->Hstar[revPoint]; + Sln[10][iPoint] = sections[iSec][0]->Hstar2[revPoint]; + Sln[11][iPoint] = sections[iSec][0]->Hk[revPoint]; + Sln[12][iPoint] = sections[iSec][0]->ctEq[revPoint]; + Sln[13][iPoint] = sections[iSec][0]->us[revPoint]; + Sln[14][iPoint] = sections[iSec][0]->delta[revPoint]; + + Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint); + --revPoint; + } + // Lower side. + for (size_t iPoint = 0; iPoint < sections[iSec][1]->mesh->getnMarkers(); ++iPoint) + { + Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0]; + Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1]; + Sln[2][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 2]; + Sln[3][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 3]; + Sln[4][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 4]; + Sln[5][nMarkersUp + iPoint] = sections[iSec][1]->deltaStar[iPoint]; + + Sln[6][nMarkersUp + iPoint] = sections[iSec][1]->Ret[iPoint]; + Sln[7][nMarkersUp + iPoint] = sections[iSec][1]->cd[iPoint]; + Sln[8][nMarkersUp + iPoint] = sections[iSec][1]->cf[iPoint]; + Sln[9][nMarkersUp + iPoint] = sections[iSec][1]->Hstar[iPoint]; + Sln[10][nMarkersUp + iPoint] = sections[iSec][1]->Hstar2[iPoint]; + Sln[11][nMarkersUp + iPoint] = sections[iSec][1]->Hk[iPoint]; + Sln[12][nMarkersUp + iPoint] = sections[iSec][1]->ctEq[iPoint]; + Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint]; + Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint]; + + Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint); + } + // Wake. + for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint) + { + Sln[0][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 0]; + Sln[1][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 1]; + Sln[2][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 2]; + Sln[3][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 3]; + Sln[4][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 4]; + Sln[5][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->deltaStar[iPoint]; + + Sln[6][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Ret[iPoint]; + Sln[7][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cd[iPoint]; + Sln[8][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cf[iPoint]; + Sln[9][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar[iPoint]; + Sln[10][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar2[iPoint]; + Sln[11][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hk[iPoint]; + Sln[12][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->ctEq[iPoint]; + Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint]; + Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint]; + + Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint); + } + + if (verbose > 2) + std::cout << "Solution structure sent." << std::endl; + return Sln; +} \ No newline at end of file diff --git a/vii/src/DDriver.h b/vii/src/DDriver.h new file mode 100644 index 0000000..0ba4854 --- /dev/null +++ b/vii/src/DDriver.h @@ -0,0 +1,74 @@ +#ifndef DDriver_H +#define DDriver_H + +#include "vii.h" +#include "wObject.h" +#include "DBoundaryLayer.h" +#include "wTimers.h" +namespace vii +{ + +/** + * @brief Boundary layer solver driver. Performs the space marching and set up time marching for each point. + * @authors Paul Dechamps + */ +class VII_API Driver : public fwk::wSharedObject +{ +public: + double Cdt; ///< Total drag coefficient. + double Cdf; ///< Total friction drag coefficient. + double Cdp; ///< Total pressure drag coefficient. + std::vector<double> Cdt_sec; ///< Sectional total drag coefficient. + std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient. + std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient. + std::vector<std::vector<BoundaryLayer *>> sections; ///< Boundary layer regions sorted by section. + unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1. + +private: + double Re; ///< Freestream Reynolds number. + double CFL0; ///< Initial CFL number for pseudo time integration. + double span; ///< Wing Span (Used only for drag computation, not used if 2D case). + std::vector<bool> stagPtMvmt; ///< Flag to account for stagnation point movements. + Solver *tSolver; ///< Pseudo-time solver. + int solverExitCode; ///< Exit code of viscous calculations. + std::vector<std::vector<std::vector<size_t>>> convergenceStatus; ///< Vector containing points that did not converged. + +protected: + fwk::Timers tms; ///< Internal timers + +public: + Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _span = 0., unsigned int verbose = 1); + ~Driver(); + int run(unsigned int const couplIter); + + // Getters. + double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; }; + double getRe() const { return Re; }; + std::vector<std::vector<double>> getSolution(size_t iSec); + + // Setters. + void setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z); + void setVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz); + void setMe(size_t iSec, size_t iRegion, std::vector<double> _Me); + void setRhoe(size_t iSec, size_t iRegion, std::vector<double> _rhoe); + void setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt); + void setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _deltaStarExt); + void setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue); + + // Extractors. + std::vector<double> extractBlowingVel(size_t iSec, size_t iRegion) const; + std::vector<double> extractxx(size_t iSec, size_t iRegion) const; + std::vector<double> extractDeltaStar(size_t iSec, size_t iRegion) const; + std::vector<double> extractxCoord(size_t iSec, size_t iRegion) const; + std::vector<double> extractUe(size_t iSec, size_t iRegion) const; + +private: + void checkWaves(size_t iPoint, BoundaryLayer *bl); + void averageTransition(size_t iPoint, BoundaryLayer *bl, int forced); + void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i); + void computeTotalDrag(); + void computeBlwVel(); + int outputStatus(double maxMach) const; +}; +} // namespace vii +#endif // DDriver_H \ No newline at end of file diff --git a/vii/src/DEdge.cpp b/vii/src/DEdge.cpp new file mode 100644 index 0000000..48bde11 --- /dev/null +++ b/vii/src/DEdge.cpp @@ -0,0 +1,59 @@ +#include "DEdge.h" + +using namespace vii; + +/** + * @brief Construct a new Edge::Edge object. + * + */ +Edge::Edge() +{ + Me.resize(0, 0.); + vt.resize(0, 0.); + rhoe.resize(0, 0.); + deltaStar.resize(0, 0.); + vtOld.resize(0, 0.); +} + +/** + * @brief Return the maximum inviscid Mach number in the region. + * + * @return double + */ +double Edge::getMaxMach() const +{ + if (Me.size() <= 0) + throw std::runtime_error("vii::Edge Invalid Mach number vector."); + double Mmax = 0.0; + for (size_t i = 0; i < Me.size(); ++i) + if (Me[i] > Mmax) + Mmax = Me[i]; + return Mmax; +} + +// Setters. +void Edge::setMe(std::vector<double> const _Me) +{ + Me.resize(_Me.size(), 0.); + Me = _Me; +} +void Edge::setVt(std::vector<double> const _vt) +{ + vt.resize(_vt.size(), 0.); + vt = _vt; +} +void Edge::setRhoe(std::vector<double> const _rhoe) +{ + rhoe.resize(_rhoe.size(), 0.); + rhoe = _rhoe; +} +void Edge::setDeltaStar(std::vector<double> const _deltaStar) +{ + deltaStar.resize(_deltaStar.size(), 0.); + deltaStar = _deltaStar; +} +void Edge::setVtOld(std::vector<double> const _vtOld) +{ + vtOld.resize(_vtOld.size(), 0.); + vtOld = _vtOld; +} diff --git a/vii/src/DEdge.h b/vii/src/DEdge.h new file mode 100644 index 0000000..0231139 --- /dev/null +++ b/vii/src/DEdge.h @@ -0,0 +1,44 @@ +#ifndef DEdge_H +#define DEdge_H + +#include "vii.h" +#include <vector> +#include <iostream> + +namespace vii +{ + +/** + * @brief Inviscid quantities at the edge of the boundary layer of a BoundaryLayer region. + * @author Paul Dechamps + */ +class VII_API Edge +{ +private: + std::vector<double> Me; ///< Mach number. + std::vector<double> vt; ///< Velocity. + std::vector<double> rhoe; ///< Density. + std::vector<double> deltaStar; ///< Displacement thickness. + std::vector<double> vtOld; ///< Previous velocity. + +public: + Edge(); + ~Edge(){}; + + // Getters. + double getMe(size_t iPoint) const { return Me[iPoint]; }; + double getVt(size_t iPoint) const { return vt[iPoint]; }; + double getRhoe(size_t iPoint) const { return rhoe[iPoint]; }; + double getDeltaStar(size_t iPoint) const { return deltaStar[iPoint]; }; + double getVtOld(size_t iPoint) const { return vtOld[iPoint]; }; + double getMaxMach() const; + + // Setters. + void setMe(std::vector<double> const _Me); + void setVt(std::vector<double> const _vt); + void setRhoe(std::vector<double> const _rhoe); + void setDeltaStar(std::vector<double> const _deltaStar); + void setVtOld(std::vector<double> const _vtOld); +}; +} // namespace vii +#endif // DEdge_H \ No newline at end of file diff --git a/vii/src/DFluxes.cpp b/vii/src/DFluxes.cpp new file mode 100644 index 0000000..66ff922 --- /dev/null +++ b/vii/src/DFluxes.cpp @@ -0,0 +1,242 @@ +#include "DFluxes.h" +#include "DBoundaryLayer.h" +#include <Eigen/Dense> + +using namespace vii; +using namespace Eigen; + +/** + * @brief Compute residual vector of the integral boundary layer equations. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + * @return VectorXd + */ +VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + std::vector<double> u(nVar, 0.); + for (size_t k = 0; k < nVar; ++k) + u[k] = bl->u[iPoint * nVar + k]; + return blLaws(iPoint, bl, u); +} + +/** + * @brief Compute Jacobian of the integral boundary layer equations. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + * @param sysRes Residual of the IBL at point iPoint. + * @param eps Pertubation from current solution used to compute the Jacobian. + * @return MatrixXd + */ +MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, VectorXd const &sysRes, double eps) +{ + size_t nVar = bl->getnVar(); + MatrixXd JacMatrix = MatrixXd::Zero(5, 5); + std::vector<double> uUp(nVar, 0.); + for (size_t k = 0; k < nVar; ++k) + uUp[k] = bl->u[iPoint * nVar + k]; + + double varSave = 0.; + for (size_t iVar = 0; iVar < nVar; ++iVar) + { + varSave = uUp[iVar]; + uUp[iVar] += eps; + JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps; + uUp[iVar] = varSave; + } + return JacMatrix; +} + +/** + * @brief Integral boundary layer equation. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + * @param u Solution vector at point iPoint. + * @return VectorXd + */ +VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u) +{ + size_t nVar = bl->getnVar(); + + MatrixXd timeMatrix = MatrixXd::Zero(5, 5); + VectorXd spaceVector = VectorXd::Zero(5); + + // Dissipation ratio. + double dissipRatio; + if (bl->name == "wake") + dissipRatio = 0.9; + else + dissipRatio = 1.; + + bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta. + bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness. + + // Boundary layer closure relations. + if (bl->regime[iPoint] == 0) + { + // Laminar closure relations. + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name); + bl->Hstar[iPoint] = lamParam[0]; + bl->Hstar2[iPoint] = lamParam[1]; + bl->Hk[iPoint] = lamParam[2]; + bl->cf[iPoint] = lamParam[3]; + bl->cd[iPoint] = lamParam[4]; + bl->delta[iPoint] = lamParam[5]; + bl->ctEq[iPoint] = lamParam[6]; + bl->us[iPoint] = lamParam[7]; + } + + else if (bl->regime[iPoint] == 1) + { + // Turbulent closure relations. + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; + } + else + { + std::cout << "Wrong regime\n" + << std::endl; + throw; + } + + // Gradients computation. + double dx = (bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1)); + double dxExt = (bl->mesh->getExt(iPoint) - bl->mesh->getExt(iPoint - 1)); + double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx; + double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx; + double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx; + double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx; + double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx; + + double dueExt_dx = (bl->blEdge->getVt(iPoint) - bl->blEdge->getVt(iPoint - 1)) / dxExt; + double ddeltaStarExt_dx = (bl->blEdge->getDeltaStar(iPoint) - bl->blEdge->getDeltaStar(iPoint - 1)) / dxExt; + + double c = 4 / (M_PI * dx); + double cExt = 4 / (M_PI * dxExt); + + // Amplification routine. + double dN_dx = 0.0; + double ax = 0.0; + if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0) + { + // No forced transition and laminar regime. + ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]); + dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx; + } + + double Mea = bl->blEdge->getMe(iPoint); + double Hstara = bl->Hstar[iPoint]; + double Hstar2a = bl->Hstar2[iPoint]; + double Hka = bl->Hk[iPoint]; + double cfa = bl->cf[iPoint]; + double cda = bl->cd[iPoint]; + double deltaa = bl->delta[iPoint]; + double ctEqa = bl->ctEq[iPoint]; + double usa = bl->us[iPoint]; + + // Space part. + spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2; + spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2; + spaceVector(2) = dN_dx - ax; + spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx; + + if (bl->regime[iPoint] == 1) + spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx); + + // Time part. + timeMatrix(0, 0) = u[1] / u[3]; + timeMatrix(0, 1) = u[0] / u[3]; + timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]); + + timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3]; + timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3]; + timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]); + timeMatrix(2, 2) = 1.; + + timeMatrix(3, 0) = -c * u[1]; + timeMatrix(3, 1) = -c * u[0]; + timeMatrix(3, 3) = 1.; + + if (bl->regime[iPoint] == 1) + { + timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]); + timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa); + } + else + timeMatrix(4, 4) = 1.0; + + return timeMatrix.partialPivLu().solve(spaceVector); +} + +/** + * @brief Amplification routines of the e^N method for transition capturing. + * + * @param Hk Kinematic shape parameter. + * @param Ret Reynolds number based on the momentum thickness. + * @param theta Momentum thickness. + * @return double + */ +double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const +{ + double dgr = 0.08; + double Hk1 = 3.5; + double Hk2 = 4.; + double Hmi = (1 / (Hk - 1)); + double logRet_crit = 2.492 * std::pow(1 / (Hk - 1.), 0.43) + 0.7 * (std::tanh(14 * Hmi - 9.24) + 1.0); // value at which the amplification starts to grow + + double ax = 0.; + if (Ret <= 0.) + Ret = 1.; + if (std::log10(Ret) < logRet_crit - dgr) + ax = 0.; + else + { + double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2 * dgr); + double Rfac = 0.; + if (rNorm <= 0) + Rfac = 0.; + if (rNorm >= 1) + Rfac = 1.; + else + Rfac = 3 * rNorm * rNorm - 2 * rNorm * rNorm * rNorm; + + // Normal envelope amplification rate + double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi + 0.1 * std::exp(-20 * Hmi); + double arg = 3.87 * Hmi - 2.52; + double dndRet = 0.028 * (Hk - 1) - 0.0345 * std::exp(-(arg * arg)); + ax = (m * dndRet / theta) * Rfac; + + // Set correction for separated profiles + if (Hk > Hk1) + { + double Hnorm = (Hk - Hk1) / (Hk2 - Hk1); + double Hfac = 0.; + if (Hnorm >= 1) + Hfac = 1.; + else + Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm; + double ax1 = ax; + double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5)); + double tnr = std::tanh(1.2 * (std::log10(Ret) - gro)); + double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1), 1.5)) / theta; + if (ax2 < 0.) + ax2 = 0.; + ax = Hfac * ax2 + (1 - Hfac) * ax1; + if (ax < 0.) + ax = 0.; + } + } + return ax; +} \ No newline at end of file diff --git a/vii/src/DFluxes.h b/vii/src/DFluxes.h new file mode 100644 index 0000000..d77449f --- /dev/null +++ b/vii/src/DFluxes.h @@ -0,0 +1,31 @@ +#ifndef DFluxes_H +#define DFluxes_H + +#include "vii.h" +#include "DBoundaryLayer.h" +#include <Eigen/Dense> + +namespace vii +{ + +/** + * @brief Residual and Jacobian computation of the unsteady integral boundary layer equations. + * @author Paul Dechamps + */ +class VII_API Fluxes +{ +public: + double Re; ///< Freestream Reynolds number. + +public: + Fluxes(double _Re) : Re(_Re){}; + ~Fluxes(){}; + Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl); + Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl, Eigen::VectorXd const &sysRes, double eps); + +private: + Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u); + double amplificationRoutine(double Hk, double Ret, double theta) const; +}; +} // namespace vii +#endif // DFluxes_H \ No newline at end of file diff --git a/vii/src/DSolver.cpp b/vii/src/DSolver.cpp new file mode 100644 index 0000000..c08c278 --- /dev/null +++ b/vii/src/DSolver.cpp @@ -0,0 +1,247 @@ +#include "DSolver.h" +#include "DBoundaryLayer.h" +#include "DFluxes.h" +#include <Eigen/Dense> + +using namespace Eigen; +using namespace tbox; +using namespace vii; + +/** + * @brief Construct a new Time Solver:: Time Solver object. + * + * @param _CFL0 Initial CFL number. + * @param _Minf Freestream Mach number. + * @param Re Freestream Reynolds number. + * @param _maxIt Maximum number of iterations. + * @param _tol Convergence tolerance. + * @param _itFrozenJac Number of iterations between which the Jacobian is frozen. + */ +Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt, double _tol, unsigned int _itFrozenJac) +{ + CFL0 = _CFL0; + maxIt = _maxIt; + tol = _tol; + itFrozenJac0 = _itFrozenJac; + Minf = std::max(_Minf, 0.1); + initSln = true; + sysMatrix = new Fluxes(Re); +} + +/** + * @brief Destroy the Time Solver:: Time Solver object. + * + */ +Solver::~Solver() +{ + delete sysMatrix; + std::cout << "~vii::Solver()\n"; +} + +/** + * @brief Impose initial condition at one point. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + */ +void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + if (initSln) + { + bl->u[iPoint * nVar + 0] = bl->u[(iPoint - 1) * nVar + 0]; + bl->u[iPoint * nVar + 1] = bl->u[(iPoint - 1) * nVar + 1]; + bl->u[iPoint * nVar + 2] = bl->u[(iPoint - 1) * nVar + 2]; + bl->u[iPoint * nVar + 3] = bl->blEdge->getVt(iPoint); + bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4]; + + bl->Ret[iPoint] = bl->Ret[iPoint - 1]; + bl->cd[iPoint] = bl->cd[iPoint - 1]; + bl->cf[iPoint] = bl->cf[iPoint - 1]; + bl->Hstar[iPoint] = bl->Hstar[iPoint - 1]; + bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1]; + bl->Hk[iPoint] = bl->Hk[iPoint - 1]; + bl->ctEq[iPoint] = bl->ctEq[iPoint - 1]; + bl->us[iPoint] = bl->us[iPoint - 1]; + bl->delta[iPoint] = bl->delta[iPoint - 1]; + bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1]; + } + if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0) + bl->u[iPoint * nVar + 4] = 0.03; +} + +/** + * @brief Pseudo-time integration. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + * @return int + */ +int Solver::integration(size_t iPoint, BoundaryLayer *bl) +{ + size_t nVar = bl->getnVar(); + + // Save initial condition. + std::vector<double> sln0(nVar, 0.); + for (size_t i = 0; i < nVar; ++i) + sln0[i] = bl->u[iPoint * nVar + i]; + + // Initialize time integration variables. + double dx = bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1); + + // Initial time step. + double CFL = CFL0; + unsigned int itFrozenJac = itFrozenJac0; + double timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint)); + MatrixXd IMatrix(5, 5); + IMatrix = MatrixXd::Identity(5, 5) / timeStep; + + // Initial system residuals. + VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl); + double normSysRes0 = sysRes.norm(); + double normSysRes = normSysRes0; + + MatrixXd JacMatrix = MatrixXd::Zero(5, 5); + VectorXd slnIncr = VectorXd::Zero(5); + + unsigned int innerIters = 0; // Inner (non-linear) iterations + + while (innerIters < maxIt) + { + // Jacobian computation. + if (innerIters % itFrozenJac == 0) + JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8); + + // Linear system. + slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-sysRes); + + // Solution increment. + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] += slnIncr(k); + bl->u[iPoint * nVar + 0] = std::max(bl->u[iPoint * nVar + 0], 1e-6); + bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1], 1.0005); + + sysRes = sysMatrix->computeResiduals(iPoint, bl); + normSysRes = (sysRes + slnIncr / timeStep).norm(); + + if (std::isnan(normSysRes)) + { + resetSolution(iPoint, bl, sln0, true); + sysRes = sysMatrix->computeResiduals(iPoint, bl); + CFL = 1.0; + timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint)); + IMatrix = MatrixXd::Identity(5, 5) / timeStep; + normSysRes = (sysRes + slnIncr / timeStep).norm(); + itFrozenJac = 1; + } + + if (normSysRes <= tol * normSysRes0) + return 0; // Successfull exit. + + // CFL adaptation. + CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1); + timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint)); + IMatrix = MatrixXd::Identity(5, 5) / timeStep; + + innerIters++; + } + + if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3) + { + resetSolution(iPoint, bl, sln0); + return -1; + } + return innerIters; +} + +/** + * @brief Set time step. + * + * @param CFL CFL number. + * @param dx Cell size. + * @param invVel Inviscid velocity. + * @return double + */ +double Solver::setTimeStep(double CFL, double dx, double invVel) const +{ + // Speed of sound. + double gradU2 = (invVel * invVel); + double soundSpeed = computeSoundSpeed(gradU2); + + // Velocity of the fastest travelling wave. + double advVel = soundSpeed + invVel; + + // Time step computation. + return CFL * dx / advVel; +} + +/** + * @brief Compute the speed of sound in a cell. + * + * @param gradU2 Inviscid velocity squared in the cell. + * @return double + */ +double Solver::computeSoundSpeed(double gradU2) const +{ + return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2)); +} + +/** + * @brief Resets the solution of the point wrt its initial condition (sln0) or the previous point. + * + * @param iPoint Marker id. + * @param bl Boundary layer region. + * @param sln0 Initial solution at the point. + * @param usePrevPoint if 1, use the solution at previous point instead of sln0. + */ +void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> sln0, bool usePrevPoint) +{ + size_t nVar = bl->getnVar(); + + // Reset solution. + if (usePrevPoint) + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k]; + else + for (size_t k = 0; k < nVar; ++k) + bl->u[iPoint * nVar + k] = sln0[k]; + + // Reset closures. + bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] * bl->u[iPoint * nVar + 0] * sysMatrix->Re, 1.0); // Reynolds number based on theta. + bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1]; // Displacement thickness. + + if (bl->regime[iPoint] == 0) + { + std::vector<double> lamParam(8, 0.); + bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name); + bl->Hstar[iPoint] = lamParam[0]; + bl->Hstar2[iPoint] = lamParam[1]; + bl->Hk[iPoint] = lamParam[2]; + bl->cf[iPoint] = lamParam[3]; + bl->cd[iPoint] = lamParam[4]; + bl->delta[iPoint] = lamParam[5]; + bl->ctEq[iPoint] = lamParam[6]; + bl->us[iPoint] = lamParam[7]; + } + else if (bl->regime[iPoint] == 1) + { + std::vector<double> turbParam(8, 0.); + bl->closSolver->turbulentClosures(turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name); + bl->Hstar[iPoint] = turbParam[0]; + bl->Hstar2[iPoint] = turbParam[1]; + bl->Hk[iPoint] = turbParam[2]; + bl->cf[iPoint] = turbParam[3]; + bl->cd[iPoint] = turbParam[4]; + bl->delta[iPoint] = turbParam[5]; + bl->ctEq[iPoint] = turbParam[6]; + bl->us[iPoint] = turbParam[7]; + } +} + +void Solver::setCFL0(double _CFL0) +{ + if (_CFL0 > 0) + CFL0 = _CFL0; + else + throw std::runtime_error("Negative CFL numbers are not allowed.\n"); +} \ No newline at end of file diff --git a/vii/src/DSolver.h b/vii/src/DSolver.h new file mode 100644 index 0000000..f3abc6d --- /dev/null +++ b/vii/src/DSolver.h @@ -0,0 +1,51 @@ +#ifndef DSolver_H +#define DSolver_H + +#include "vii.h" +#include "DFluxes.h" + +namespace vii +{ + +/** + * @brief Pseudo-time integration. + * @author Paul Dechamps + */ +class VII_API Solver +{ + +private: + double CFL0; ///< Initial CFL number. + unsigned int maxIt; ///< Maximum number of iterations. + double tol; ///< Convergence tolerance. + unsigned int itFrozenJac0; ///< Number of iterations between which the Jacobian is frozen. + bool initSln; ///< Flag to initialize the solution at the point. + double const gamma = 1.4; ///< Air heat capacity ratio. + double Minf; ///< Freestream Mach number. + + Fluxes *sysMatrix; + +public: + Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100, double _tol = 1e-8, unsigned int _itFrozenJac = 5); + ~Solver(); + void initialCondition(size_t iPoint, BoundaryLayer *bl); + int integration(size_t iPoint, BoundaryLayer *bl); + + // Getters. + double getCFL0() const { return CFL0; }; + unsigned int getmaxIt() const { return maxIt; }; + unsigned int getitFrozenJac() const { return itFrozenJac0; }; + + // Setters. + void setCFL0(double _CFL0); + void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; }; + void setitFrozenJac(unsigned int _itFrozenJac) { itFrozenJac0 = _itFrozenJac; }; + void setinitSln(bool _initSln) { initSln = _initSln; }; + +private: + double setTimeStep(double CFL, double dx, double invVel) const; + double computeSoundSpeed(double gradU2) const; + void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0, bool usePrevPoint = false); +}; +} // namespace vii +#endif // DSolver_H \ No newline at end of file diff --git a/vii/src/vii.h b/vii/src/vii.h new file mode 100644 index 0000000..300376a --- /dev/null +++ b/vii/src/vii.h @@ -0,0 +1,53 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +// Global header of the "vii" module. + +#ifndef vii_H +#define vii_H + +#if defined(WIN32) +#ifdef vii_EXPORTS +#define VII_API __declspec(dllexport) +#else +#define VII_API __declspec(dllimport) +#endif +#else +#define VII_API +#endif + +#include "tbox.h" + +/** + * @brief Namespace for vii module + */ +namespace vii +{ + +// Solvers +class Driver; +class Solver; + +// Data structures +class BoundaryLayer; +class Discretization; +class Edge; + +// Other +class Closures; +} // namespace vii + +#endif // VII_H \ No newline at end of file diff --git a/vii/tests/bli2D.py b/vii/tests/bli2D.py new file mode 100644 index 0000000..34a8247 --- /dev/null +++ b/vii/tests/bli2D.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +# Copyright 2022 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. + + +# @author Paul Dechamps +# @date 2022 +# Test the vii implementation. The test case is a compressible attached transitional flow. +# Tested functionalities; +# - Time integration. +# - Two time-marching techniques (agressive and safe). +# - Transition routines. +# - Compressible flow routines. +# - Laminar and turbulent flow. +# +# Untested functionalities. +# - Separation routines. +# - Multiple failure case at one iteration. +# - Response to unconverged inviscid solver. + +# Imports. +import dart.utils as invUtils +import dart.default as invDefault + +import vii.coupler as viiCoupler +import vii.utils as viscUtils + +import tbox.utils as tboxU +import fwk +from fwk.testing import * +from fwk.coloring import ccolors + +import math + +def main(): + # Timer. + tms = fwk.Timers() + tms['total'].start() + + # Flow variables. + Re = 1e7 + alpha = 2.*math.pi/180 + M_inf = 0.2 + # Numerical parameters. + CFL0 = 1 + + # Mesh the geometry. + print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) + tms['msh'].start() + + pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01} + msh, gmshWritter = invDefault.mesh(2, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream']) + tms['msh'].stop() + + # Set the problem and add medium, initial/boundary conditions. + tms['pre'].start() + pbl, _, _, bnd, blw = invDefault.problem(msh, 2, alpha, 0., M_inf, 1.0, 1.0, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) + tms['pre'].stop() + + # Solve coupled problem. + print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) + tms['solver'].start() + + vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1) + iSolverAPI = viscUtils.initDart(pbl, blw, vSolver) + Coupler = viiCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=20, _couplTol=1e-4, _iterPrint=1, _resetInv=False) + Coupler.run() + + # Save results. + iSolverAPI.save(gmshWritter) + Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolverAPI.solver.Cp) + tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '') + + print('') + # Display results. + print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) + print(' Re M alpha Cl Cd Cdp Cdf Cm') + print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm())) + + # Write results to file. + vSolution = viscUtils.getSolution(vSolver) + viscUtils.writeFile(vSolution, vSolver.getRe()) + + # Display timers. + tms['total'].stop() + print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) + print('CPU statistics') + print(tms) + print('SOLVERS statistics') + print(Coupler.tms) + + # Check results. + # Test for n0012 (le=0.01, te=0.01), alpha=2, M=.2, Re=1e7. + print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) + tests = CTests() + tests.add(CTest('Cl', iSolverAPI.getCl(), 0.234, 5e-3)) # XFOIL 0.2325 + tests.add(CTest('Cd', vSolution['Cdt'], 0.0055, 1e-3, forceabs=True)) # XFOIL 0.00531 + tests.add(CTest('Cdp', vSolution['Cdp'], 0.0012, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447 + tests.add(CTest('xtr_top', vSolution['xtrT'], 0.21, 0.03, forceabs=True)) # XFOIL 0.1877 + tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.50, 0.03, forceabs=True)) # XFOIL 0.4998 + tests.run() + + # Plot results + tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolverAPI.getCl(), vSolver.Cdt, iSolverAPI.getCm(), 4), True) + tboxU.plot(vSolution['x/c'], vSolution['Cf'], '$x/c$', '$c_f$', 'Cdt = {0:.{3}f}, Cdf = {1:.{3}f}, Cdp = {2:.{3}f}'.format(vSolver.Cdt, vSolver.Cdf, vSolver.Cdp, 4)) + # eof. + print('') + +if __name__ == "__main__": + main() diff --git a/vii/utils.py b/vii/utils.py new file mode 100644 index 0000000..b855767 --- /dev/null +++ b/vii/utils.py @@ -0,0 +1,143 @@ +import vii +from fwk.wutils import parseargs +from fwk.coloring import ccolors +import numpy as np + +def initBL(Re, Minf, CFL0, nSections, span=0, verb=None): + """ Initialize boundary layer solver. + + Params + ------ + - Re: Flow Reynolds number. + - Minf: Freestream Mach number. + - CFL0: Initial CFL number for time integration. + - nSections: Number of sections in the domain. + - span: Wing span (not used for 2D computations. + - verb: Verbosity level of the solver. + + Return + ------ + - solver: vii::Driver class. + """ + if Re<=0.: + print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET) + raise RuntimeError("Invalid parameter") + if Minf<0.: + print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Mach number.", Minf, ccolors.ANSI_RESET) + raise RuntimeError("Invalid parameter") + elif Minf>=1.: + print(ccolors.ANSI_YELLOW, "dart::vUtils Warning : (Super)sonic freestream Mach number.", Minf, ccolors.ANSI_RESET) + if nSections < 0: + print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET) + raise RuntimeError("Invalid parameter") + if verb is None: + args = parseargs() + _verbose = args.verb + else: + if not(0<=verb<=3): + print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", _verbose, ccolors.ANSI_RESET) + raise RuntimeError("Invalid parameter") + else: + _verbose = verb + solver = vii.Driver(Re, Minf, CFL0, nSections, span, verbose=_verbose) + return solver + +def initDart(pbl, blw, vSolver, iters = 25): + """Initialize Newton solver and create the interface object + between dart and the viscous solver. + + Parameters + ---------- + pbl (dart.Problem object): Problem formulation. + blw (np.ndarray(dart.Blowing, dart.Blowing)): Blowing boundaries for (airfoil, wake). + vSolver (vii.Driver object): Viscous solver. + + Return + ------ + solverAPI (DartInterface object): Interface between Dart and the viscous solver. + """ + import dart + from vii.interfaces.dart.DartInterface import DartInterface as dInterface + from tbox.solvers import LinearSolver + + args = parseargs() + dartSolver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb, mIt=iters) + solverAPI = dInterface(dartSolver, vSolver, blw) + return solverAPI + +def mesh(file, pars): + """Initialize mesh and mesh writer + + Parameters + ---------- + file : str + file contaning mesh (w/o extention) + pars : dict + parameters for mesh + """ + import tbox.gmsh as tmsh + # Load the mesh. + msh = tmsh.MeshLoader(file,__file__).execute(**pars) + return msh + + +def getSolution(vSolver, iSec=0): + """Extract viscous solution. + Args + ---- + - vSolver: vii::Driver class. Drive of the viscous calculations + - iSec (int): Marker of the section (default: 0) + """ + if iSec<0: + raise RuntimeError("dart::viscU Invalid section number", iSec) + + solverOutput = vSolver.getSolution(iSec) + + sln = {'theta' : solverOutput[0], + 'H' : solverOutput[1], + 'N' : solverOutput[2], + 'ue' : solverOutput[3], + 'Ct' : solverOutput[4], + 'deltaStar': solverOutput[5], + 'Ret' : solverOutput[6], + 'Cd' : solverOutput[7], + 'Cf' : np.asarray(solverOutput[8])*np.asarray(solverOutput[3])**2, + 'Hstar' : solverOutput[9], + 'Hstar2' : solverOutput[10], + 'Hk' : solverOutput[11], + 'Cteq' : solverOutput[12], + 'Us' : solverOutput[13], + 'delta' : solverOutput[14], + 'x/c' : solverOutput[15], + 'xtrT' : vSolver.getxtr(iSec, 0), + 'xtrB' : vSolver.getxtr(iSec, 1), + 'Cdt' : vSolver.Cdt, + 'Cdf' : vSolver.Cdf, + 'Cdp' : vSolver.Cdp} + return sln + +def writeFile(wData, Re, toW=['deltaStar', 'H', 'Cf', 'Ct', 'ue', 'delta']): + """Write the results in airfoil_viscous.dat + """ + # Write + print('Writing file: airfoil_viscous.dat...', end = '') + f = open('airfoil_viscous.dat', 'w+') + + f.write('$Aerodynamic coefficients\n') + f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n') + f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt'], wData['Cdp'],wData['Cdf'], wData['xtrT'], wData['xtrB'])) + f.write('$Boundary layer variables\n') + f.write(' x/c') + + for s in toW: + f.write(' {0:>15s}'.format(s)) + f.write('\n') + + for i in range(len(wData['x/c'])): + f.write('{0:15.6f}'.format(wData['x/c'][i])) + for s in toW: + f.write(' {0:15.6f}'.format(wData[s][i])) + f.write('\n') + + f.close() + print('done.') \ No newline at end of file -- GitLab