Skip to content
Snippets Groups Projects
Commit 641fd0ca authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon: Committed by Adrien Crovato
Browse files

Added vii to dart

parent b3dec627
No related branches found
No related tags found
No related merge requests found
Pipeline #8613 passed
Stage: build
    Stage: test
      Expand all Hide whitespace changes
      Inline Side-by-side
      Showing
      with 1660 additions and 1728 deletions
      CMakeLists.txt
      + 1
      − 0
      View file @ 641fd0ca
      ......@@ -114,6 +114,7 @@ ENDIF()
      # -- Sub directories
      ADD_SUBDIRECTORY( ext )
      ADD_SUBDIRECTORY( dart )
      ADD_SUBDIRECTORY( vii )
      # -- FINAL
      MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}")
      ......
      dart/CMakeLists.txt
      + 0
      − 1
      View file @ 641fd0ca
      ......@@ -24,5 +24,4 @@ INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
      ${CMAKE_CURRENT_LIST_DIR}/utils.py
      DESTINATION dart)
      INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
      ${CMAKE_CURRENT_LIST_DIR}/viscous
      DESTINATION dart)
      dart/api/core.py
      + 266
      − 236
      View file @ 641fd0ca
      #!/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
      }
      dart/api/internal/polar.py
      + 137
      − 131
      View file @ 641fd0ca
      #!/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', '')
      dart/api/internal/trim.py
      + 137
      − 131
      View file @ 641fd0ca
      #!/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)
      dart/api/mphys.py
      + 573
      − 573
      View file @ 641fd0ca
      This diff is collapsed.
      dart/viscous/coupler.py deleted 100644 → 0
      + 0
      − 82
      View file @ b3dec627
      #!/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')
      dart/viscous/newton.py deleted 100644 → 0
      + 0
      − 75
      View file @ b3dec627
      #!/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)
      dart/viscous/solver.py deleted 100644 → 0
      + 0
      − 489
      View file @ b3dec627
      This diff is collapsed.
      vii/CMakeLists.txt 0 → 100644
      + 27
      − 0
      View file @ 641fd0ca
      # 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)
      vii/__init__.py 0 → 100644
      + 21
      − 0
      View file @ 641fd0ca
      # -*- 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 *
      vii/_src/CMakeLists.txt 0 → 100644
      + 35
      − 0
      View file @ 641fd0ca
      # 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})
      vii/_src/viiw.h 0 → 100644
      + 17
      − 0
      View file @ 641fd0ca
      /*
      * 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
      vii/_src/viiw.i 0 → 100644
      + 57
      − 0
      View file @ 641fd0ca
      /*
      * 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
      dart/viscous/__init__.py → vii/api/__init__.py
      + 0
      − 0
      View file @ 641fd0ca
      File moved
      vii/api/core.py 0 → 100644
      + 110
      − 0
      View file @ 641fd0ca
      #!/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}
      vii/coupler.py 0 → 100644
      + 97
      − 0
      View file @ 641fd0ca
      #!/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
      vii/interfaces/dart/DartInterface.py 0 → 100644
      + 167
      − 0
      View file @ 641fd0ca
      #!/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
      vii/interfaces/dart/__init__.py 0 → 100644
      + 1
      − 0
      View file @ 641fd0ca
      # -*- coding: utf-8 -*-
      \ No newline at end of file
      dart/viscous/airfoil.py → vii/interfaces/viscous/DAirfoil.py 100644 → 100755
      + 14
      − 10
      View file @ 641fd0ca
      ......@@ -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
      ......
      0% or .
      You are about to add 0 people to the discussion. Proceed with caution.
      Finish editing this message first!
      Please register or to comment