Skip to content
Snippets Groups Projects
Verified Commit 5a78e517 authored by Paul Dechamps's avatar Paul Dechamps :speech_balloon:
Browse files

(chore) Removed vii from DARTFLO 

Use BLASTER (gitlab.uliege.be/am-dept/blaster) instead.
parent 6f16bb57
No related branches found
No related tags found
No related merge requests found
Pipeline #49834 passed
Stage: build
    Stage: test
      Hide whitespace changes
      Inline Side-by-side
      Showing
      with 0 additions and 1501 deletions
      CMakeLists.txt
      + 0
      − 1
      View file @ 5a78e517
      ......@@ -114,7 +114,6 @@ ENDIF()
      # -- Sub directories
      ADD_SUBDIRECTORY( ext )
      ADD_SUBDIRECTORY( dart )
      ADD_SUBDIRECTORY( vii )
      # -- FINAL
      MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}")
      ......
      dart/README.md
      + 0
      − 1
      View file @ 5a78e517
      ......@@ -2,6 +2,5 @@
      * [tests](dart/tests): simple tests to be run in battery checking the basic functionnalities of the solver
      * [validation](dart/validation): large-scale tests used to validate the solver results
      * [benchmark](dart/benchmark): large-scale tests used to monitor the performance
      * [viscous](dart/viscous): viscous-inviscid interaction module
      * [api](dart/api): APIs for DART
      * [cases](dart/cases): sample cases meant to be run using the APIs
      vii/CMakeLists.txt deleted 100644 → 0
      + 0
      − 29
      View file @ 6f16bb57
      # 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}/coupler.py
      ${CMAKE_CURRENT_LIST_DIR}/utils.py
      DESTINATION vii)
      INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
      ${CMAKE_CURRENT_LIST_DIR}/interfaces
      DESTINATION vii)
      vii/__init__.py deleted 100644 → 0
      + 0
      − 21
      View file @ 6f16bb57
      # -*- 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 deleted 100644 → 0
      + 0
      − 35
      View file @ 6f16bb57
      # 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 ${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 deleted 100644 → 0
      + 0
      − 17
      View file @ 6f16bb57
      /*
      * 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 deleted 100644 → 0
      + 0
      − 57
      View file @ 6f16bb57
      /*
      * 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
      vii/api/__init__.py deleted 100644 → 0
      + 0
      − 1
      View file @ 6f16bb57
      # -*- coding: utf-8 -*-
      vii/api/core.py deleted 100644 → 0
      + 0
      − 117
      View file @ 6f16bb57
      #!/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.
      - xtrF ([float, float]): Forced transition location [upper, lower].
      - 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 init_dart
      iSolverObjects = init_dart(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
      _xtrF = [-1, -1]
      if 'xtrF' in cfg:
      for i, xtr in enumerate(cfg['xtrF']):
      if not(0 <= xtr <= 1) and xtr != -1:
      raise RuntimeError("Incorrect forced transition location.")
      _xtrF[i] = xtr
      _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, _xtrF[0], _xtrF[1], _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 deleted 100644 → 0
      + 0
      − 97
      View file @ 6f16bb57
      #!/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(sfx='_Inviscid')
      # 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}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(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}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>6s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
      print(' ----------------------------------------------------------')
      if couplIter % self.iterPrint == 0:
      print(' {:>4.0f}| {:>7.4f} {:>7.4f} {:>7.4f} | {:>6.4f} {:>6.4f} | {:>6.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(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}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>8s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
      print(' ----------------------------------------------------------')
      print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
      return aeroCoeffs
      \ No newline at end of file
      vii/interfaces/dart/DartInterface.py deleted 100644 → 0
      + 0
      − 170
      View file @ 6f16bb57
      #!/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, sfx=''):
      """ Write Cp distribution around the airfoil on file.
      """
      Cp = floU.extract(self.group[0].boundary.tag.elems, self.solver.Cp)
      tboxU.write(Cp, 'Cp'+sfx+'.dat', '%1.4e', ',', 'x, y, 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 getCd(self):
      return self.solver.Cd
      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 deleted 100644 → 0
      + 0
      − 1
      View file @ 6f16bb57
      # -*- coding: utf-8 -*-
      \ No newline at end of file
      vii/interfaces/viscous/DAirfoil.py deleted 100755 → 0
      + 0
      − 135
      View file @ 6f16bb57
      #!/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.
      ## Airfoil around which the boundary layer is computed
      # Amaury Bilocq
      from .DBoundary import Boundary
      import numpy as np
      class Airfoil(Boundary):
      def __init__(self, _boundary):
      Boundary.__init__(self, _boundary)
      self.name = 'airfoil' # type of boundary
      self.T0 = 0 # initial condition for the momentum thickness
      self.H0 = 0
      self.n0 = 0
      self.ct0 = 0
      self.TEnd = [0,0]
      self.HEnd = [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))
      else:
      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
      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
      data = np.zeros((self.boundary.nodes.size(), 10))
      i = 0
      for n in self.boundary.nodes:
      data[i,0] = n.no
      data[i,1] = n.pos[0]
      data[i,2] = n.pos[1]
      data[i,3] = n.pos[2]
      data[i,4] = self.v[i,0]
      data[i,5] = self.v[i,1]
      data[i,6] = self.Me[i]
      data[i,7] = self.rhoe[i]
      data[i,8] = self.deltaStar[i]
      data[i,9] = self.xx[i]
      i += 1
      # Table containing the element and its nodes
      eData = np.zeros((self.nE,3), dtype=int)
      for i in range(0, self.nE):
      eData[i,0] = self.boundary.tag.elems[i].no
      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
      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
      idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
      idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
      ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
      ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
      if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
      upperGlobTE = int(data[idxTE[0],0])
      lowerGlobTE = int(data[idxTE[1],0])
      else:
      upperGlobTE = int(data[idxTE[1],0])
      lowerGlobTE = int(data[idxTE[0],0])
      # Connectivity
      connectListElems[0] = np.where(eData[:,1] == globStag)[0]
      N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
      connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
      i = 1
      upperTE = 0
      lowerTE = 0
      # Sort the suction part
      while upperTE == 0:
      N1[i] = eData[connectListElems[i-1],2] # Second node of the element
      connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
      connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
      if eData[connectListElems[i],2] == upperGlobTE:
      upperTE = 1
      i += 1
      # Sort the pressure side
      connectListElems[i] = np.where(eData[:,2] == globStag)[0]
      connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
      N1[i] = eData[connectListElems[i],2]
      while lowerTE == 0:
      N1[i+1] = eData[connectListElems[i],1] # First node of the element
      connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
      connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
      if eData[connectListElems[i+1],1] == lowerGlobTE:
      lowerTE = 1
      i += 1
      connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
      data[:,0] = data[connectListNodes,0]
      data[:,1] = data[connectListNodes,1]
      data[:,2] = data[connectListNodes,2]
      data[:,3] = data[connectListNodes,3]
      data[:,4] = data[connectListNodes,4]
      data[:,5] = data[connectListNodes,5]
      data[:,6] = data[connectListNodes,6]
      data[:,7] = data[connectListNodes,7]
      data[:,8] = data[connectListNodes,8]
      data[:,9] = data[connectListNodes,9]
      # Separated upper and lower part
      data = np.delete(data,0,1)
      uData = data[0:np.argmax(data[:,0])+1]
      lData = data[np.argmax(data[:,0])+1:None]
      lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
      return connectListNodes, connectListElems,[uData, lData]
      vii/interfaces/viscous/DBoundary.py deleted 100755 → 0
      + 0
      − 36
      View file @ 6f16bb57
      #!/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.
      ## Base class representing a physical boundary
      # Amaury Bilocq
      import numpy as np
      class Boundary:
      def __init__(self, _boundary):
      self.boundary = _boundary # boundary of interest
      self.nN = len(self.boundary.nodes) # number of nodes
      self.nE = len(self.boundary.tag.elems) # number of elements
      self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
      self.u = np.zeros(self.nE) # blowing velocity
      self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
      self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
      self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
      self.connectListElems = np.zeros(self.nE) # connectivity for elements
      self.deltaStar = np.zeros(self.nN) # displacement thickness
      self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law
      vii/interfaces/viscous/DGroupBuilder.py deleted 100755 → 0
      + 0
      − 169
      View file @ 6f16bb57
      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
      vii/interfaces/viscous/DWake.py deleted 100755 → 0
      + 0
      − 77
      View file @ 6f16bb57
      #!/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.
      ## Wake behind airfoil (around which the boundary layer is computed)
      # Amaury Bilocq & Paul Dechamps
      from .DBoundary import Boundary
      import numpy as np
      class Wake(Boundary):
      def __init__(self, _boundary):
      Boundary.__init__(self, _boundary)
      self.name = 'wake' # type of boundary
      self.T0 = 0 # initial condition for the momentum thickness
      self.H0 = 0
      self.n0 = 9 # wake is always turbulent
      self.ct0 = 0
      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
      data = np.zeros((self.boundary.nodes.size(), 10))
      i = 0
      for n in self.boundary.nodes:
      data[i,0] = n.no
      data[i,1] = n.pos[0]
      data[i,2] = n.pos[1]
      data[i,3] = n.pos[2]
      data[i,4] = self.v[i,0]
      data[i,5] = self.v[i,1]
      data[i,6] = self.Me[i]
      data[i,7] = self.rhoe[i]
      data[i,8] = self.deltaStar[i]
      data[i,9] = self.xx[i]
      i += 1
      # Table containing the element and its nodes
      eData = np.zeros((self.nE,3), dtype=int)
      for i in range(0, self.nE):
      eData[i,0] = self.boundary.tag.elems[i].no
      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] = 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]
      data[:,1] = data[connectListNodes,1]
      data[:,2] = data[connectListNodes,2]
      data[:,3] = data[connectListNodes,3]
      data[:,4] = data[connectListNodes,4]
      data[:,5] = data[connectListNodes,5]
      data[:,6] = data[connectListNodes,6]
      data[:,7] = data[connectListNodes,7]
      data[:,8] = data[connectListNodes,8]
      data[:,9] = data[connectListNodes,9]
      # Separated upper and lower part
      data = np.delete(data,0,1)
      return connectListNodes, connectListElems, data
      vii/src/CMakeLists.txt deleted 100644 → 0
      + 0
      − 27
      View file @ 6f16bb57
      # 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)")
      vii/src/DBoundaryLayer.cpp deleted 100644 → 0
      + 0
      − 185
      View file @ 6f16bb57
      #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;
      xoctr = 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
      vii/src/DBoundaryLayer.h deleted 100644 → 0
      + 0
      − 69
      View file @ 6f16bb57
      #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.
      double xoctr; ///< Transition location in % of the chord.
      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
      vii/src/DClosures.cpp deleted 100644 → 0
      + 0
      − 256
      View file @ 6f16bb57
      #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
      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