From 5a78e5175e1bc450ae28f8c80a284d3e90fa438f Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Fri, 22 Nov 2024 10:05:52 +0100
Subject: [PATCH] (chore) Removed vii from DARTFLO
Use BLASTER (gitlab.uliege.be/am-dept/blaster) instead.
---
CMakeLists.txt | 1 -
dart/README.md | 1 -
vii/CMakeLists.txt | 29 -
vii/__init__.py | 21 -
vii/_src/CMakeLists.txt | 35 --
vii/_src/viiw.h | 17 -
vii/_src/viiw.i | 57 --
vii/api/__init__.py | 1 -
vii/api/core.py | 117 ----
vii/coupler.py | 97 ---
vii/interfaces/dart/DartInterface.py | 170 -----
vii/interfaces/dart/__init__.py | 1 -
vii/interfaces/viscous/DAirfoil.py | 135 ----
vii/interfaces/viscous/DBoundary.py | 36 --
vii/interfaces/viscous/DGroupBuilder.py | 169 -----
vii/interfaces/viscous/DWake.py | 77 ---
vii/src/CMakeLists.txt | 27 -
vii/src/DBoundaryLayer.cpp | 185 ------
vii/src/DBoundaryLayer.h | 69 --
vii/src/DClosures.cpp | 256 --------
vii/src/DClosures.h | 30 -
vii/src/DDiscretization.cpp | 76 ---
vii/src/DDiscretization.h | 54 --
vii/src/DDriver.cpp | 799 ------------------------
vii/src/DDriver.h | 76 ---
vii/src/DEdge.cpp | 59 --
vii/src/DEdge.h | 44 --
vii/src/DFluxes.cpp | 242 -------
vii/src/DFluxes.h | 31 -
vii/src/DSolver.cpp | 247 --------
vii/src/DSolver.h | 51 --
vii/src/vii.h | 53 --
vii/tests/bli2D.py | 126 ----
vii/utils.py | 149 -----
34 files changed, 3538 deletions(-)
delete mode 100644 vii/CMakeLists.txt
delete mode 100644 vii/__init__.py
delete mode 100644 vii/_src/CMakeLists.txt
delete mode 100644 vii/_src/viiw.h
delete mode 100644 vii/_src/viiw.i
delete mode 100644 vii/api/__init__.py
delete mode 100644 vii/api/core.py
delete mode 100644 vii/coupler.py
delete mode 100644 vii/interfaces/dart/DartInterface.py
delete mode 100644 vii/interfaces/dart/__init__.py
delete mode 100755 vii/interfaces/viscous/DAirfoil.py
delete mode 100755 vii/interfaces/viscous/DBoundary.py
delete mode 100755 vii/interfaces/viscous/DGroupBuilder.py
delete mode 100755 vii/interfaces/viscous/DWake.py
delete mode 100644 vii/src/CMakeLists.txt
delete mode 100644 vii/src/DBoundaryLayer.cpp
delete mode 100644 vii/src/DBoundaryLayer.h
delete mode 100644 vii/src/DClosures.cpp
delete mode 100644 vii/src/DClosures.h
delete mode 100644 vii/src/DDiscretization.cpp
delete mode 100644 vii/src/DDiscretization.h
delete mode 100644 vii/src/DDriver.cpp
delete mode 100644 vii/src/DDriver.h
delete mode 100644 vii/src/DEdge.cpp
delete mode 100644 vii/src/DEdge.h
delete mode 100644 vii/src/DFluxes.cpp
delete mode 100644 vii/src/DFluxes.h
delete mode 100644 vii/src/DSolver.cpp
delete mode 100644 vii/src/DSolver.h
delete mode 100644 vii/src/vii.h
delete mode 100644 vii/tests/bli2D.py
delete mode 100644 vii/utils.py
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f7e8db5..ec130f6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -114,7 +114,6 @@ ENDIF()
# -- Sub directories
ADD_SUBDIRECTORY( ext )
ADD_SUBDIRECTORY( dart )
-ADD_SUBDIRECTORY( vii )
# -- FINAL
MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}")
diff --git a/dart/README.md b/dart/README.md
index 71a9c79..b1cf6b8 100644
--- a/dart/README.md
+++ b/dart/README.md
@@ -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
diff --git a/vii/CMakeLists.txt b/vii/CMakeLists.txt
deleted file mode 100644
index fba4e43..0000000
--- a/vii/CMakeLists.txt
+++ /dev/null
@@ -1,29 +0,0 @@
-# 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)
diff --git a/vii/__init__.py b/vii/__init__.py
deleted file mode 100644
index 43be8da..0000000
--- a/vii/__init__.py
+++ /dev/null
@@ -1,21 +0,0 @@
-# -*- coding: utf-8 -*-
-# dart MODULE initialization file
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-import fwk
-import tbox
-from viiw import *
diff --git a/vii/_src/CMakeLists.txt b/vii/_src/CMakeLists.txt
deleted file mode 100644
index e00df4b..0000000
--- a/vii/_src/CMakeLists.txt
+++ /dev/null
@@ -1,35 +0,0 @@
-# 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})
diff --git a/vii/_src/viiw.h b/vii/_src/viiw.h
deleted file mode 100644
index ec45463..0000000
--- a/vii/_src/viiw.h
+++ /dev/null
@@ -1,17 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-#include "DDriver.h"
\ No newline at end of file
diff --git a/vii/_src/viiw.i b/vii/_src/viiw.i
deleted file mode 100644
index cc8f9da..0000000
--- a/vii/_src/viiw.i
+++ /dev/null
@@ -1,57 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-// SWIG input file of the 'dart' module
-
-%feature("autodoc","1");
-
-%module(docstring=
-"'viiw' module: Viscous inviscid interaction for dart 2021-2022
-(c) ULiege - A&M",
-directors="0",
-threads="1"
-) viiw
-%{
-
-#include "vii.h"
-
-#include "fwkw.h"
-#include "tboxw.h"
-#include "viiw.h"
-
-%}
-
-
-%include "fwkw.swg"
-
-// ----------- MODULES UTILISES ------------
-%import "tboxw.i"
-
-// ----------- VII CLASSES ----------------
-%include "vii.h"
-
-%shared_ptr(vii::Driver);
-
-%immutable vii::Driver::Re; // read-only variable (no setter)
-%immutable vii::Driver::Cdt;
-%immutable vii::Driver::Cdt_sec;
-%immutable vii::Driver::Cdf;
-%immutable vii::Driver::Cdf_sec;
-%immutable vii::Driver::Cdp;
-%immutable vii::Driver::Cdp_sec;
-%immutable vii::Driver::CFL0;
-%immutable vii::Driver::verbose;
-%include "DDriver.h"
\ No newline at end of file
diff --git a/vii/api/__init__.py b/vii/api/__init__.py
deleted file mode 100644
index 40a96af..0000000
--- a/vii/api/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-# -*- coding: utf-8 -*-
diff --git a/vii/api/core.py b/vii/api/core.py
deleted file mode 100644
index ff23fe5..0000000
--- a/vii/api/core.py
+++ /dev/null
@@ -1,117 +0,0 @@
-#!/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}
diff --git a/vii/coupler.py b/vii/coupler.py
deleted file mode 100644
index c60a22a..0000000
--- a/vii/coupler.py
+++ /dev/null
@@ -1,97 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# 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
diff --git a/vii/interfaces/dart/DartInterface.py b/vii/interfaces/dart/DartInterface.py
deleted file mode 100644
index c347fdb..0000000
--- a/vii/interfaces/dart/DartInterface.py
+++ /dev/null
@@ -1,170 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# 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
diff --git a/vii/interfaces/dart/__init__.py b/vii/interfaces/dart/__init__.py
deleted file mode 100644
index 7c68785..0000000
--- a/vii/interfaces/dart/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-# -*- coding: utf-8 -*-
\ No newline at end of file
diff --git a/vii/interfaces/viscous/DAirfoil.py b/vii/interfaces/viscous/DAirfoil.py
deleted file mode 100755
index 79d9494..0000000
--- a/vii/interfaces/viscous/DAirfoil.py
+++ /dev/null
@@ -1,135 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## 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]
diff --git a/vii/interfaces/viscous/DBoundary.py b/vii/interfaces/viscous/DBoundary.py
deleted file mode 100755
index f3ff33a..0000000
--- a/vii/interfaces/viscous/DBoundary.py
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## 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
diff --git a/vii/interfaces/viscous/DGroupBuilder.py b/vii/interfaces/viscous/DGroupBuilder.py
deleted file mode 100755
index d4fcbca..0000000
--- a/vii/interfaces/viscous/DGroupBuilder.py
+++ /dev/null
@@ -1,169 +0,0 @@
-import numpy as np
-
-def mergeVectors(dataIn, mapMesh, nFactor):
- """Merges 3 groups upper/lower sides and wake.
- """
- for k in range(len(dataIn)):
- if len(dataIn[k]) == 2: # Airfoil case.
- xx_up, dv_up, vtExt_up, _, alpha_up = getBLcoordinates(dataIn[k][0])
- xx_lw, dv_lw, vtExt_lw, _, alpha_lw = getBLcoordinates(dataIn[k][1])
- else: # Wake case
- xx_wk, dv_wk, vtExt_wk, _, alpha_wk = getBLcoordinates(dataIn[k])
-
- if mapMesh:
- # Save initial mesh.
- cMesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][1:, 0], dataIn[1][:, 0]))
- cMeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][1:, 0])]
- cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
- cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
- cxxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][1:, 8], dataIn[1][:, 8]))
-
- # Create fine mesh.
- fMeshUp = createFineMesh(dataIn[0][0][:, 0], nFactor)
- fMeshLw = createFineMesh(dataIn[0][1][:, 0], nFactor)
- fMeshWk = createFineMesh(dataIn[1][:, 0], nFactor)
- fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
- fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
-
- fxx_up = interpolateToFineMesh(xx_up, fMeshUp, nFactor)
- fxx_lw = interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
- fxx_wk = interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
- fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
-
- fvtExt_up = interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
- fvtExt_lw = interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
- fvtExt_wk = interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
- fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
-
- fMe_up = interpolateToFineMesh(dataIn[0][0][:, 5], fMeshUp, nFactor)
- fMe_lw = interpolateToFineMesh(dataIn[0][1][:, 5], fMeshLw, nFactor)
- fMe_wk = interpolateToFineMesh(dataIn[1][:, 5], fMeshWk, nFactor)
- fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
-
- frho_up = interpolateToFineMesh(dataIn[0][0][:, 6], fMeshUp, nFactor)
- frho_lw = interpolateToFineMesh(dataIn[0][1][:, 6], fMeshLw, nFactor)
- frho_wk = interpolateToFineMesh(dataIn[1][:, 6], fMeshWk, nFactor)
- frho = np.concatenate((frho_up, frho_lw[1:], frho_wk))
-
- fdStarExt_up = interpolateToFineMesh(dataIn[0][0][:, 7], fMeshUp, nFactor)
- fdStarExt_lw = interpolateToFineMesh(dataIn[0][1][:, 7], fMeshLw, nFactor)
- fdStarExt_wk = interpolateToFineMesh(dataIn[1][:, 7], fMeshWk, nFactor)
-
- fxxExt_up = interpolateToFineMesh(dataIn[0][0][:, 8], fMeshUp, nFactor)
- fxxExt_lw = interpolateToFineMesh(dataIn[0][1][:, 8], fMeshLw, nFactor)
- fxxExt_wk = interpolateToFineMesh(dataIn[1][:, 8], fMeshWk, nFactor)
-
- fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
- fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
- fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
-
- # Create dictionnaries for the fine and coarse meshes.
- fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
- 'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
-
- cMe = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][1:, 5], dataIn[1][:, 5]))
- cRho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][1:, 6], dataIn[1][:, 6]))
- cdStarExt = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][1:, 7], dataIn[1][:, 7]))
- cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
- 'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
-
- dataUpper = np.zeros((len(fMeshUp), 0))
- dataLower = np.zeros((len(fMeshLw), 0))
- dataWake = np.zeros((len(fMeshWk), 0))
- for iParam in range(len(dataIn[0][0][0, :])):
- interpParamUp = interpolateToFineMesh(dataIn[0][0][:, iParam], fMeshUp, nFactor)
- dataUpper = np.column_stack((dataUpper, interpParamUp))
- interpParamLw = interpolateToFineMesh(dataIn[0][1][:, iParam], fMeshLw, nFactor)
- dataLower = np.column_stack((dataLower, interpParamLw))
- interpParamWk = interpolateToFineMesh(dataIn[1][:, iParam], fMeshWk, nFactor)
- dataWake = np.column_stack((dataWake, interpParamWk))
- # Remove stagnation point doublon.
- dataLower = np.delete(dataLower, 0, axis=0)
- else:
- Mesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][:, 0], dataIn[1][:, 0]))
- Meshy = np.concatenate((dataIn[0][0][:, 1], dataIn[0][1][:, 1], dataIn[1][:, 1]))
- Meshz = np.concatenate((dataIn[0][0][:, 2], dataIn[0][1][:, 2], dataIn[1][:, 2]))
- MeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]) + len(dataIn[1][:, 0])]
-
- # Concatenate all parameters (upper side, lower side, wake)
- xx = np.concatenate((xx_up, xx_lw, xx_wk))
- vtExt = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk))
- vx = np.concatenate((dataIn[0][0][:, 3], dataIn[0][1][:, 3], dataIn[1][:, 3]))
- vy = np.concatenate((dataIn[0][0][:, 4], dataIn[0][1][:, 4], dataIn[1][:, 4]))
- vz = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5]))
- alpha = np.concatenate((alpha_up, alpha_lw, alpha_wk))
- dv = np.concatenate((dv_up, dv_lw, dv_wk))
- Me = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5]))
- rho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][:, 6], dataIn[1][:, 6]))
- dStarext = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][:, 7], dataIn[1][:, 7]))
- xxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][:, 8], dataIn[1][:, 8]))
-
- # Create dictionnaries for the fine and coarse meshes which are equivalent if meshMap is disabled.
- fMeshDict = {'globCoord': Mesh, 'yGlobCoord': Meshy, 'zGlobCoord': Meshz, 'bNodes': MeshbNodes, 'locCoord': xx,
- 'vtExt': vtExt, 'vx': vx, 'vy': vy, 'vz': vz, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
- cMeshDict = fMeshDict.copy()
- dataUpper = dataIn[0][0]
- dataLower = dataIn[0][1][1:, :]
- dataWake = dataIn[1]
-
- data = np.concatenate((dataUpper, dataLower, dataWake), axis=0)
-
- return fMeshDict, cMeshDict, data
-
-def getBLcoordinates(data):
- """Transform the reference coordinates into airfoil coordinates.
- """
- nN = len(data[:, 0])
- nE = nN-1 # Nbr of element along the surface.
- vt = np.zeros(nN) # Outer flow velocity.
- xx = np.zeros(nN) # Position along the surface of the airfoil.
- dx = np.zeros(nE) # Distance along two nodes.
- dv = np.zeros(nE) # Speed difference according to element.
- alpha = np.zeros(nE) # Angle of the element for the rotation matrix.
-
- for i in range(0, nE):
- alpha[i] = np.arctan2((data[i+1, 1] - data[i, 1]), (data[i+1, 0]-data[i, 0]))
- dx[i] = np.sqrt((data[i+1, 0] - data[i, 0]) **2 + (data[i+1, 1]-data[i, 1])**2)
- xx[i+1] = dx[i] + xx[i]
- # Tangential speed at node 1 element i
- vt[i] = (data[i, 3] * np.cos(alpha[i]) + data[i, 4] * np.sin(alpha[i]))
- # Tangential speed at node 2 element i
- vt[i+1] = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i]))
- # Velocity gradient according to element i.
- dv[i] = (vt[i+1] - vt[i]) / dx[i]
- return xx, dv, vt, nN, alpha
-
-def interpolateToFineMesh(cVector, fMesh, nFactor):
- """ Interpolate variables of cVector on the fine mesh fMesh.
-
- Params
- ------
- - cVector (numpy.ndarray): Distribution to be interpolated.
- - fMesh (numpy.ndarray): Local tangential coordinate (xi) of the fine mesh.
- - nFactor (int): Number of times each cell was split when creating the fine mesh.
- """
- fVector = np.zeros(len(fMesh)-1)
- for cPoint in range(len(cVector) - 1):
- dv = cVector[cPoint + 1] - cVector[cPoint]
- for fPoint in range(nFactor):
- fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + \
- fPoint * dv / nFactor
- fVector = np.append(fVector, cVector[-1])
- return fVector
-
-def createFineMesh(cMesh, nFactor):
- """ Create fine mesh distribution in the local frame of reference.
-
- Params
- ------
- - cMesh (numpy.ndarray): Initial mesh (used for the inviscid calculations.
- - nFactor (int): Number of times each cell will be split to create the fine mesh.
- """
- fMesh = []
- for iPoint in range(len(cMesh) - 1):
- dx = cMesh[iPoint + 1] - cMesh[iPoint]
- for iRefinedPoint in range(nFactor):
- fMesh = np.append(
- fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
- fMesh = np.append(fMesh, cMesh[-1])
- return fMesh
diff --git a/vii/interfaces/viscous/DWake.py b/vii/interfaces/viscous/DWake.py
deleted file mode 100755
index 9e770a9..0000000
--- a/vii/interfaces/viscous/DWake.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-## 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
diff --git a/vii/src/CMakeLists.txt b/vii/src/CMakeLists.txt
deleted file mode 100644
index 1e17fcd..0000000
--- a/vii/src/CMakeLists.txt
+++ /dev/null
@@ -1,27 +0,0 @@
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-# CMake input file of dart.so
-
-FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
-
-ADD_LIBRARY(vii SHARED ${SRCS})
-MACRO_DebugPostfix(vii)
-TARGET_INCLUDE_DIRECTORIES(vii PUBLIC ${PROJECT_SOURCE_DIR}/vii/src)
-
-TARGET_LINK_LIBRARIES(vii tbox)
-
-INSTALL(TARGETS vii DESTINATION ${CMAKE_INSTALL_PREFIX})
-
-SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
diff --git a/vii/src/DBoundaryLayer.cpp b/vii/src/DBoundaryLayer.cpp
deleted file mode 100644
index a690d71..0000000
--- a/vii/src/DBoundaryLayer.cpp
+++ /dev/null
@@ -1,185 +0,0 @@
-#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
diff --git a/vii/src/DBoundaryLayer.h b/vii/src/DBoundaryLayer.h
deleted file mode 100644
index 5d3853c..0000000
--- a/vii/src/DBoundaryLayer.h
+++ /dev/null
@@ -1,69 +0,0 @@
-#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
diff --git a/vii/src/DClosures.cpp b/vii/src/DClosures.cpp
deleted file mode 100644
index badb021..0000000
--- a/vii/src/DClosures.cpp
+++ /dev/null
@@ -1,256 +0,0 @@
-#include "DClosures.h"
-#include <cmath>
-#include <iostream>
-
-using namespace vii;
-
-/**
- * @brief Laminar closure relations at a given station.
- *
- * @param closureVars Vector where the computed variables will be stored.
- * @param theta Momentum thickness.
- * @param H Shape factor of the boundary layer.
- * @param Ret Reynolds number based on the momentum thickness.
- * @param Me Local Mach number.
- * @param name Name of the region.
- */
-void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
- double H, double Ret, const double Me,
- const std::string name)
-{
- H = std::max(H, 1.00005);
-
- // Kinematic shape factor H*.
- double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
- if (name == "wake")
- Hk = std::max(Hk, 1.00005);
- else
- Hk = std::max(Hk, 1.05000);
-
- // Density shape parameter.
- double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
-
- // Boundary layer thickness.
- double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
-
- double Hstar = 0.;
- double Hks = Hk - 4.35;
- if (Hk <= 4.35)
- Hstar = 0.0111 * Hks * Hks / (Hk + 1) - 0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 - 0.0002 * (Hks * Hk) * (Hks * Hk);
- else
- Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
-
- // Whitfield's minor additional compressibility correction.
- Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
-
- // Friction coefficient.
- double tmp = 0.;
- double cf = 0.;
- if (Hk < 5.5)
- {
- tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
- cf = (-0.07 + 0.0727 * tmp) / Ret;
- }
- else if (Hk >= 5.5)
- {
- tmp = 1.0 - 1.0 / (Hk - 4.5);
- cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
- }
-
- // Dissipation coefficient.
- double Cd = 0.;
- if (Hk < 4)
- Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
- else
- {
- double HkCd = (Hk - 4) * (Hk - 4);
- double denCd = 1 + 0.02 * HkCd;
- Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
- }
-
- // Wake relations.
- if (name == "wake")
- {
- Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2 * Ret));
- cf = 0.;
- }
-
- double us = 0.;
- double ctEq = 0.;
-
- closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
-}
-
-/**
- * @brief Turbulent closure relations at a given station.
- *
- * @param closureVars Vector where the computed variables will be stored.
- * @param theta Momentum thickness.
- * @param H Shape factor of the boundary layer.
- * @param Ct Shear stress coefficient.
- * @param Ret Reynolds number based on the momentum thickness.
- * @param Me Local Mach number.
- * @param name Name of the region.
- */
-void Closures::turbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name)
-{
- H = std::max(H, 1.00005);
- Ct = std::min(Ct, 0.3);
- // Ct = std::max(std::min(Ct, 0.30), 0.0000001);
-
- double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
- if (name == "wake")
- Hk = std::max(Hk, 1.00005);
- else
- Hk = std::max(Hk, 1.05000);
-
- double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
- double gamma = 1.4 - 1.;
- double Fc = std::sqrt(1 + 0.5 * gamma * Me * Me);
-
- double H0 = 0.;
- if (Ret <= 400)
- H0 = 4.0;
- else
- H0 = 3 + (400 / Ret);
- if (Ret <= 200)
- Ret = 200;
-
- double Hstar = 0.;
- if (Hk <= H0)
- Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
- else
- Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
-
- // Whitfield's minor additional compressibility correction.
- Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
-
- double logRt = std::log(Ret / Fc);
- logRt = std::max(logRt, 3.0);
- double arg = std::max(-1.33 * Hk, -20.0);
-
- // Equivalent normalized wall slip velocity.
- double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H));
-
- // Boundary layer thickness.
- double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
-
- double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
- double Hkc = 0.;
- double Cdw = 0.;
- double Cdd = 0.;
- double Cdl = 0.;
-
- double Hmin = 0.;
- double Fl = 0.;
- double Dfac = 0.;
-
- double cf = 0.;
- double Cd = 0.;
-
- double n = 1.0;
-
- double ctEq = 0.;
-
- if (name == "wake")
- {
- if (us > 0.99995)
- us = 0.99995;
- n = 2.0; // two wake halves
- cf = 0.0; // no friction inside the wake
- Hkc = Hk - 1;
- Cdw = 0.0; // Wall contribution.
- Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution.
- Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; // Laminar stress contribution to outer layer.
- ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5.
- }
- else
- {
- if (us > 0.95)
- us = 0.98;
- n = 1.0;
- cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4 - (Hk / 0.875)) - 1));
- Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
- // Dissipation coefficient.
- Hmin = 1 + 2.1 / std::log(Ret);
- Fl = (Hk - 1) / (Hmin - 1);
- Dfac = 0.5 + 0.5 * std::tanh(Fl);
- Cdw = 0.5 * (cf * us) * Dfac;
- Cdd = (0.995 - us) * Ct * Ct;
- Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret;
- ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
- }
- if (n * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0)
- std::cout << "Negative sqrt encountered " << std::endl;
-
- // Dissipation coefficient.
- Cd = n * (Cdw + Cdd + Cdl);
- closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
-}
-
-/**
- * @brief Turbulent closure relations at a given station.
- *
- * @param closureVars Computed equilibrium shear stress coefficient.
- * @param theta Momentum thickness.
- * @param H Shape factor of the boundary layer.
- * @param Ct Shear stress coefficient.
- * @param Ret Reynolds number based on the momentum thickness.
- * @param Me Local Mach number.
- * @param name Name of the region.
- */
-void Closures::turbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name)
-{
- H = std::max(H, 1.00005);
- Ct = std::min(Ct, 0.3);
-
- double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
- if (name == "wake")
- Hk = std::max(Hk, 1.00005);
- else
- Hk = std::max(Hk, 1.05000);
-
- double H0 = 0.;
- if (Ret <= 400)
- H0 = 4.0;
- else
- H0 = 3 + (400 / Ret);
- if (Ret <= 200)
- Ret = 200;
-
- double Hstar = 0.;
- if (Hk <= H0)
- Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
- else
- Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
-
- // Whitfield's minor additional compressibility correction.
- Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
-
- // Equivalent normalized wall slip velocity.
- double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H));
-
- double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
-
- double Hkc = 0.;
- double ctEq = 0.;
-
- if (name == "wake")
- {
- if (us > 0.99995)
- us = 0.99995;
- Hkc = Hk - 1;
- ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
- }
- else
- {
- if (us > 0.95)
- us = 0.98;
- Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
- ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
- }
- if (Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0)
- std::cout << "Negative sqrt encountered " << std::endl;
-
- closureVars = ctEq;
-}
\ No newline at end of file
diff --git a/vii/src/DClosures.h b/vii/src/DClosures.h
deleted file mode 100644
index fa4aca2..0000000
--- a/vii/src/DClosures.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifndef DCLOSURES_H
-#define DCLOSURES_H
-
-#include "vii.h"
-#include <vector>
-#include <string>
-
-namespace vii
-{
-
-/**
- * @brief Boundary layer closure relations.
- * @author Paul Dechamps
- */
-class VII_API Closures
-{
-
-public:
- static void laminarClosures(std::vector<double> &closureVars, double theta,
- double H, double Ret, const double Me,
- const std::string name);
- static void turbulentClosures(std::vector<double> &closureVars, double theta,
- double H, double Ct, double Ret, double Me,
- const std::string name);
- static void turbulentClosures(double &closureVars, double theta,
- double H, double Ct, double Ret, const double Me,
- const std::string name);
-};
-} // namespace vii
-#endif // DCLOSURES_H
diff --git a/vii/src/DDiscretization.cpp b/vii/src/DDiscretization.cpp
deleted file mode 100644
index 872c555..0000000
--- a/vii/src/DDiscretization.cpp
+++ /dev/null
@@ -1,76 +0,0 @@
-#include "DDiscretization.h"
-#include <cmath>
-#include <algorithm>
-
-using namespace vii;
-
-Discretization::Discretization()
-{
- nMarkers = 0;
- nElm = 0;
- prevnMarkers = 0;
-}
-
-Discretization::~Discretization() {}
-
-void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z)
-{
- if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size())
- throw std::runtime_error("vii::Discretization Wrong mesh sizes.\n");
-
- nMarkers = _x.size();
- x.resize(nMarkers, 0.);
- xoc.resize(nMarkers, 0.);
- y.resize(nMarkers, 0.);
- z.resize(nMarkers, 0.);
- nElm = nMarkers - 1;
- x = _x;
- y = _y;
- z = _z;
-
- computeBLcoord();
- computeOCcoord();
-}
-
-/**
- * @brief Compute x/chord coordinate.
- */
-void Discretization::computeOCcoord()
-{
- // Find min and max of array x
- auto min_it = std::min_element(x.begin(), x.end());
- double xMin = *min_it;
- auto max_it = std::max_element(x.begin(), x.end());
- double xMax = *max_it;
-
- // Compute xoc
- chord = xMax - xMin;
- if (chord <= 0)
- throw;
- for (size_t iPoint = 0; iPoint < x.size(); ++iPoint)
- xoc[iPoint] = (x[iPoint] - xMin) / chord;
-}
-
-/**
- * @brief Compute nodes coordinates in the local frame of reference.
- */
-void Discretization::computeBLcoord()
-{
- loc.resize(nMarkers, 0.);
- alpha.resize(nElm, 0.);
- dx.resize(nElm, 0.);
-
- for (size_t iPoint = 0; iPoint < nElm; ++iPoint)
- {
- alpha[iPoint] = std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]);
- // dx = sqrt(delta x ^2 + delta y ^2).
- dx[iPoint] = std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint]));
- loc[iPoint + 1] = loc[iPoint] + dx[iPoint];
- }
-}
-
-void Discretization::setExt(std::vector<double> extM)
-{
- ext.resize(extM.size(), 0.);
- ext = extM;
-};
diff --git a/vii/src/DDiscretization.h b/vii/src/DDiscretization.h
deleted file mode 100644
index 1f40cb8..0000000
--- a/vii/src/DDiscretization.h
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef DDISCRETIZATION_H
-#define DDISCRETIZATION_H
-
-#include "vii.h"
-#include "DBoundaryLayer.h"
-
-namespace vii
-{
-
-/**
- * @brief Boundary layer mesh.
- * @author Paul Dechamps
- */
-class VII_API Discretization
-{
-private:
- std::vector<double> x; ///< x coordinate in the global frame of reference.
- std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
- std::vector<double> y; ///< y coordinate in the global frame of reference.
- std::vector<double> z; ///< z coordinate in the global frame of reference.
- std::vector<double> loc; ///< xi coordinate in the local frame of reference.
- std::vector<double> ext; ///< xi coordinate in the local frame of reference, at the edge of the boundary layer.
- std::vector<double> dx; ///< Cell size.
- std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
- size_t nMarkers; ///< Number of points in the domain.
- size_t nElm; ///< Number of cells in the domain.
- double chord; ///< Chord of the section.
-
-public:
- size_t prevnMarkers; ///< Number of points in the domain at the previous iteration.
-
- Discretization();
- ~Discretization();
-
- // Getters.
- size_t getnMarkers() const { return nMarkers; };
- double getLoc(size_t iPoint) const { return loc[iPoint]; };
- double getx(size_t iPoint) const { return x[iPoint]; };
- double getxoc(size_t iPoint) const { return xoc[iPoint]; };
- double gety(size_t iPoint) const { return y[iPoint]; };
- double getz(size_t iPoint) const { return z[iPoint]; };
- double getExt(size_t iPoint) const { return ext[iPoint]; };
- double getAlpha(size_t iElm) const { return alpha[iElm]; };
- size_t getDiffnElm() const { return nMarkers - prevnMarkers; };
- double getChord() const { return chord; };
-
- // Setters.
- void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z);
- void setExt(std::vector<double> extM);
- void computeBLcoord();
- void computeOCcoord();
-};
-} // namespace vii
-#endif // WDISCRETIZATION_H
diff --git a/vii/src/DDriver.cpp b/vii/src/DDriver.cpp
deleted file mode 100644
index 0828044..0000000
--- a/vii/src/DDriver.cpp
+++ /dev/null
@@ -1,799 +0,0 @@
-#include "DDriver.h"
-#include "DBoundaryLayer.h"
-#include "DSolver.h"
-#include <iomanip>
-
-#define ANSI_COLOR_RED "\x1b[1;31m"
-#define ANSI_COLOR_GREEN "\x1b[1;32m"
-#define ANSI_COLOR_YELLOW "\x1b[1;33m"
-#define ANSI_COLOR_RESET "\x1b[0m"
-
-using namespace vii;
-
-/**
- * @brief Driver of the viscous solver.
- *
- * @param _Re Freestream Reynolds number.
- * @param _Minf Freestream Mach number.
- * @param _CFL0 Initial CFL number for pseudo-time integration.
- * @param nSections Number of sections in the computation (1 for 2D cases).
- * @param xtrF_top Forced transition location on the upper side.
- * @param xtrF_bot Forced transition location on the lower side.
- * @param _Span Span of the considered wing (only used for drag computation).
- * @param _verbose Verbosity level of the solver 0<=verbose<=1.
- */
-Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
- double xtrF_top, double xtrF_bot, double _span, unsigned int _verbose)
-{
- // Freestream parameters.
- Re = _Re;
- CFL0 = _CFL0;
- verbose = _verbose;
-
- // Vector of forced transition {upper, lower, wake}.
- std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.};
-
- // Initialzing boundary layer
- sections.resize(nSections);
- convergenceStatus.resize(nSections);
- for (size_t iSec = 0; iSec < nSections; ++iSec)
- {
- convergenceStatus[iSec].resize(3);
- for (size_t i = 0; i < 3; ++i)
- {
- BoundaryLayer *region = new BoundaryLayer(xtrF[i]);
- sections[iSec].push_back(region);
- }
- sections[iSec][0]->name = "upper";
- sections[iSec][1]->name = "lower";
- sections[iSec][2]->name = "wake";
- }
-
- Cdt_sec.resize(sections.size(), 0.);
- Cdf_sec.resize(sections.size(), 0.);
- Cdp_sec.resize(sections.size(), 0.);
- Cdt = 0.;
- Cdf = 0.;
- Cdp = 0.;
- span = _span;
-
- // Pre/post processing flags.
- stagPtMvmt.resize(sections.size(), false);
-
- // Time solver initialization.
- tSolver = new Solver(_CFL0, _Minf, Re);
-
- // Numerical parameters.
- std::cout << "--- Viscous solver setup ---\n"
- << "Reynolds number: " << Re << " \n"
- << "Freestream Mach number: " << _Minf << " \n"
- << "Initial CFL number: " << CFL0 << " \n"
- << std::endl;
- std::cout << "Boundary layer initialized." << std::endl;
-}
-
-/**
- * @brief Driver and subsequent dependencies destruction.
- */
-Driver::~Driver()
-{
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- for (size_t i = 0; i < sections[iSec].size(); ++i)
- delete sections[iSec][i];
- delete tSolver;
- std::cout << "~vii::Driver()\n";
-} // end ~Driver
-
-/**
- * @brief Runs viscous operations. Temporal solver parametrisation is first adapted,
- * Solutions are initialized than time marched towards steady state successively for each points.
- *
- * @param couplIter Current coupling iteration.
- * @return int Solver exit code.
- */
-int Driver::run(unsigned int const couplIter)
-{
- bool lockTrans = false; // Flagging activation of transition routines.
- int pointExitCode; // Output of pseudo time integration (0: converged; -1: unsuccessful, failed; >0: unsuccessful nb iter).
-
- double maxMach = 0.;
-
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- {
- tms["0-CheckIn"].start();
- if (verbose > 1)
- std::cout << "Sec " << iSec << " ";
-
- // Check for stagnation point movement.
- if (couplIter != 0 && (sections[iSec][0]->mesh->getDiffnElm()))
- {
- stagPtMvmt[iSec] = true;
- if (verbose > 2)
- std::cout << "Stagnation point moved by " << sections[iSec][0]->mesh->getDiffnElm() << " elements." << std::endl;
- }
- else
- stagPtMvmt[iSec] = false;
-
- // Set boundary conditions.
- sections[iSec][0]->xtr = 1.; // Upper side initial xtr.
- sections[iSec][1]->xtr = 1.; // Lower side initial xtr.
- sections[iSec][2]->xtr = 0.; // Wake initial xtr (fully turbulent).
- sections[iSec][0]->setStagBC(Re);
- sections[iSec][1]->setStagBC(Re);
- tms["0-CheckIn"].stop();
-
- // Loop over the regions (upper, lower and wake).
- for (size_t iRegion = 0; iRegion < sections[iSec].size(); ++iRegion)
- {
- convergenceStatus[iSec][iRegion].resize(0);
-
- // Maximum Mach number in the region.
- if (sections[iSec][iRegion]->blEdge->getMaxMach() > maxMach)
- maxMach = sections[iSec][iRegion]->blEdge->getMaxMach();
-
- if (verbose > 1)
- std::cout << sections[iSec][iRegion]->name << " ";
-
- // Check if safe mode is required.
- if (couplIter == 0 || sections[iSec][iRegion]->mesh->getDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
- {
- tSolver->setCFL0(1);
- tSolver->setitFrozenJac(1);
- tSolver->setinitSln(true);
-
- if (sections[iSec][iRegion]->mesh->getDiffnElm())
- {
- std::vector<double> xxExtCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xxExtCurr.size(); ++iPoint)
- xxExtCurr[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
- sections[iSec][iRegion]->mesh->setExt(xxExtCurr);
-
- std::vector<double> DeltaStarZeros(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- sections[iSec][iRegion]->blEdge->setDeltaStar(DeltaStarZeros);
-
- std::vector<double> ueOld(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- for (size_t iPoint = 0; iPoint < ueOld.size(); ++iPoint)
- ueOld[iPoint] = sections[iSec][iRegion]->blEdge->getVt(iPoint);
- sections[iSec][iRegion]->blEdge->setVtOld(ueOld);
- }
-
- // Keyword annoncing iteration safety level.
- if (verbose > 1)
- std::cout << "restart. ";
- }
- else
- {
- tSolver->setCFL0(CFL0);
- tSolver->setitFrozenJac(5);
- tSolver->setinitSln(false);
-
- // Keyword annoncing iteration safety level.
- if (verbose > 1)
- std::cout << "update. ";
- }
-
- // Save number of points in each regions.
- sections[iSec][iRegion]->mesh->prevnMarkers = sections[iSec][iRegion]->mesh->getnMarkers();
-
- // Reset regime.
- lockTrans = false;
- if (iRegion < 2) // Airfoil
- for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
- sections[iSec][iRegion]->regime[k] = 0;
- else // Wake
- {
- for (size_t k = 0; k < sections[iSec][iRegion]->mesh->getnMarkers(); k++)
- sections[iSec][iRegion]->regime[k] = 1;
-
- // Impose wake boundary condition.
- std::vector<double> upperConditions(sections[iSec][iRegion]->getnVar(), 0.);
- std::vector<double> lowerConditions(sections[iSec][iRegion]->getnVar(), 0.);
- for (size_t k = 0; k < upperConditions.size(); ++k)
- {
- upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->mesh->getnMarkers() - 1) * sections[iSec][0]->getnVar() + k];
- lowerConditions[k] = sections[iSec][1]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k];
- }
- sections[iSec][iRegion]->setWakeBC(Re, upperConditions, lowerConditions);
- lockTrans = true;
- }
-
- // Loop over points
- for (size_t iPoint = 1; iPoint < sections[iSec][iRegion]->mesh->getnMarkers(); ++iPoint)
- {
- // Initialize solution at point
- tSolver->initialCondition(iPoint, sections[iSec][iRegion]);
- // Solve equations
- tms["1-Solver"].start();
- pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]);
- // Force transition if required by the user
- if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) >= sections[iSec][iRegion]->xtrF)
- sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] = 9.;
- tms["1-Solver"].stop();
-
- if (pointExitCode != 0)
- convergenceStatus[iSec][iRegion].push_back(iPoint);
-
- // Transition
- tms["2-Transition"].start();
- if (!lockTrans)
- {
- // Check if perturbation waves are growing or not.
- // Avoided in the case of a forced transition.
- if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) <= sections[iSec][iRegion]->xtrF)
- checkWaves(iPoint, sections[iSec][iRegion]);
-
- // Amplification factor is compared to critical amplification factor.
- if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit())
- {
- averageTransition(iPoint, sections[iSec][iRegion], 0);
- if (verbose > 1)
- {
- std::cout << std::fixed;
- std::cout << std::setprecision(2);
- std::cout << sections[iSec][iRegion]->xoctr
- << " (" << sections[iSec][iRegion]->transMarker << ") ";
- }
- lockTrans = true;
- }
- }
- tms["2-Transition"].stop();
- }
- tms["3-Blowing"].start();
- sections[iSec][iRegion]->computeBlwVel();
- tms["3-Blowing"].stop();
- }
- if (verbose > 1)
- std::cout << std::endl;
- }
-
- for (size_t i = 0; i < sections.size(); ++i)
- computeSectionalDrag(sections[i], i);
- computeTotalDrag();
-
- // Output information.
- solverExitCode = outputStatus(maxMach);
-
- return solverExitCode; // exit code
-}
-
-/**
- * @brief Check whether or not instabilities are growing in the laminar boundary layer.
- *
- * @param iPoint Marker id.
- * @param bl BoundaryLayer region.
- */
-void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl)
-{
-
- double logRet_crit = 2.492 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
-
- if (bl->Ret[iPoint] > 0) // Avoid log of negative number.
- {
- if (std::log10(bl->Ret[iPoint]) <= logRet_crit - 0.08)
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
- }
- else
- bl->u[iPoint * bl->getnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
-}
-
-/**
- * @brief Computes turbulent solution @ the transition point and averages solutions.
- *
- * @param iPoint Marker id.
- * @param bl BoundaryLayer region.
- * @param forced Flag for forced transition.
- */
-void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
-{
- // Averages solution on transition marker.
- size_t nVar = bl->getnVar();
-
- // Save laminar solution.
- std::vector<double> lamSol(nVar, 0.0);
- for (size_t k = 0; k < lamSol.size(); ++k)
- lamSol[k] = bl->u[iPoint * nVar + k];
-
- // Set up turbulent portion boundary condition.
- bl->transMarker = iPoint; // Save transition marker.
-
- // Loop over remaining points in the region.
- for (size_t k = iPoint; k < bl->mesh->getnMarkers(); ++k)
- bl->regime[k] = 1; // Set Turbulent regime.
-
- // Compute transition location.
- bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
- bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
-
- // Percentage of the subsection corresponding to a turbulent flow.
- double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
- double avLam = 1. - avTurb;
-
- // Impose boundary condition.
- double Cteq_trans;
- bl->closSolver->turbulentClosures(Cteq_trans, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
- bl->ctEq[iPoint - 1] = Cteq_trans;
- bl->u[(iPoint - 1) * nVar + 4] = 0.7 * bl->ctEq[iPoint - 1];
-
- // Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above).
- // These initial guess values do not influence the converged solution.
- bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
- bl->u[iPoint * nVar + 1] = 1.515; // Because we expect a significant drop of the shape factor.
-
- // Solve point in turbulent configuration.
- int exitCode = tSolver->integration(iPoint, bl);
-
- if (exitCode != 0 && verbose > 1)
- std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
-
- // Average both solutions (Now solution @ iPoint is the turbulent solution).
- bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
- bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
- bl->u[iPoint * nVar + 2] = 9.; // N.
- bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
- bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
-
- // Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
- std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->blEdge->getMe(iPoint - 1), bl->name);
- bl->Hstar[iPoint - 1] = lamParam[0];
- bl->Hstar2[iPoint - 1] = lamParam[1];
- bl->Hk[iPoint - 1] = lamParam[2];
- bl->cf[iPoint - 1] = lamParam[3];
- bl->cd[iPoint - 1] = lamParam[4];
- bl->delta[iPoint - 1] = lamParam[5];
- bl->ctEq[iPoint - 1] = lamParam[6];
- bl->us[iPoint - 1] = lamParam[7];
-
- std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, bl->u[(iPoint)*nVar + 0], bl->u[(iPoint)*nVar + 1], bl->u[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
- bl->Hstar[iPoint] = turbParam[0];
- bl->Hstar2[iPoint] = turbParam[1];
- bl->Hk[iPoint] = turbParam[2];
- bl->cf[iPoint] = turbParam[3];
- bl->cd[iPoint] = turbParam[4];
- bl->delta[iPoint] = turbParam[5];
- bl->ctEq[iPoint] = turbParam[6];
- bl->us[iPoint] = turbParam[7];
-}
-
-/**
- * @brief Friction, pressure and total drag computation qt each station of the configuration.
- *
- * @tparam Cdt = theta * Ue^((H+5)/2).
- * @tparam Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
- * @tparam Cdp = Cdtot - Cdf.
- *
- * @param bl BoundaryLayer region.
- * @param i Considered section.
- */
-void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i)
-{
- size_t nVar = bl[0]->getnVar();
- size_t lastWkPt = (bl[2]->mesh->getnMarkers() - 1) * nVar;
-
- // Normalize friction and dissipation coefficients.
- std::vector<double> frictioncoeff_upper(bl[0]->mesh->getnMarkers(), 0.0);
- for (size_t k = 0; k < frictioncoeff_upper.size(); ++k)
- frictioncoeff_upper[k] = bl[0]->blEdge->getVt(k) * bl[0]->blEdge->getVt(k) * bl[0]->cf[k];
-
- std::vector<double> frictioncoeff_lower(bl[1]->mesh->getnMarkers(), 0.0);
- for (size_t k = 0; k < frictioncoeff_lower.size(); ++k)
- frictioncoeff_lower[k] = bl[1]->blEdge->getVt(k) * bl[1]->blEdge->getVt(k) * bl[1]->cf[k];
-
- // Compute friction drag coefficient.
- double Cdf_upper = 0.0;
- for (size_t k = 1; k < bl[0]->mesh->getnMarkers(); ++k)
- Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->getx(k) - bl[0]->mesh->getx(k - 1)) / 2;
- double Cdf_lower = 0.0;
- for (size_t k = 1; k < bl[1]->mesh->getnMarkers(); ++k)
- Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->getx(k) - bl[1]->mesh->getx(k - 1)) / 2;
-
- // Friction drag coefficient.
- Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
- // Total drag coefficient.
- Cdt_sec[i] = 2 * bl[2]->u[lastWkPt + 0] * pow(bl[2]->u[lastWkPt + 3], (bl[2]->u[lastWkPt + 1] + 5) / 2);
- // Pressure drag coefficient.
- Cdp_sec[i] = Cdt_sec[i] - Cdf_sec[i];
-}
-
-/**
- * @brief Compute total drag coefficient on the airfoil/wing.
- * Performs the sectional drag integration for 3D computations.
- *
- */
-void Driver::computeTotalDrag()
-{
- Cdt = 0.;
- Cdf = 0.;
- Cdp = 0.;
-
- if (sections.size() == 1 || span == 0.)
- {
- Cdt = Cdt_sec[0];
- Cdf = Cdf_sec[0];
- Cdp = Cdp_sec[0];
- }
- else
- {
- Cdt += (sections[0][0]->mesh->getz(0) - 0) * Cdt_sec[0];
- Cdp += (sections[0][0]->mesh->getz(0) - 0) * Cdp_sec[0];
- Cdf += (sections[0][0]->mesh->getz(0) - 0) * Cdf_sec[0];
- double dz = 0.;
- for (size_t k = 1; k < sections.size(); ++k)
- {
- dz = (sections[k][0]->mesh->getz(0) - sections[k - 1][0]->mesh->getz(0));
- Cdt += dz * Cdt_sec[k];
- Cdp += dz * Cdp_sec[k];
- Cdf += dz * Cdf_sec[k];
- }
- Cdt /= span;
- Cdp /= span;
- Cdf /= span;
- }
-}
-
-/**
- * @brief Outputs solver status depending on the verbosity level (verbose).
- *
- * @param maxMach Maximum Mach number identified on the configuration.
- * @return int Solver exit code (0 converged, !=0 not converged).
- */
-int Driver::outputStatus(double maxMach) const
-{
- int solverExitCode = 0;
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg)
- if (convergenceStatus[iSec][iReg].size() != 0)
- {
- solverExitCode = -1;
- break;
- }
-
- if (verbose > 2 && solverExitCode != 0)
- {
- // Output unconverged points.
- std::cout << "Unconverged points" << std::endl;
- std::cout << std::setw(10) << "Section"
- << std::setw(12) << "Region"
- << std::setw(16) << "Markers" << std::endl;
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- {
- unsigned int count = 0;
- for (size_t iReg = 0; iReg < sections[iSec].size(); ++iReg)
- {
- std::cout << std::setw(10) << iSec;
- std::cout << std::setw(12) << iReg << " ";
- if (convergenceStatus[iSec][iReg].size() != 0)
- for (size_t iPoint = 0; iPoint < convergenceStatus[iSec][iReg].size(); ++iPoint)
- {
- std::cout << convergenceStatus[iSec][iReg][iPoint] << " (" << sections[iSec][iReg]->mesh->getx(convergenceStatus[iSec][iReg][iPoint]) << "), ";
- ++count;
- }
- std::cout << "\n";
- }
- std::cout << "total: " << count << "\n";
- }
- std::cout << "" << std::endl;
- }
-
- if (verbose > 0)
- {
- // Compute mean transition locations.
- std::vector<double> meanTransitions(2, 0.);
- for (size_t iSec = 0; iSec < sections.size(); ++iSec)
- for (size_t iReg = 0; iReg < 2; ++iReg)
- meanTransitions[iReg] += sections[iSec][iReg]->xoctr;
- for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg)
- meanTransitions[iReg] /= sections.size();
-
- std::cout << std::fixed << std::setprecision(2)
- << "Viscous solver for Re " << Re / 1e6 << "e6, Mach max "
- << maxMach
- << std::endl;
-
- std::cout << std::setw(10) << "xTr Upper"
- << std::setw(12) << "xTr Lower"
- << std::setw(12) << "Cd" << std::endl;
- std::cout << std::fixed << std::setprecision(2);
- std::cout << std::setw(10) << meanTransitions[0]
- << std::setw(12) << meanTransitions[1];
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(12) << Cdt << "\n";
-
- if (verbose > 0)
- {
- if (solverExitCode == 0)
- std::cout << ANSI_COLOR_GREEN << "Viscous solver converged." << ANSI_COLOR_RESET << std::endl;
- else
- std::cout << ANSI_COLOR_YELLOW << "Viscous solver not converged." << ANSI_COLOR_RESET << std::endl;
- }
-
- if (verbose > 2)
- {
- std::cout << "Viscous solver CPU" << std::endl
- << tms;
- }
- }
-
- if (verbose > 0)
- std::cout << "" << std::endl;
- return solverExitCode;
-}
-
-/**
- * @brief Set nodes coordinates in the global frame of reference.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _x Vector of x coordinates of the nodes in the region.
- * @param _y Vector of y coordinates of the nodes in the region.
- * @param _z Vector of z coordinates of the nodes in the region.
- */
-void Driver::setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
-{
- sections[iSec][iRegion]->setMesh(_x, _y, _z);
-}
-
-/**
- * @brief Set the velocity at the nodes in the global frame of reference.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param vx Vector of velocities projected on the x axis of the global frame of reference.
- * @param vy Vector of velocities projected on the y axis of the global frame of reference.
- * @param vz Vector of velocities projected on the z axis of the global frame of reference.
- */
-void Driver::setVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
-{
- if (vx.size() != vy.size() || vy.size() != vz.size() || vx.size() != vz.size())
- throw std::runtime_error("dart::Driver Wrong velocity vector sizes.");
- if (vx.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("dart::Driver Velocity vector is not coherent with the mesh.");
-
- std::vector<double> vt(vx.size(), 0.0);
-
- for (size_t iPoint = 0; iPoint < vx.size() - 1; ++iPoint)
- {
- vt[iPoint] = vx[iPoint] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
- vt[iPoint + 1] = vx[iPoint + 1] * std::cos(sections[iSec][iRegion]->mesh->getAlpha(iPoint)) + vy[iPoint + 1] * std::sin(sections[iSec][iRegion]->mesh->getAlpha(iPoint));
- }
-
- sections[iSec][iRegion]->blEdge->setVt(vt);
-}
-
-/**
- * @brief Set the Mach number at the nodes.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _Me Vector of Mach numbers at the nodes.
- */
-void Driver::setMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
-{
- if (_Me.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("dart::Driver Mach number vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setMe(_Me);
-}
-
-/**
- * @brief Set the density at the nodes.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _Rhoe Vector of density at the nodes.
- */
-void Driver::setRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
-{
- if (_Rhoe.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("dart::Driver Density vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setRhoe(_Rhoe);
-}
-
-/**
- * @brief Set the nodes coordinates in the local frame of reference at the edge of the boundary layer.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _xxExt Vector of coordinates of the nodes.
- */
-void Driver::setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
-{
- if (_xxExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("dart::Driver External mesh vector is not coherent with the mesh.");
- sections[iSec][iRegion]->mesh->setExt(_xxExt);
-}
-
-/**
- * @brief Set the displacement thickness (previous iteration).
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _DeltaStarExt Vector of displacement thicknesses at the nodes.
- */
-void Driver::setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
-{
- if (_DeltaStarExt.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- {
- std::cout << "size DeltaStar: " << _DeltaStarExt.size() << ". nMarkers: " << sections[iSec][iRegion]->mesh->getnMarkers() << std::endl;
- throw std::runtime_error("dart::Driver External delta star vector is not coherent with the mesh.");
- }
- sections[iSec][iRegion]->blEdge->setDeltaStar(_DeltaStarExt);
-}
-
-/**
- * @brief Returns x coordinates of the nodes in the local frame of reference in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractxCoord(size_t iSec, size_t iRegion) const
-{
- std::vector<double> xCoord(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xCoord.size(); ++iPoint)
- xCoord[iPoint] = sections[iSec][iRegion]->mesh->getx(iPoint);
- return xCoord;
-}
-
-/**
- * @brief Returns blowing velocities (defined on the cells' cg) in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractBlowingVel(size_t iSec, size_t iRegion) const
-{
- return sections[iSec][iRegion]->blowingVelocity;
-}
-std::vector<double> Driver::extractxx(size_t iSec, size_t iRegion) const
-{
- std::vector<double> xx(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
- for (size_t iPoint = 0; iPoint < xx.size(); ++iPoint)
- xx[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
- return xx;
-}
-
-/**
- * @brief Returns the displacement thickness in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractDeltaStar(size_t iSec, size_t iRegion) const
-{
- return sections[iSec][iRegion]->deltaStar;
-}
-
-/**
- * @brief Returns the velocity of the interaction law in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @return std::vector<double>
- */
-std::vector<double> Driver::extractUe(size_t iSec, size_t iRegion) const
-{
- std::vector<double> ueCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.0);
- for (size_t i = 0; i < ueCurr.size(); ++i)
- ueCurr[i] = sections[iSec][iRegion]->u[i * sections[iSec][iRegion]->getnVar() + 3];
- return ueCurr;
-}
-
-/**
- * @brief Set the velocity at the edge of the BL at the previous iteration in the considered region.
- *
- * @param iSec id of the considered section.
- * @param iRegion id of the considered region.
- * @param _ue Velocity vector.
- * @return std::vector<double>
- */
-void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
-{
- if (_ue.size() != sections[iSec][iRegion]->mesh->getnMarkers())
- throw std::runtime_error("dart::Driver External velocity vector is not coherent with the mesh.");
- sections[iSec][iRegion]->blEdge->setVtOld(_ue);
-}
-
-/**
- * @brief Returns the integral boundary layer solution in a section (sorted from upper Te to lower Te than wake).
- *
- * @param iSec id of the considered section.
- * @return std::vector<std::vector<double>>
- */
-std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
-{
- size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers();
- size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers(); // No stagnation point doublon.
- size_t nVar = sections[iSec][0]->getnVar();
-
- std::vector<std::vector<double>> Sln(20);
- for (size_t iVector = 0; iVector < Sln.size(); ++iVector)
- Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.);
- // Upper side.
- size_t revPoint = nMarkersUp - 1;
- for (size_t iPoint = 0; iPoint < nMarkersUp; ++iPoint)
- {
- Sln[0][iPoint] = sections[iSec][0]->u[revPoint * nVar + 0];
- Sln[1][iPoint] = sections[iSec][0]->u[revPoint * nVar + 1];
- Sln[2][iPoint] = sections[iSec][0]->u[revPoint * nVar + 2];
- Sln[3][iPoint] = sections[iSec][0]->u[revPoint * nVar + 3];
- Sln[4][iPoint] = sections[iSec][0]->u[revPoint * nVar + 4];
- Sln[5][iPoint] = sections[iSec][0]->deltaStar[revPoint];
-
- Sln[6][iPoint] = sections[iSec][0]->Ret[revPoint];
- Sln[7][iPoint] = sections[iSec][0]->cd[revPoint];
- Sln[8][iPoint] = sections[iSec][0]->cf[revPoint];
- Sln[9][iPoint] = sections[iSec][0]->Hstar[revPoint];
- Sln[10][iPoint] = sections[iSec][0]->Hstar2[revPoint];
- Sln[11][iPoint] = sections[iSec][0]->Hk[revPoint];
- Sln[12][iPoint] = sections[iSec][0]->ctEq[revPoint];
- Sln[13][iPoint] = sections[iSec][0]->us[revPoint];
- Sln[14][iPoint] = sections[iSec][0]->delta[revPoint];
-
- Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint);
- Sln[16][iPoint] = sections[iSec][0]->mesh->getxoc(revPoint);
- Sln[17][iPoint] = sections[iSec][0]->mesh->gety(revPoint);
- Sln[18][iPoint] = sections[iSec][0]->mesh->getz(revPoint);
- Sln[19][iPoint] = sections[iSec][0]->blEdge->getVt(revPoint);
-
- --revPoint;
- }
- // Lower side.
- for (size_t iPoint = 0; iPoint < sections[iSec][1]->mesh->getnMarkers(); ++iPoint)
- {
- Sln[0][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 0];
- Sln[1][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 1];
- Sln[2][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 2];
- Sln[3][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 3];
- Sln[4][nMarkersUp + iPoint] = sections[iSec][1]->u[iPoint * nVar + 4];
- Sln[5][nMarkersUp + iPoint] = sections[iSec][1]->deltaStar[iPoint];
-
- Sln[6][nMarkersUp + iPoint] = sections[iSec][1]->Ret[iPoint];
- Sln[7][nMarkersUp + iPoint] = sections[iSec][1]->cd[iPoint];
- Sln[8][nMarkersUp + iPoint] = sections[iSec][1]->cf[iPoint];
- Sln[9][nMarkersUp + iPoint] = sections[iSec][1]->Hstar[iPoint];
- Sln[10][nMarkersUp + iPoint] = sections[iSec][1]->Hstar2[iPoint];
- Sln[11][nMarkersUp + iPoint] = sections[iSec][1]->Hk[iPoint];
- Sln[12][nMarkersUp + iPoint] = sections[iSec][1]->ctEq[iPoint];
- Sln[13][nMarkersUp + iPoint] = sections[iSec][1]->us[iPoint];
- Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint];
-
- Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint);
- Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getxoc(iPoint);
- Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->mesh->gety(iPoint);
- Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getz(iPoint);
- Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->blEdge->getVt(iPoint);
- }
- // Wake.
- for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint)
- {
- Sln[0][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 0];
- Sln[1][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 1];
- Sln[2][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 2];
- Sln[3][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 3];
- Sln[4][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->u[iPoint * nVar + 4];
- Sln[5][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->deltaStar[iPoint];
-
- Sln[6][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Ret[iPoint];
- Sln[7][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cd[iPoint];
- Sln[8][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->cf[iPoint];
- Sln[9][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar[iPoint];
- Sln[10][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hstar2[iPoint];
- Sln[11][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->Hk[iPoint];
- Sln[12][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->ctEq[iPoint];
- Sln[13][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->us[iPoint];
- Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint];
-
- Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint);
- Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getxoc(iPoint);
- Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->gety(iPoint);
- Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getz(iPoint);
- Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->blEdge->getVt(iPoint);
- }
-
- if (verbose > 2)
- std::cout << "Solution structure sent." << std::endl;
- return Sln;
-}
\ No newline at end of file
diff --git a/vii/src/DDriver.h b/vii/src/DDriver.h
deleted file mode 100644
index eddf310..0000000
--- a/vii/src/DDriver.h
+++ /dev/null
@@ -1,76 +0,0 @@
-#ifndef DDRIVER_H
-#define DDRIVER_H
-
-#include "vii.h"
-#include "wObject.h"
-#include "DBoundaryLayer.h"
-#include "wTimers.h"
-namespace vii
-{
-
-/**
- * @brief Boundary layer solver driver. Performs the space marching and set up time marching for each point.
- * @authors Paul Dechamps
- */
-class VII_API Driver : public fwk::wSharedObject
-{
-public:
- double Cdt; ///< Total drag coefficient.
- double Cdf; ///< Total friction drag coefficient.
- double Cdp; ///< Total pressure drag coefficient.
- std::vector<double> Cdt_sec; ///< Sectional total drag coefficient.
- std::vector<double> Cdf_sec; ///< Sectional friction drag coefficient.
- std::vector<double> Cdp_sec; ///< Sectional pressure drag coefficient.
- std::vector<std::vector<BoundaryLayer *>> sections; ///< Boundary layer regions sorted by section.
- unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1.
-
-private:
- double Re; ///< Freestream Reynolds number.
- double CFL0; ///< Initial CFL number for pseudo time integration.
- double span; ///< Wing Span (Used only for drag computation, not used if 2D case).
- std::vector<bool> stagPtMvmt; ///< Flag to account for stagnation point movements.
- Solver *tSolver; ///< Pseudo-time solver.
- int solverExitCode; ///< Exit code of viscous calculations.
- std::vector<std::vector<std::vector<size_t>>> convergenceStatus; ///< Vector containing points that did not converged.
-
-protected:
- fwk::Timers tms; ///< Internal timers
-
-public:
- Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
- double xtrF_top = -1., double xtrF_bot = -1., double _span = 0., unsigned int verbose = 1);
- ~Driver();
- int run(unsigned int const couplIter);
-
- // Getters.
- double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
- double getRe() const { return Re; };
- double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
- std::vector<std::vector<double>> getSolution(size_t iSec);
-
- // Setters.
- void setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
- void setVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz);
- void setMe(size_t iSec, size_t iRegion, std::vector<double> _Me);
- void setRhoe(size_t iSec, size_t iRegion, std::vector<double> _rhoe);
- void setxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt);
- void setDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _deltaStarExt);
- void setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue);
-
- // Extractors.
- std::vector<double> extractBlowingVel(size_t iSec, size_t iRegion) const;
- std::vector<double> extractxx(size_t iSec, size_t iRegion) const;
- std::vector<double> extractDeltaStar(size_t iSec, size_t iRegion) const;
- std::vector<double> extractxCoord(size_t iSec, size_t iRegion) const;
- std::vector<double> extractUe(size_t iSec, size_t iRegion) const;
-
-private:
- void checkWaves(size_t iPoint, BoundaryLayer *bl);
- void averageTransition(size_t iPoint, BoundaryLayer *bl, int forced);
- void computeSectionalDrag(std::vector<BoundaryLayer *> const &bl, size_t i);
- void computeTotalDrag();
- void computeBlwVel();
- int outputStatus(double maxMach) const;
-};
-} // namespace vii
-#endif // DDRIVER_H
diff --git a/vii/src/DEdge.cpp b/vii/src/DEdge.cpp
deleted file mode 100644
index 48bde11..0000000
--- a/vii/src/DEdge.cpp
+++ /dev/null
@@ -1,59 +0,0 @@
-#include "DEdge.h"
-
-using namespace vii;
-
-/**
- * @brief Construct a new Edge::Edge object.
- *
- */
-Edge::Edge()
-{
- Me.resize(0, 0.);
- vt.resize(0, 0.);
- rhoe.resize(0, 0.);
- deltaStar.resize(0, 0.);
- vtOld.resize(0, 0.);
-}
-
-/**
- * @brief Return the maximum inviscid Mach number in the region.
- *
- * @return double
- */
-double Edge::getMaxMach() const
-{
- if (Me.size() <= 0)
- throw std::runtime_error("vii::Edge Invalid Mach number vector.");
- double Mmax = 0.0;
- for (size_t i = 0; i < Me.size(); ++i)
- if (Me[i] > Mmax)
- Mmax = Me[i];
- return Mmax;
-}
-
-// Setters.
-void Edge::setMe(std::vector<double> const _Me)
-{
- Me.resize(_Me.size(), 0.);
- Me = _Me;
-}
-void Edge::setVt(std::vector<double> const _vt)
-{
- vt.resize(_vt.size(), 0.);
- vt = _vt;
-}
-void Edge::setRhoe(std::vector<double> const _rhoe)
-{
- rhoe.resize(_rhoe.size(), 0.);
- rhoe = _rhoe;
-}
-void Edge::setDeltaStar(std::vector<double> const _deltaStar)
-{
- deltaStar.resize(_deltaStar.size(), 0.);
- deltaStar = _deltaStar;
-}
-void Edge::setVtOld(std::vector<double> const _vtOld)
-{
- vtOld.resize(_vtOld.size(), 0.);
- vtOld = _vtOld;
-}
diff --git a/vii/src/DEdge.h b/vii/src/DEdge.h
deleted file mode 100644
index 2eadb17..0000000
--- a/vii/src/DEdge.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef DEDGE_H
-#define DEDGE_H
-
-#include "vii.h"
-#include <vector>
-#include <iostream>
-
-namespace vii
-{
-
-/**
- * @brief Inviscid quantities at the edge of the boundary layer of a BoundaryLayer region.
- * @author Paul Dechamps
- */
-class VII_API Edge
-{
-private:
- std::vector<double> Me; ///< Mach number.
- std::vector<double> vt; ///< Velocity.
- std::vector<double> rhoe; ///< Density.
- std::vector<double> deltaStar; ///< Displacement thickness.
- std::vector<double> vtOld; ///< Previous velocity.
-
-public:
- Edge();
- ~Edge(){};
-
- // Getters.
- double getMe(size_t iPoint) const { return Me[iPoint]; };
- double getVt(size_t iPoint) const { return vt[iPoint]; };
- double getRhoe(size_t iPoint) const { return rhoe[iPoint]; };
- double getDeltaStar(size_t iPoint) const { return deltaStar[iPoint]; };
- double getVtOld(size_t iPoint) const { return vtOld[iPoint]; };
- double getMaxMach() const;
-
- // Setters.
- void setMe(std::vector<double> const _Me);
- void setVt(std::vector<double> const _vt);
- void setRhoe(std::vector<double> const _rhoe);
- void setDeltaStar(std::vector<double> const _deltaStar);
- void setVtOld(std::vector<double> const _vtOld);
-};
-} // namespace vii
-#endif // DEDGE_H
diff --git a/vii/src/DFluxes.cpp b/vii/src/DFluxes.cpp
deleted file mode 100644
index 66ff922..0000000
--- a/vii/src/DFluxes.cpp
+++ /dev/null
@@ -1,242 +0,0 @@
-#include "DFluxes.h"
-#include "DBoundaryLayer.h"
-#include <Eigen/Dense>
-
-using namespace vii;
-using namespace Eigen;
-
-/**
- * @brief Compute residual vector of the integral boundary layer equations.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @return VectorXd
- */
-VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl)
-{
- size_t nVar = bl->getnVar();
- std::vector<double> u(nVar, 0.);
- for (size_t k = 0; k < nVar; ++k)
- u[k] = bl->u[iPoint * nVar + k];
- return blLaws(iPoint, bl, u);
-}
-
-/**
- * @brief Compute Jacobian of the integral boundary layer equations.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @param sysRes Residual of the IBL at point iPoint.
- * @param eps Pertubation from current solution used to compute the Jacobian.
- * @return MatrixXd
- */
-MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, VectorXd const &sysRes, double eps)
-{
- size_t nVar = bl->getnVar();
- MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
- std::vector<double> uUp(nVar, 0.);
- for (size_t k = 0; k < nVar; ++k)
- uUp[k] = bl->u[iPoint * nVar + k];
-
- double varSave = 0.;
- for (size_t iVar = 0; iVar < nVar; ++iVar)
- {
- varSave = uUp[iVar];
- uUp[iVar] += eps;
- JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps;
- uUp[iVar] = varSave;
- }
- return JacMatrix;
-}
-
-/**
- * @brief Integral boundary layer equation.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @param u Solution vector at point iPoint.
- * @return VectorXd
- */
-VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
-{
- size_t nVar = bl->getnVar();
-
- MatrixXd timeMatrix = MatrixXd::Zero(5, 5);
- VectorXd spaceVector = VectorXd::Zero(5);
-
- // Dissipation ratio.
- double dissipRatio;
- if (bl->name == "wake")
- dissipRatio = 0.9;
- else
- dissipRatio = 1.;
-
- bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta.
- bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness.
-
- // Boundary layer closure relations.
- if (bl->regime[iPoint] == 0)
- {
- // Laminar closure relations.
- std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
- bl->Hstar[iPoint] = lamParam[0];
- bl->Hstar2[iPoint] = lamParam[1];
- bl->Hk[iPoint] = lamParam[2];
- bl->cf[iPoint] = lamParam[3];
- bl->cd[iPoint] = lamParam[4];
- bl->delta[iPoint] = lamParam[5];
- bl->ctEq[iPoint] = lamParam[6];
- bl->us[iPoint] = lamParam[7];
- }
-
- else if (bl->regime[iPoint] == 1)
- {
- // Turbulent closure relations.
- std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
- bl->Hstar[iPoint] = turbParam[0];
- bl->Hstar2[iPoint] = turbParam[1];
- bl->Hk[iPoint] = turbParam[2];
- bl->cf[iPoint] = turbParam[3];
- bl->cd[iPoint] = turbParam[4];
- bl->delta[iPoint] = turbParam[5];
- bl->ctEq[iPoint] = turbParam[6];
- bl->us[iPoint] = turbParam[7];
- }
- else
- {
- std::cout << "Wrong regime\n"
- << std::endl;
- throw;
- }
-
- // Gradients computation.
- double dx = (bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1));
- double dxExt = (bl->mesh->getExt(iPoint) - bl->mesh->getExt(iPoint - 1));
- double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx;
- double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx;
- double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx;
- double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx;
- double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
-
- double dueExt_dx = (bl->blEdge->getVt(iPoint) - bl->blEdge->getVt(iPoint - 1)) / dxExt;
- double ddeltaStarExt_dx = (bl->blEdge->getDeltaStar(iPoint) - bl->blEdge->getDeltaStar(iPoint - 1)) / dxExt;
-
- double c = 4 / (M_PI * dx);
- double cExt = 4 / (M_PI * dxExt);
-
- // Amplification routine.
- double dN_dx = 0.0;
- double ax = 0.0;
- if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0)
- {
- // No forced transition and laminar regime.
- ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
- dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx;
- }
-
- double Mea = bl->blEdge->getMe(iPoint);
- double Hstara = bl->Hstar[iPoint];
- double Hstar2a = bl->Hstar2[iPoint];
- double Hka = bl->Hk[iPoint];
- double cfa = bl->cf[iPoint];
- double cda = bl->cd[iPoint];
- double deltaa = bl->delta[iPoint];
- double ctEqa = bl->ctEq[iPoint];
- double usa = bl->us[iPoint];
-
- // Space part.
- spaceVector(0) = dt_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2;
- spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * cda + Hstara * cfa / 2;
- spaceVector(2) = dN_dx - ax;
- spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
-
- if (bl->regime[iPoint] == 1)
- spaceVector(4) = (2 * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
-
- // Time part.
- timeMatrix(0, 0) = u[1] / u[3];
- timeMatrix(0, 1) = u[0] / u[3];
- timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
-
- timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
- timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
- timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
- timeMatrix(2, 2) = 1.;
-
- timeMatrix(3, 0) = -c * u[1];
- timeMatrix(3, 1) = -c * u[0];
- timeMatrix(3, 3) = 1.;
-
- if (bl->regime[iPoint] == 1)
- {
- timeMatrix(4, 3) = 2 * deltaa / (usa * u[3] * u[3]);
- timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * usa);
- }
- else
- timeMatrix(4, 4) = 1.0;
-
- return timeMatrix.partialPivLu().solve(spaceVector);
-}
-
-/**
- * @brief Amplification routines of the e^N method for transition capturing.
- *
- * @param Hk Kinematic shape parameter.
- * @param Ret Reynolds number based on the momentum thickness.
- * @param theta Momentum thickness.
- * @return double
- */
-double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
-{
- double dgr = 0.08;
- double Hk1 = 3.5;
- double Hk2 = 4.;
- double Hmi = (1 / (Hk - 1));
- double logRet_crit = 2.492 * std::pow(1 / (Hk - 1.), 0.43) + 0.7 * (std::tanh(14 * Hmi - 9.24) + 1.0); // value at which the amplification starts to grow
-
- double ax = 0.;
- if (Ret <= 0.)
- Ret = 1.;
- if (std::log10(Ret) < logRet_crit - dgr)
- ax = 0.;
- else
- {
- double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2 * dgr);
- double Rfac = 0.;
- if (rNorm <= 0)
- Rfac = 0.;
- if (rNorm >= 1)
- Rfac = 1.;
- else
- Rfac = 3 * rNorm * rNorm - 2 * rNorm * rNorm * rNorm;
-
- // Normal envelope amplification rate
- double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi + 0.1 * std::exp(-20 * Hmi);
- double arg = 3.87 * Hmi - 2.52;
- double dndRet = 0.028 * (Hk - 1) - 0.0345 * std::exp(-(arg * arg));
- ax = (m * dndRet / theta) * Rfac;
-
- // Set correction for separated profiles
- if (Hk > Hk1)
- {
- double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
- double Hfac = 0.;
- if (Hnorm >= 1)
- Hfac = 1.;
- else
- Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm;
- double ax1 = ax;
- double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5));
- double tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
- double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1), 1.5)) / theta;
- if (ax2 < 0.)
- ax2 = 0.;
- ax = Hfac * ax2 + (1 - Hfac) * ax1;
- if (ax < 0.)
- ax = 0.;
- }
- }
- return ax;
-}
\ No newline at end of file
diff --git a/vii/src/DFluxes.h b/vii/src/DFluxes.h
deleted file mode 100644
index 54f2920..0000000
--- a/vii/src/DFluxes.h
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef DFLUXES_H
-#define DFLUXES_H
-
-#include "vii.h"
-#include "DBoundaryLayer.h"
-#include <Eigen/Dense>
-
-namespace vii
-{
-
-/**
- * @brief Residual and Jacobian computation of the unsteady integral boundary layer equations.
- * @author Paul Dechamps
- */
-class VII_API Fluxes
-{
-public:
- double Re; ///< Freestream Reynolds number.
-
-public:
- Fluxes(double _Re) : Re(_Re){};
- ~Fluxes(){};
- Eigen::VectorXd computeResiduals(size_t iPoint, BoundaryLayer *bl);
- Eigen::MatrixXd computeJacobian(size_t iPoint, BoundaryLayer *bl, Eigen::VectorXd const &sysRes, double eps);
-
-private:
- Eigen::VectorXd blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u);
- double amplificationRoutine(double Hk, double Ret, double theta) const;
-};
-} // namespace vii
-#endif // DFLUXES_H
diff --git a/vii/src/DSolver.cpp b/vii/src/DSolver.cpp
deleted file mode 100644
index c08c278..0000000
--- a/vii/src/DSolver.cpp
+++ /dev/null
@@ -1,247 +0,0 @@
-#include "DSolver.h"
-#include "DBoundaryLayer.h"
-#include "DFluxes.h"
-#include <Eigen/Dense>
-
-using namespace Eigen;
-using namespace tbox;
-using namespace vii;
-
-/**
- * @brief Construct a new Time Solver:: Time Solver object.
- *
- * @param _CFL0 Initial CFL number.
- * @param _Minf Freestream Mach number.
- * @param Re Freestream Reynolds number.
- * @param _maxIt Maximum number of iterations.
- * @param _tol Convergence tolerance.
- * @param _itFrozenJac Number of iterations between which the Jacobian is frozen.
- */
-Solver::Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt, double _tol, unsigned int _itFrozenJac)
-{
- CFL0 = _CFL0;
- maxIt = _maxIt;
- tol = _tol;
- itFrozenJac0 = _itFrozenJac;
- Minf = std::max(_Minf, 0.1);
- initSln = true;
- sysMatrix = new Fluxes(Re);
-}
-
-/**
- * @brief Destroy the Time Solver:: Time Solver object.
- *
- */
-Solver::~Solver()
-{
- delete sysMatrix;
- std::cout << "~vii::Solver()\n";
-}
-
-/**
- * @brief Impose initial condition at one point.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- */
-void Solver::initialCondition(size_t iPoint, BoundaryLayer *bl)
-{
- size_t nVar = bl->getnVar();
- if (initSln)
- {
- bl->u[iPoint * nVar + 0] = bl->u[(iPoint - 1) * nVar + 0];
- bl->u[iPoint * nVar + 1] = bl->u[(iPoint - 1) * nVar + 1];
- bl->u[iPoint * nVar + 2] = bl->u[(iPoint - 1) * nVar + 2];
- bl->u[iPoint * nVar + 3] = bl->blEdge->getVt(iPoint);
- bl->u[iPoint * nVar + 4] = bl->u[(iPoint - 1) * nVar + 4];
-
- bl->Ret[iPoint] = bl->Ret[iPoint - 1];
- bl->cd[iPoint] = bl->cd[iPoint - 1];
- bl->cf[iPoint] = bl->cf[iPoint - 1];
- bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
- bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
- bl->Hk[iPoint] = bl->Hk[iPoint - 1];
- bl->ctEq[iPoint] = bl->ctEq[iPoint - 1];
- bl->us[iPoint] = bl->us[iPoint - 1];
- bl->delta[iPoint] = bl->delta[iPoint - 1];
- bl->deltaStar[iPoint] = bl->deltaStar[iPoint - 1];
- }
- if (bl->regime[iPoint] == 1 && bl->u[iPoint * nVar + 4] <= 0)
- bl->u[iPoint * nVar + 4] = 0.03;
-}
-
-/**
- * @brief Pseudo-time integration.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @return int
- */
-int Solver::integration(size_t iPoint, BoundaryLayer *bl)
-{
- size_t nVar = bl->getnVar();
-
- // Save initial condition.
- std::vector<double> sln0(nVar, 0.);
- for (size_t i = 0; i < nVar; ++i)
- sln0[i] = bl->u[iPoint * nVar + i];
-
- // Initialize time integration variables.
- double dx = bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1);
-
- // Initial time step.
- double CFL = CFL0;
- unsigned int itFrozenJac = itFrozenJac0;
- double timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
- MatrixXd IMatrix(5, 5);
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
-
- // Initial system residuals.
- VectorXd sysRes = sysMatrix->computeResiduals(iPoint, bl);
- double normSysRes0 = sysRes.norm();
- double normSysRes = normSysRes0;
-
- MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
- VectorXd slnIncr = VectorXd::Zero(5);
-
- unsigned int innerIters = 0; // Inner (non-linear) iterations
-
- while (innerIters < maxIt)
- {
- // Jacobian computation.
- if (innerIters % itFrozenJac == 0)
- JacMatrix = sysMatrix->computeJacobian(iPoint, bl, sysRes, 1e-8);
-
- // Linear system.
- slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-sysRes);
-
- // Solution increment.
- for (size_t k = 0; k < nVar; ++k)
- bl->u[iPoint * nVar + k] += slnIncr(k);
- bl->u[iPoint * nVar + 0] = std::max(bl->u[iPoint * nVar + 0], 1e-6);
- bl->u[iPoint * nVar + 1] = std::max(bl->u[iPoint * nVar + 1], 1.0005);
-
- sysRes = sysMatrix->computeResiduals(iPoint, bl);
- normSysRes = (sysRes + slnIncr / timeStep).norm();
-
- if (std::isnan(normSysRes))
- {
- resetSolution(iPoint, bl, sln0, true);
- sysRes = sysMatrix->computeResiduals(iPoint, bl);
- CFL = 1.0;
- timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
- normSysRes = (sysRes + slnIncr / timeStep).norm();
- itFrozenJac = 1;
- }
-
- if (normSysRes <= tol * normSysRes0)
- return 0; // Successfull exit.
-
- // CFL adaptation.
- CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
- timeStep = setTimeStep(CFL, dx, bl->blEdge->getVt(iPoint));
- IMatrix = MatrixXd::Identity(5, 5) / timeStep;
-
- innerIters++;
- }
-
- if (std::isnan(normSysRes) || normSysRes / normSysRes0 > 1e-3)
- {
- resetSolution(iPoint, bl, sln0);
- return -1;
- }
- return innerIters;
-}
-
-/**
- * @brief Set time step.
- *
- * @param CFL CFL number.
- * @param dx Cell size.
- * @param invVel Inviscid velocity.
- * @return double
- */
-double Solver::setTimeStep(double CFL, double dx, double invVel) const
-{
- // Speed of sound.
- double gradU2 = (invVel * invVel);
- double soundSpeed = computeSoundSpeed(gradU2);
-
- // Velocity of the fastest travelling wave.
- double advVel = soundSpeed + invVel;
-
- // Time step computation.
- return CFL * dx / advVel;
-}
-
-/**
- * @brief Compute the speed of sound in a cell.
- *
- * @param gradU2 Inviscid velocity squared in the cell.
- * @return double
- */
-double Solver::computeSoundSpeed(double gradU2) const
-{
- return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
-}
-
-/**
- * @brief Resets the solution of the point wrt its initial condition (sln0) or the previous point.
- *
- * @param iPoint Marker id.
- * @param bl Boundary layer region.
- * @param sln0 Initial solution at the point.
- * @param usePrevPoint if 1, use the solution at previous point instead of sln0.
- */
-void Solver::resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> sln0, bool usePrevPoint)
-{
- size_t nVar = bl->getnVar();
-
- // Reset solution.
- if (usePrevPoint)
- for (size_t k = 0; k < nVar; ++k)
- bl->u[iPoint * nVar + k] = bl->u[(iPoint - 1) * nVar + k];
- else
- for (size_t k = 0; k < nVar; ++k)
- bl->u[iPoint * nVar + k] = sln0[k];
-
- // Reset closures.
- bl->Ret[iPoint] = std::max(bl->u[iPoint * nVar + 3] * bl->u[iPoint * nVar + 0] * sysMatrix->Re, 1.0); // Reynolds number based on theta.
- bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1]; // Displacement thickness.
-
- if (bl->regime[iPoint] == 0)
- {
- std::vector<double> lamParam(8, 0.);
- bl->closSolver->laminarClosures(lamParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
- bl->Hstar[iPoint] = lamParam[0];
- bl->Hstar2[iPoint] = lamParam[1];
- bl->Hk[iPoint] = lamParam[2];
- bl->cf[iPoint] = lamParam[3];
- bl->cd[iPoint] = lamParam[4];
- bl->delta[iPoint] = lamParam[5];
- bl->ctEq[iPoint] = lamParam[6];
- bl->us[iPoint] = lamParam[7];
- }
- else if (bl->regime[iPoint] == 1)
- {
- std::vector<double> turbParam(8, 0.);
- bl->closSolver->turbulentClosures(turbParam, bl->u[iPoint * nVar + 0], bl->u[iPoint * nVar + 1], bl->u[iPoint * nVar + 4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
- bl->Hstar[iPoint] = turbParam[0];
- bl->Hstar2[iPoint] = turbParam[1];
- bl->Hk[iPoint] = turbParam[2];
- bl->cf[iPoint] = turbParam[3];
- bl->cd[iPoint] = turbParam[4];
- bl->delta[iPoint] = turbParam[5];
- bl->ctEq[iPoint] = turbParam[6];
- bl->us[iPoint] = turbParam[7];
- }
-}
-
-void Solver::setCFL0(double _CFL0)
-{
- if (_CFL0 > 0)
- CFL0 = _CFL0;
- else
- throw std::runtime_error("Negative CFL numbers are not allowed.\n");
-}
\ No newline at end of file
diff --git a/vii/src/DSolver.h b/vii/src/DSolver.h
deleted file mode 100644
index 435a7b7..0000000
--- a/vii/src/DSolver.h
+++ /dev/null
@@ -1,51 +0,0 @@
-#ifndef DSOLVER_H
-#define DSOLVER_H
-
-#include "vii.h"
-#include "DFluxes.h"
-
-namespace vii
-{
-
-/**
- * @brief Pseudo-time integration.
- * @author Paul Dechamps
- */
-class VII_API Solver
-{
-
-private:
- double CFL0; ///< Initial CFL number.
- unsigned int maxIt; ///< Maximum number of iterations.
- double tol; ///< Convergence tolerance.
- unsigned int itFrozenJac0; ///< Number of iterations between which the Jacobian is frozen.
- bool initSln; ///< Flag to initialize the solution at the point.
- double const gamma = 1.4; ///< Air heat capacity ratio.
- double Minf; ///< Freestream Mach number.
-
- Fluxes *sysMatrix;
-
-public:
- Solver(double _CFL0, double _Minf, double Re, unsigned int _maxIt = 100, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
- ~Solver();
- void initialCondition(size_t iPoint, BoundaryLayer *bl);
- int integration(size_t iPoint, BoundaryLayer *bl);
-
- // Getters.
- double getCFL0() const { return CFL0; };
- unsigned int getmaxIt() const { return maxIt; };
- unsigned int getitFrozenJac() const { return itFrozenJac0; };
-
- // Setters.
- void setCFL0(double _CFL0);
- void setmaxIt(unsigned int _maxIt) { maxIt = _maxIt; };
- void setitFrozenJac(unsigned int _itFrozenJac) { itFrozenJac0 = _itFrozenJac; };
- void setinitSln(bool _initSln) { initSln = _initSln; };
-
-private:
- double setTimeStep(double CFL, double dx, double invVel) const;
- double computeSoundSpeed(double gradU2) const;
- void resetSolution(size_t iPoint, BoundaryLayer *bl, std::vector<double> Sln0, bool usePrevPoint = false);
-};
-} // namespace vii
-#endif // DSOLVER_H
diff --git a/vii/src/vii.h b/vii/src/vii.h
deleted file mode 100644
index ac4d1e5..0000000
--- a/vii/src/vii.h
+++ /dev/null
@@ -1,53 +0,0 @@
-/*
- * Copyright 2020 University of Liège
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-// Global header of the "vii" module.
-
-#ifndef VII_H
-#define VII_H
-
-#if defined(WIN32)
-#ifdef vii_EXPORTS
-#define VII_API __declspec(dllexport)
-#else
-#define VII_API __declspec(dllimport)
-#endif
-#else
-#define VII_API
-#endif
-
-#include "tbox.h"
-
-/**
- * @brief Namespace for vii module
- */
-namespace vii
-{
-
-// Solvers
-class Driver;
-class Solver;
-
-// Data structures
-class BoundaryLayer;
-class Discretization;
-class Edge;
-
-// Other
-class Closures;
-} // namespace vii
-
-#endif // VII_H
diff --git a/vii/tests/bli2D.py b/vii/tests/bli2D.py
deleted file mode 100644
index 6d05345..0000000
--- a/vii/tests/bli2D.py
+++ /dev/null
@@ -1,126 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2022 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# @author Paul Dechamps
-# @date 2022
-# Test the vii implementation. The test case is a compressible attached transitional flow.
-# Tested functionalities;
-# - Time integration.
-# - Two time-marching techniques (agressive and safe).
-# - Transition routines.
-# - Compressible flow routines.
-# - Laminar and turbulent flow.
-# - Forced transition.
-#
-# Untested functionalities.
-# - Separation routines.
-# - Multiple failure case at one iteration.
-# - Response to unconverged inviscid solver.
-
-# Imports.
-import dart.utils as invUtils
-import dart.default as invDefault
-
-import vii.coupler as viiCoupler
-import vii.utils as viscUtils
-
-import tbox.utils as tboxU
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-import math
-
-def main():
- # Timer.
- tms = fwk.Timers()
- tms['total'].start()
-
- # Flow variables.
- Re = 1e7
- alpha = 2.*math.pi/180
- M_inf = 0.2
- xtrF = [-1, 0.40]
- # Numerical parameters.
- CFL0 = 1
-
- # Mesh the geometry.
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
-
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
- msh, gmshWritter = invDefault.mesh(2, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
- tms['msh'].stop()
-
- # Set the problem and add medium, initial/boundary conditions.
- tms['pre'].start()
- pbl, _, _, bnd, blw = invDefault.problem(msh, 2, alpha, 0., M_inf, 1.0, 1.0, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
- tms['pre'].stop()
-
- # Solve coupled problem.
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
-
- vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1, xtrF=xtrF)
- iSolverAPI = viscUtils.initDart(pbl, blw, vSolver)
- Coupler = viiCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=20, _couplTol=1e-4, _iterPrint=1, _resetInv=False)
- Coupler.run()
-
- # Save results.
- iSolverAPI.save(gmshWritter)
- Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolverAPI.solver.Cp)
- tboxU.write(Cp, 'Cp.dat', '%1.4e', ',', 'x, y, cp', '')
-
- print('')
- # Display results.
- print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' Re M alpha Cl Cd Cdp Cdf Cm')
- print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
-
- # Write results to file.
- vSolution = viscUtils.getSolution(vSolver)
- viscUtils.writeFile(vSolution, vSolver.getRe())
-
- # Display timers.
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
- print('SOLVERS statistics')
- print(Coupler.tms)
-
- # Plot results.
- invDefault.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(iSolverAPI.getCl(), vSolver.Cdt, iSolverAPI.getCm(), 4), 'invert': True})
- invDefault.plot(vSolution['x/c'], vSolution['Cf'], {'xlabel': '$x/c$', 'ylabel': '$c_f$', 'title': 'Cdt = {0:.{3}f}, Cdf = {1:.{3}f}, Cdp = {2:.{3}f}'.format(vSolver.Cdt, vSolver.Cdf, vSolver.Cdp, 4), 'xlim': [0, 1]})
-
- # Check results.
- # Test for n0012 (le=0.01, te=0.01), alpha=2, M=.2, Re=1e7.
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.235, 5e-3)) # XFOIL 0.2325
- tests.add(CTest('Cd', vSolution['Cdt'], 0.0051, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cdp', vSolution['Cdp'], 0.0005, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
- tests.add(CTest('xtr_top', vSolution['xtrT'], 0.21, 0.03, forceabs=True)) # XFOIL 0.1877
- tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/vii/utils.py b/vii/utils.py
deleted file mode 100644
index c090dc2..0000000
--- a/vii/utils.py
+++ /dev/null
@@ -1,149 +0,0 @@
-import vii
-from fwk.wutils import parseargs
-from fwk.coloring import ccolors
-import numpy as np
-
-def initBL(Re, Minf, CFL0, nSections, xtrF=[-1, -1], span=0, verb=None):
- """ Initialize boundary layer solver.
-
- Params
- ------
- - Re: Flow Reynolds number.
- - Minf: Freestream Mach number.
- - CFL0: Initial CFL number for time integration.
- - nSections: Number of sections in the domain.
- - xtrF: Forced transition location [upper, lower].
- - span: Wing span (not used for 2D computations.
- - verb: Verbosity level of the solver.
-
- Return
- ------
- - solver: vii::Driver class.
- """
- if Re<=0.:
- print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
- raise RuntimeError("Invalid parameter")
- if Minf<0.:
- print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Mach number.", Minf, ccolors.ANSI_RESET)
- raise RuntimeError("Invalid parameter")
- elif Minf>=1.:
- print(ccolors.ANSI_YELLOW, "dart::vUtils Warning : (Super)sonic freestream Mach number.", Minf, ccolors.ANSI_RESET)
- if nSections < 0:
- print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
- raise RuntimeError("Invalid parameter")
- if verb is None:
- args = parseargs()
- _verbose = args.verb
- else:
- if not(0<=verb<=3):
- print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", _verbose, ccolors.ANSI_RESET)
- raise RuntimeError("Invalid parameter")
- else:
- _verbose = verb
- _xtrF = [-1, -1]
- for i, xtr in enumerate(xtrF):
- if not(0 <= xtr <= 1) and xtr != -1:
- raise RuntimeError("Incorrect forced transition location.")
- _xtrF[i] = xtr
- solver = vii.Driver(Re, Minf, CFL0, nSections, xtrF[0], xtrF[1], span, verbose=_verbose)
- return solver
-
-def initDart(pbl, blw, vSolver, iters = 25):
- """Initialize Newton solver and create the interface object
- between dart and the viscous solver.
-
- Parameters
- ----------
- pbl (dart.Problem object): Problem formulation.
- blw (np.ndarray(dart.Blowing, dart.Blowing)): Blowing boundaries for (airfoil, wake).
- vSolver (vii.Driver object): Viscous solver.
-
- Return
- ------
- solverAPI (DartInterface object): Interface between Dart and the viscous solver.
- """
- import dart
- from vii.interfaces.dart.DartInterface import DartInterface as dInterface
- from tbox.solvers import LinearSolver
-
- args = parseargs()
- dartSolver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb, mIt=iters)
- solverAPI = dInterface(dartSolver, vSolver, blw)
- return solverAPI
-
-def mesh(file, pars):
- """Initialize mesh and mesh writer
-
- Parameters
- ----------
- file : str
- file contaning mesh (w/o extention)
- pars : dict
- parameters for mesh
- """
- import tbox.gmsh as tmsh
- # Load the mesh.
- msh = tmsh.MeshLoader(file,__file__).execute(**pars)
- return msh
-
-
-def getSolution(vSolver, iSec=0):
- """Extract viscous solution.
- Args
- ----
- - vSolver: vii::Driver class. Drive of the viscous calculations
- - iSec (int): Marker of the section (default: 0)
- """
- if iSec<0:
- raise RuntimeError("dart::viscU Invalid section number", iSec)
-
- solverOutput = vSolver.getSolution(iSec)
-
- sln = {'theta' : solverOutput[0],
- 'H' : solverOutput[1],
- 'N' : solverOutput[2],
- 'ue' : solverOutput[3],
- 'Ct' : solverOutput[4],
- 'deltaStar': solverOutput[5],
- 'Ret' : solverOutput[6],
- 'Cd' : solverOutput[7],
- 'Cf' : np.asarray(solverOutput[8])*np.asarray(solverOutput[3])**2,
- 'Hstar' : solverOutput[9],
- 'Hstar2' : solverOutput[10],
- 'Hk' : solverOutput[11],
- 'Cteq' : solverOutput[12],
- 'Us' : solverOutput[13],
- 'delta' : solverOutput[14],
- 'x/c' : solverOutput[15],
- 'xtrT' : vSolver.getxtr(iSec, 0),
- 'xtrB' : vSolver.getxtr(iSec, 1),
- 'Cdt' : vSolver.Cdt,
- 'Cdf' : vSolver.Cdf,
- 'Cdp' : vSolver.Cdp}
- return sln
-
-def writeFile(wData, Re, toW=['deltaStar', 'H', 'Cf', 'Ct', 'ue', 'delta']):
- """Write the results in airfoil_viscous.dat
- """
- # Write
- print('Writing file: airfoil_viscous.dat...', end = '')
- f = open('airfoil_viscous.dat', 'w+')
-
- f.write('$Aerodynamic coefficients\n')
- f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cdt'], wData['Cdp'],wData['Cdf'], wData['xtrT'], wData['xtrB']))
- f.write('$Boundary layer variables\n')
- f.write(' x/c')
-
- for s in toW:
- f.write(' {0:>15s}'.format(s))
- f.write('\n')
-
- for i in range(len(wData['x/c'])):
- f.write('{0:15.6f}'.format(wData['x/c'][i]))
- for s in toW:
- f.write(' {0:15.6f}'.format(wData[s][i]))
- f.write('\n')
-
- f.close()
- print('done.')
\ No newline at end of file
--
GitLab