diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h index 5f7c560002d6f9da2d0b186223c1edf1c035170e..54d3a31e17266bd06bf7acd18aac25e805049383 100644 --- a/dart/_src/dartw.h +++ b/dart/_src/dartw.h @@ -26,11 +26,3 @@ #include "wProblem.h" #include "wSolver.h" #include "wWake.h" -#include "wViscSolver.h" -#include "wBLRegion.h" -#include "wTimeSolver.h" -#include "wViscFluxes.h" -#include "wClosures.h" -#include "wInterpolator.h" -#include "wViscMesh.h" -#include "wInvBnd.h" diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i index f0583cfe0c0288c4a995913ce920b011cbcdd856..6847f36fbe3da3589edd5ee4e941da3ed3b3945b 100644 --- a/dart/_src/dartw.i +++ b/dart/_src/dartw.i @@ -58,15 +58,6 @@ threads="1" %shared_ptr(dart::Newton); %shared_ptr(dart::Adjoint); -%include "wViscSolver.h" -%include "wBLRegion.h" -%include "wTimeSolver.h" -%include "wViscFluxes.h" -%include "wClosures.h" -%include "wInterpolator.h" -%include "wViscMesh.h" -%include "wInvBnd.h" - %include "wFluid.h" %include "wAssign.h" %include "wFreestream.h" diff --git a/dart/src/dart.h b/dart/src/dart.h index f242f898daa4f09861f38f5b6d778a6d3b4ec5a2..ad5628cf805a2b095c93bf24e4a2f6224dc59c79 100644 --- a/dart/src/dart.h +++ b/dart/src/dart.h @@ -65,16 +65,6 @@ class Newton; class FpeLSFunction; class Adjoint; -// Viscous -class ViscSolver; -class BLRegion; -class TimeSolver; -class ViscSolver; -class Closures; -class Interpolator; -class ViscMesh; -class InvBnd; - // General class F0Ps; class F0El; diff --git a/dart/src/wInterpolator.cpp b/dart/src/wInterpolator.cpp deleted file mode 100644 index afef795e64dbf40310acbb880b08c92f1748a29c..0000000000000000000000000000000000000000 --- a/dart/src/wInterpolator.cpp +++ /dev/null @@ -1,18 +0,0 @@ -#include "wInterpolator.h" -#include "wClosures.h" -#include "wViscMesh.h" -#include "wInvBnd.h" -#include <iostream> -#include <vector> -#include <cmath> - -using namespace dart; - -Interpolator::Interpolator() -{ -} // end Interpolator - -Interpolator::~Interpolator() -{ -} // end ~Interpolator - diff --git a/dart/src/wInterpolator.h b/dart/src/wInterpolator.h deleted file mode 100644 index 95453bb05a790eb5097d294c2e1664d95137cc83..0000000000000000000000000000000000000000 --- a/dart/src/wInterpolator.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef WINTERPOLATOR_H -#define WINTERPOLATOR_H - -#include "dart.h" -#include "wClosures.h" -#include "wViscMesh.h" -#include <vector> -#include <string> - -namespace dart -{ - -class DART_API Interpolator -{ - ViscMesh *vmesh; - -public: - Interpolator(); - ~Interpolator(); -}; -} // namespace dart -#endif //WINTERPOLATOR_H \ No newline at end of file diff --git a/dart/tests/bli2.py b/dart/tests/bli2.py deleted file mode 100644 index e16ea38b4f02e135b1a46398fa90d0f0b776e5b4..0000000000000000000000000000000000000000 --- a/dart/tests/bli2.py +++ /dev/null @@ -1,158 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -#import dart.viscUtils as viscU - -import dart.pyVII.vCoupler as floVC -import dart.pyVII.vUtils as viscU -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors -from matplotlib import pyplot as plt - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 2*math.pi/180 - M_inf = 0.73 - - CFL0 = 1 - - nSections = 1 - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - # 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.005} - msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream']) - tms['msh'].stop() - - # set the problem and add medium, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - vsolver = viscU.initBL(Re, M_inf, CFL0, nSections) - config = {'nDim': dim, 'nSections':nSections, 'nPoints':[600, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10} - - coupler = floVC.Coupler(isolver, vsolver) - coupler.CreateProblem(blw[0], blw[1]) - coupler.Run() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # Write results to file - vSol = viscU.GetSolution(0, vsolver) - plt.plot(vSol['x/c'], vSol['Cf']) - plt.show() - viscU.WriteFile(vSol, 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) - - xtr=vsolver.Getxtr() - - """# visualize solution and plot results - #floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - blScalars, blVectors = viscU.ExtractVars(vsolver) - wData = viscU.Dictionary(blScalars, blVectors, xtr) - viscU.PlotVars(wData, plotVar)""" - - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True)) - if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True)) - else: - raise Exception('Test not defined for this flow') - - tests.run() - - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/bli_Polar.py b/dart/tests/bli_Polar.py deleted file mode 100644 index 3e08f1380f598bfee417e25bd366612cb4bfb3a0..0000000000000000000000000000000000000000 --- a/dart/tests/bli_Polar.py +++ /dev/null @@ -1,124 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -import dart.VII.viscUtils as viscU - -import dart.VII.vCoupler as floVC -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha=0 - alphaSeq = np.arange(-5,5.5,5) - M_inf = 0. - - meshFactor = 1 - CFL0 = 1 - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - # mesh the geometry - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.005, 'msLe' : 0.001} - msh, gmshWriter = floD.mesh(dim, '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 = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor) - - coupler = floVC.Coupler(isolver, vsolver) - - coupler.CreateProblem(blw[0], blw[1]) - - polar = coupler.RunPolar(alphaSeq) - #coupler.Run() - tms['solver'].stop() - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - print('SOLVERS statistics') - print(coupler.tms) - - #viscU.PlotPolar(polar) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - tests.add(CTest('Cl1', polar[1,0], -0.430901, 5e-3)) - tests.add(CTest('Cd1', polar[2,0], 0.005728, 1e-3, forceabs=True)) - tests.add(CTest('Cl2', polar[1,-1], 0.430797, 5e-3)) - tests.add(CTest('Cd2', polar[2,-1], 0.005789, 1e-3, forceabs=True)) - - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/bli_Sep.py b/dart/tests/bli_Sep.py deleted file mode 100644 index cd35876dc54c6ecb9af7e497c4167792e9728dc2..0000000000000000000000000000000000000000 --- a/dart/tests/bli_Sep.py +++ /dev/null @@ -1,149 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -import dart.VII.viscUtils as viscU - -import dart.VII.vCoupler as floVC -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors -import math - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 15*math.pi/180 - #alphaSeq = np.arange(-5, 5.5, 0.5) - M_inf = 0.2 - - meshFactor = 2 - CFL0 = 1 - - plotVar = 'H' - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - # 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.0001} - msh, gmshWriter = floD.mesh(dim, '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 = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor) - - coupler = floVC.Coupler(isolver, vsolver) - - coupler.CreateProblem(blw[0], blw[1]) - - #coupler.RunPolar(alphaSeq) - coupler.Run() - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - xtr=vsolver.Getxtr() - - - """# visualize solution and plot results - #floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - blScalars, blVectors = viscU.ExtractVars(vsolver) - xtr=vsolver.Getxtr() - wData = viscU.Dictionary(blScalars, blVectors, xtr) - viscU.PlotVars(wData, plotVar)""" - - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 1.61416, 5e-2)) - tests.add(CTest('Cd', vsolver.Cdt, 0.01452, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0105, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.01, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 1.00, 0.03, forceabs=True)) - else: - raise Exception('Test not defined for this flow') - - tests.run() - - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/bli_Trans.py b/dart/tests/bli_Trans.py deleted file mode 100644 index 9e27a0a12ea13501c48c16e1432b01297f9c490a..0000000000000000000000000000000000000000 --- a/dart/tests/bli_Trans.py +++ /dev/null @@ -1,148 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -import dart.VII.viscUtils as viscU - -import dart.VII.vCoupler as floVC -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 2.*math.pi/180 - #alphaSeq = np.arange(-5, 5.5, 0.5) - M_inf = 0.715 - - meshFactor = 2 - CFL0 = 1 - - plotVar = 'H' - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - # 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.005} - msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream']) - tms['msh'].stop() - - # set the problem and add medium, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor) - - coupler = floVC.Coupler(isolver, vsolver) - - coupler.CreateProblem(blw[0], blw[1]) - - #coupler.RunPolar(alphaSeq) - coupler.resetInv = True - coupler.Run() - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - - - """# visualize solution and plot results - #floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - blScalars, blVectors = viscU.ExtractVars(vsolver) - xtr=vsolver.Getxtr() - wData = viscU.Dictionary(blScalars, blVectors, xtr) - viscU.PlotVars(wData, plotVar)""" - - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58 - tests.add(CTest('Cd', vsolver.Cdt, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617 - tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176 - tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510 - tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473 - else: - raise Exception('Test not defined for this flow') - - tests.run() - - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py deleted file mode 100644 index 10c1ae6b5d204a67fb997400c3175ce4d84d16d4..0000000000000000000000000000000000000000 --- a/dart/tests/bli_lowLift.py +++ /dev/null @@ -1,195 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD -#import dart.viscUtils as viscU - -import dart.pyVII.vCoupler as floVC -import dart.pyVII.vUtils as viscU -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors -from matplotlib import pyplot as plt - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 5*math.pi/180 - M_inf = 0. - - CFL0 = 1 - - nSections = 1 - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - nRel = 1 - - nomDeVariableDePorc_x = [None for _ in range(nRel)] - nomDeVariableDePorc_var = [None for _ in range(nRel)] - nomDeVariableDePorc_Uin = [None for _ in range(nRel)] - - for i in range(nRel): - # 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.001} - msh, gmshWriter = floD.mesh(dim, '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 = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - vsolver = viscU.initBL(Re, M_inf, CFL0, nSections) - config = {'nDim': dim, 'nSections':nSections, 'nPoints':[600, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10} - - coupler = floVC.Coupler(isolver, vsolver) - coupler.CreateProblem(blw[0], blw[1]) - #fig, axs = plt.subplots(2) - aeroCoeff = coupler.Run() - #isolver.reset() - - nomDeVariableDePorc_x[i] = coupler.sln[15].copy() - nomDeVariableDePorc_var[i] = coupler.sln[1].copy() - nomDeVariableDePorc_Uin[i] = coupler.fMeshDict['vx'].copy() - - del isolver - del vsolver - del coupler - del pbl - del bnd - del blw - del msh - del gmshWriter - - print(len(nomDeVariableDePorc_x), len(nomDeVariableDePorc_var)) - for i in range(len(nomDeVariableDePorc_Uin)): - plt.plot(nomDeVariableDePorc_Uin[i] - nomDeVariableDePorc_Uin[0], 'x-') - plt.show() - for i in range(len(nomDeVariableDePorc_var)): - plt.plot(nomDeVariableDePorc_x[i], nomDeVariableDePorc_var[i] - nomDeVariableDePorc_var[0], 'x-') - plt.show() - - for i in range(len(nomDeVariableDePorc_var)): - plt.plot(nomDeVariableDePorc_x[i], nomDeVariableDePorc_var[i], '-') - plt.show() - - """sln = np.asarray(vsolver.GetSolution(0)) - - print(sln[8][sln[8]>0.1]) - plt.plot(sln[16], sln[1]) - plt.xlim([0, 2]) - plt.show()""" - #coupler.RunPolar(alphaSeq) - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - print('SOLVERS statistics') - print(coupler.tms) - - xtr=vsolver.Getxtr() - - """# visualize solution and plot results - #floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - blScalars, blVectors = viscU.ExtractVars(vsolver) - wData = viscU.Dictionary(blScalars, blVectors, xtr) - viscU.PlotVars(wData, plotVar)""" - - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True)) - if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True)) - else: - raise Exception('Test not defined for this flow') - - tests.run() - - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/bli_lowLift2.py b/dart/tests/bli_lowLift2.py deleted file mode 100644 index 05ab0abc10bb6735a221e0f5c48b6f569b2502d0..0000000000000000000000000000000000000000 --- a/dart/tests/bli_lowLift2.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. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -from email.errors import NonPrintableDefect -import dart.utils as floU -import dart.default as floD -#import dart.viscUtils as viscU - -import dart.pyVII.vCoupler2 as viscC # Rbf coupler -import dart.pyVII.vInterpolator as Interpol # Rbf interpolator -import dart.pyVII.vUtils as viscU -import numpy as np - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors -from matplotlib import pyplot as plt - -import math - -def main(): - print('Start') - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 5*math.pi/180 - M_inf = 0. - - CFL0 = 1 - - nSections = 1 - verbose = 0 - - # ======================================================================================== - - c_ref = 1 - dim = 2 - - nPv = [200] - aeroCoeff = [np.zeros(2) for i in range(len(nPv))] - - - for i in range(len(nPv)): - - # 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.001} - msh, sortElm, gmshWriter = floD.mesh(dim, '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 = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - - isolver = floD.newton(pbl) - if verbose == 0: - isolver.printOn = False - vsolver = viscU.initBL(Re, M_inf, CFL0, nSections, verbose) - config = {'nDim': dim, 'nSections':nSections, 'nPoints':[nPv[i], 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 15} - - iSolverAPI = Interpol.Interpolator(isolver, sortElm, blw, config) - - coupler = viscC.Coupler(iSolverAPI, vsolver) - aeroCoeff[i] = coupler.Run() - - """del isolver - del vsolver - del iSolverAPI - del pbl - del bnd - del blw - del msh - del gmshWriter""" - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - print('SOLVERS statistics') - print(coupler.tms) - - xtr=[vsolver.Getxtr(0,0),vsolver.Getxtr(0,1)] - - # visualize solution and plot results - #floD.initViewer(pbl) - #tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.2208, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.00531, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True)) - if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.5488, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.0062, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp,0.0018, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.06, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.74, 0.03, forceabs=True)) - else: - raise Exception('Test not defined for this flow') - - tests.run() - - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/blipython.py b/dart/tests/blipython.py deleted file mode 100755 index 8f40e54b6a12f47ae5f2a5e2b7138ece58af3091..0000000000000000000000000000000000000000 --- a/dart/tests/blipython.py +++ /dev/null @@ -1,157 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -# Copyright 2020 University of Liège -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - - -## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012 -# Amaury Bilocq -# -# Test the viscous-inviscid interaction scheme -# Reference to the master's thesis: http://hdl.handle.net/2268/252195 -# Reference test cases with Naca0012 (different from master's thesis): -# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741 -# -> msLe = 0.001, xtrTop = 0.061 -# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01 -# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740 -# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001 -# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1 -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import dart.utils as floU -import dart.default as floD - -import dart.viscous.Master.DCoupler as floVC -import dart.viscous.Drivers.DDriver as floVSMaster - -import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std -import dart.viscous.Solvers.Steady.StdSolver as floVS_Std - -import dart.viscous.Physics.DValidation as validation - -import tbox -import tbox.utils as tboxU -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -import math - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - Re = 1e7 - alpha = 12.*math.pi/180 - M_inf = 0. - - #user_xtr=[0.41,0.41] - #user_xtr=[0,None] - user_xtr=[None,None] - user_Ti=None - - mapMesh = 0 - nFactor = 2 - - # Time solver parameters - Time_Domain = False # True for unsteady solver, False for steady solver - CFL0 = 1 - - # ======================================================================================== - - U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1 - c_ref = 1 - dim = 2 - tol = 1e-4 # tolerance for the VII (usually one drag count) - case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case - - # 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.005} - #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001} - msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream']) - tms['msh'].stop() - - # set the problem and add medium, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc = True, dbc=True) - tms['pre'].stop() - - # solve viscous problem - - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - isolver = floD.newton(pbl) - - # Choose between steady and unsteady solver - if Time_Domain is True: # Run unsteady solver - vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor) - coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd) - - elif Time_Domain is False: # Run steady solver - vsolver=floVS_Std.Solver(Re) - coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter) - coupler.run() - tms['solver'].stop() - - # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - print('SOLVERS statistics') - print(coupler.tms) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if Re == 1e7 and M_inf == 0. and alpha == 2*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.2197, 5e-2)) # XFOIL 0.58 - tests.add(CTest('Cd', vsolver.Cd, 0.00529, 1e-3, forceabs=True)) # XFOIL 0.00617 - tests.add(CTest('Cdp', vsolver.Cdp, 0.0008, 1e-3, forceabs=True)) # XFOIL 0.00176 - tests.add(CTest('xtr_top', vsolver.xtr[0], 0.20395, 0.03, forceabs=True)) # XFOIL 0.0510 - tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.49996, 0.03, forceabs=True)) # XFOIL 0.7473 - else: - raise Exception('Test not defined for this flow') - - tests.run() - - # visualize solution and plot results - floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True) - - val=validation.Validation(case) - val.main(isolver.Cl,vsolver.wData) - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/vii/CMakeLists.txt b/vii/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..8de83175a84e209ca53d06ddf2b4a04a7702697f --- /dev/null +++ b/vii/CMakeLists.txt @@ -0,0 +1,28 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Add source dir +ADD_SUBDIRECTORY( src ) +ADD_SUBDIRECTORY( _src ) + +# Add test dir +MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests) + +# Add to install +INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py + ${CMAKE_CURRENT_LIST_DIR}/utils.py + DESTINATION vii) +INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api + ${CMAKE_CURRENT_LIST_DIR}/viscous + DESTINATION vii) diff --git a/vii/__init__.py b/vii/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..43be8dad7f573d13449b986919b7ba733479323e --- /dev/null +++ b/vii/__init__.py @@ -0,0 +1,21 @@ +# -*- coding: utf-8 -*- +# dart MODULE initialization file + +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +import fwk +import tbox +from viiw import * diff --git a/vii/_src/CMakeLists.txt b/vii/_src/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..18dc1584c0b273c5b28563a2eb1991565627aeb2 --- /dev/null +++ b/vii/_src/CMakeLists.txt @@ -0,0 +1,35 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# CMake input file of the SWIG wrapper around "viiw.so" + +INCLUDE(${SWIG_USE_FILE}) + +FILE(GLOB SRCS *.h *.cpp *.inl *.swg) +FILE(GLOB ISRCS *.i) + +SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON) +SET(CMAKE_SWIG_FLAGS "-interface" "_viiw") # avoids "import _module_d" with MSVC/Debug +SWIG_ADD_LIBRARY(viiw LANGUAGE python SOURCES ${ISRCS} ${SRCS}) +SET_PROPERTY(TARGET viiw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON) +MACRO_DebugPostfix(viiw) + +TARGET_INCLUDE_DIRECTORIES(viiw PRIVATE ${PROJECT_SOURCE_DIR}/vii/_src + ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src + ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src + ${PYTHON_INCLUDE_PATH}) +TARGET_LINK_LIBRARIES(viiw PRIVATE vii tbox fwk ${PYTHON_LIBRARIES}) + +INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/viiw.py DESTINATION ${CMAKE_INSTALL_PREFIX}) +INSTALL(TARGETS viiw DESTINATION ${CMAKE_INSTALL_PREFIX}) diff --git a/vii/_src/viiw.h b/vii/_src/viiw.h new file mode 100644 index 0000000000000000000000000000000000000000..806068743093274c2a5ba04e883db833fca3158f --- /dev/null +++ b/vii/_src/viiw.h @@ -0,0 +1,24 @@ +/* + * 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 "wViscSolver.h" +#include "wBLRegion.h" +#include "wTimeSolver.h" +#include "wViscFluxes.h" +#include "wClosures.h" +#include "wViscMesh.h" +#include "wInvBnd.h" +#include "wCoupler.h" \ No newline at end of file diff --git a/vii/_src/viiw.i b/vii/_src/viiw.i new file mode 100644 index 0000000000000000000000000000000000000000..2d788b18118ef354402b03998f6ca3fc4181d2ca --- /dev/null +++ b/vii/_src/viiw.i @@ -0,0 +1,53 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +// 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" + +// ----------- FLOW CLASSES ---------------- +%include "vii.h" + +%include "wViscSolver.h" +%include "wBLRegion.h" +%include "wTimeSolver.h" +%include "wViscFluxes.h" +%include "wClosures.h" +%include "wViscMesh.h" +%include "wInvBnd.h" +%include "wCoupler.h" \ No newline at end of file diff --git a/vii/pyVII/DAirfoil.py b/vii/pyVII/DAirfoil.py new file mode 100755 index 0000000000000000000000000000000000000000..e8986c86036c040ebac6860fbd7b67afbeacd81c --- /dev/null +++ b/vii/pyVII/DAirfoil.py @@ -0,0 +1,140 @@ +#!/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 vii.pyVII.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 + else: + idxStag = self.idxStag + #idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2)) + 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 + print(np.shape(connectListNodes)) + quit() + return connectListNodes, connectListElems,[uData, lData] diff --git a/vii/pyVII/DBoundary.py b/vii/pyVII/DBoundary.py new file mode 100755 index 0000000000000000000000000000000000000000..f3ff33aa19dcc02f214a3271010cd51f2490a284 --- /dev/null +++ b/vii/pyVII/DBoundary.py @@ -0,0 +1,36 @@ +#!/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/pyVII/DGroupConstructor.py b/vii/pyVII/DGroupConstructor.py new file mode 100755 index 0000000000000000000000000000000000000000..036d60481048bfa13871b2b838f4fb486aaace95 --- /dev/null +++ b/vii/pyVII/DGroupConstructor.py @@ -0,0 +1,180 @@ +from matplotlib import pyplot as plt +import numpy as np + +class GroupConstructor(): + """ ---------- Group Constructor ---------- + + Builds single vectors elements out of the 3 groups given in input. + """ + def __init__(self) -> None: + pass + + def mergeVectors(self,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 = self.__getBLcoordinates(dataIn[k][0]) + xx_lw, dv_lw, vtExt_lw, _, alpha_lw = self.__getBLcoordinates(dataIn[k][1]) + + else: # Wake case + + xx_wk, dv_wk, vtExt_wk, _, alpha_wk = self.__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 = self.createFineMesh(dataIn[0][0][:,0], nFactor) + fMeshLw = self.createFineMesh(dataIn[0][1][:,0], nFactor) + fMeshWk = self.createFineMesh(dataIn[1][:,0] , nFactor) + fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk)) + fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])] + + fxx_up = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor) + fxx_lw = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor) + fxx_wk = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor) + fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk)) + + fvtExt_up = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor) + fvtExt_lw = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor) + fvtExt_wk = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor) + fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk)) + + fMe_up = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor) + fMe_lw = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor) + fMe_wk = self.interpolateToFineMesh(dataIn[1][:,5] , fMeshWk, nFactor) + fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk)) + + frho_up = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor) + frho_lw = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor) + frho_wk = self.interpolateToFineMesh(dataIn[1][:,6] , fMeshWk, nFactor) + frho = np.concatenate((frho_up, frho_lw[1:], frho_wk)) + + fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor) + fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor) + fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7] , fMeshWk, nFactor) + + fxxExt_up = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor) + fxxExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor) + fxxExt_wk = self.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)) + + 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 = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor) + dataUpper = np.column_stack((dataUpper, interpParamUp)) + interpParamLw = self.interpolateToFineMesh(dataIn[0][1][:,iParam], fMeshLw, nFactor) + dataLower = np.column_stack((dataLower, interpParamLw)) + interpParamWk = self.interpolateToFineMesh(dataIn[1][:,iParam] , fMeshWk, nFactor) + dataWake = np.column_stack((dataWake, interpParamWk)) + + dataLower = np.delete(dataLower, 0, axis = 0) # Remove stagnation point doublon. + + 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])] + + 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])) + + 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(self,data): + """Transform the reference coordinates into airfoil coordinates + """ + nN = len(data[:,0]) + nE = nN-1 # Nbr of element along the surface. + vt = np.zeros(nN) # 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] + vt[i] = (data[i, 3] * np.cos(alpha[i]) + data[i ,4] * np.sin(alpha[i])) # Tangential speed at node 1 element i + vt[i+1] = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i])) # Tangential speed at node 2 element i + dv[i] = (vt[i+1] - vt[i]) / dx[i] # Velocity gradient according to element i. + # central scheme with half point + + return xx, dv, vt, nN, alpha + + def interpolateToFineMesh(self, cVector, fMesh, nFactor): + + 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(self, cMesh, nFactor): + + 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/pyVII/DWake.py b/vii/pyVII/DWake.py new file mode 100755 index 0000000000000000000000000000000000000000..d9f054032f68095b33ca3c8933b2bae38be0aee0 --- /dev/null +++ b/vii/pyVII/DWake.py @@ -0,0 +1,77 @@ +#!/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 + +from vii.pyVII.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] + 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/dart/pyVII/vCoupler.py b/vii/pyVII/vCoupler.py similarity index 89% rename from dart/pyVII/vCoupler.py rename to vii/pyVII/vCoupler.py index f74cd4c8427571fa3458a3f89b431d071b05087d..78d78c9237cb2fe0fd4530c789cf9d240b5a925f 100644 --- a/dart/pyVII/vCoupler.py +++ b/vii/pyVII/vCoupler.py @@ -22,9 +22,9 @@ import dart.utils as floU import tbox.utils as tboxU -from dart.viscous.Master.DAirfoil import Airfoil -from dart.viscous.Master.DWake import Wake -from dart.viscous.Drivers.DGroupConstructor import GroupConstructor +from vii.pyVII.DAirfoil import Airfoil +from vii.pyVII.DWake import Wake +from vii.pyVII.DGroupConstructor import GroupConstructor from fwk.coloring import ccolors import math @@ -35,15 +35,13 @@ import numpy as np from matplotlib import pyplot as plt class Coupler: - def __init__(self, iSolver, vSolver) -> None: + def __init__(self, iSolver, vSolver, _iterPrint=1) -> None: self.iSolver=iSolver self.vSolver=vSolver - self.vSolver.verbose = self.iSolver.verbose - #self.vSolver.verbose = 0 - self.maxCouplIter = 150 self.couplTol = 1e-4 + self.iterPrint = _iterPrint self.tms=fwk.Timers() self.resetInv = False @@ -61,18 +59,17 @@ class Coupler: # Initilize meshes in c++ objects - #self.iSolver.printOn = False - #self.vSolver.printOn = False - - dragIter = np.zeros((0,2)) + aeroCoeffs = np.zeros((0,2)) self.nElmUpperPrev = 0 couplIter = 0 cdPrev = 0.0 - """print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6' - .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/1e6)) - print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))""" + + print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6' + .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.GetRe()/1e6)) + print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error')) print(' --------------------------------------------------------') + while couplIter < self.maxCouplIter: # Run inviscid solver @@ -96,10 +93,6 @@ class Coupler: self.tms['viscous'].start() exitCode = self.vSolver.Run(couplIter) self.tms['viscous'].stop() - # print('Viscous ops terminated with exit code ', exitCode) - - #plt.xlim([0, 2]) - #plt.show() self.tms['processing'].start() self.__ImposeBlwVel(couplIter) @@ -107,35 +100,29 @@ class Coupler: error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1 - dragIter = np.vstack((dragIter, [self.iSolver.Cl, self.vSolver.Cdt])) + aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolver.Cl, self.vSolver.Cdt])) if error <= self.couplTol: print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error')) print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, - self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET) - return dragIter + self.iSolver.Cl, self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET) + return aeroCoeffs cdPrev = self.vSolver.Cdt - if couplIter%1==0: - xtr = [1, 1] - xtr[0]=self.vSolver.Getxtr(0,0) - xtr[1]=self.vSolver.Getxtr(0,1) + if couplIter%self.iterPrint==0: print(' {:>4.0f}| {:>10.4f} {:>10.4f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter, - self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error))) + self.iSolver.Cl, self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error))) + + if self.iSolver.verbose != 0 or self.vSolver.verbose != 0: + print('') couplIter += 1 self.tms['processing'].stop() - - self.sln = np.asarray(self.vSolver.GetSolution(0)) - """if couplIter==1: - #plt.plot(sln[16], sln[1], 'x-') - #plt.plot(self.fMeshDict['globCoord'][:self.fMeshDict['bNodes'][1]], self.fMeshDict['vx'][:self.fMeshDict['bNodes'][1]], 'x-') - return dragIter""" else: print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET) - return dragIter + return aeroCoeffs def RunPolar(self, alphaSeq): diff --git a/dart/pyVII/vCoupler2.py b/vii/pyVII/vCoupler2.py similarity index 100% rename from dart/pyVII/vCoupler2.py rename to vii/pyVII/vCoupler2.py diff --git a/dart/pyVII/vInterpolator.py b/vii/pyVII/vInterpolator.py similarity index 100% rename from dart/pyVII/vInterpolator.py rename to vii/pyVII/vInterpolator.py diff --git a/dart/pyVII/vRbf.py b/vii/pyVII/vRbf.py similarity index 100% rename from dart/pyVII/vRbf.py rename to vii/pyVII/vRbf.py diff --git a/dart/pyVII/vUtils.py b/vii/pyVII/vUtils.py similarity index 96% rename from dart/pyVII/vUtils.py rename to vii/pyVII/vUtils.py index 7009f8dca7918734920f917da062ba7a9e7427a5..16b8868e798b19e51ab5a448c5bd392c03cd4948 100644 --- a/dart/pyVII/vUtils.py +++ b/vii/pyVII/vUtils.py @@ -1,4 +1,5 @@ -import dart +import vii +from fwk.wutils import parseargs from fwk.coloring import ccolors def initBL(Re, Minf, CFL0, nSections, verbose=1): @@ -16,8 +17,9 @@ def initBL(Re, Minf, CFL0, nSections, verbose=1): if not(0<=verbose<=3): print(ccolors.ANSI_RED, "dart::vUtils Fatal error : verbose not in valid range.", verbose, ccolors.ANSI_RESET) raise RuntimeError("Invalid parameter") - - solver = dart.ViscSolver(Re, Minf, CFL0, nSections, verbose) + + args = parseargs() + solver = vii.ViscSolver(Re, Minf, CFL0, nSections, verbose=args.verb+1) return solver diff --git a/vii/src/CMakeLists.txt b/vii/src/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..1e17fcd45b8cc9c4cde0f806fb90db42ab6b26a0 --- /dev/null +++ b/vii/src/CMakeLists.txt @@ -0,0 +1,27 @@ +# Copyright 2020 University of Liège +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# CMake input file of dart.so + +FILE(GLOB SRCS *.h *.cpp *.inl *.hpp) + +ADD_LIBRARY(vii SHARED ${SRCS}) +MACRO_DebugPostfix(vii) +TARGET_INCLUDE_DIRECTORIES(vii PUBLIC ${PROJECT_SOURCE_DIR}/vii/src) + +TARGET_LINK_LIBRARIES(vii tbox) + +INSTALL(TARGETS vii DESTINATION ${CMAKE_INSTALL_PREFIX}) + +SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)") diff --git a/vii/src/vii.h b/vii/src/vii.h new file mode 100644 index 0000000000000000000000000000000000000000..74bd31de513b4ddf0e99eaf2d6582330699887ef --- /dev/null +++ b/vii/src/vii.h @@ -0,0 +1,50 @@ +/* + * 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 DART_API __declspec(dllexport) +#else +#define DART_API __declspec(dllimport) +#endif +#else +#define DART_API +#endif + +#include "tbox.h" + +/** + * @brief Namespace for dart module + */ +namespace vii +{ +// Viscous +class ViscSolver; +class BLRegion; +class TimeSolver; +class ViscSolver; +class Closures; +class Interpolator; +class ViscMesh; +class InvBnd; +} // namespace vii + +#endif //VII_H diff --git a/dart/src/wBLRegion.cpp b/vii/src/wBLRegion.cpp similarity index 99% rename from dart/src/wBLRegion.cpp rename to vii/src/wBLRegion.cpp index f87aadb8d4e0dce4519d78a82be4fa27e77159b5..0f83eccf4342444ca174461ca471b6967f776a3f 100644 --- a/dart/src/wBLRegion.cpp +++ b/vii/src/wBLRegion.cpp @@ -9,7 +9,7 @@ #include <iomanip> #include <deque> -using namespace dart; +using namespace vii; BLRegion::BLRegion(double _xtrF) { diff --git a/dart/src/wBLRegion.h b/vii/src/wBLRegion.h similarity index 98% rename from dart/src/wBLRegion.h rename to vii/src/wBLRegion.h index afd73d57469d8a477552492778cb5a899e3e027b..b109fdf7ce271c48f47458e8f1199b5f648b7a6d 100644 --- a/dart/src/wBLRegion.h +++ b/vii/src/wBLRegion.h @@ -1,14 +1,14 @@ #ifndef WBLREGION_H #define WBLREGION_H -#include "dart.h" +#include "vii.h" #include "wClosures.h" #include "wViscMesh.h" #include "wInvBnd.h" #include <vector> #include <string> -namespace dart +namespace vii { class DART_API BLRegion diff --git a/dart/src/wClosures.cpp b/vii/src/wClosures.cpp similarity index 99% rename from dart/src/wClosures.cpp rename to vii/src/wClosures.cpp index 3e299761ea9600f81e3f6b47e803b66acdacb7e3..10b2349952d107b8005be1a4273453f8f58a9a70 100644 --- a/dart/src/wClosures.cpp +++ b/vii/src/wClosures.cpp @@ -4,7 +4,7 @@ #include <string> using namespace tbox; -using namespace dart; +using namespace vii; // This is an empty constructor Closures::Closures() {} diff --git a/dart/src/wClosures.h b/vii/src/wClosures.h similarity index 88% rename from dart/src/wClosures.h rename to vii/src/wClosures.h index e86fc183f05d9abc9aee05b67934d17f123c17ad..23f92d9e1754406f5d238b4a69fc2e729d8a9617 100644 --- a/dart/src/wClosures.h +++ b/vii/src/wClosures.h @@ -1,11 +1,11 @@ #ifndef WCLOSURES_H #define WCLOSURES_H -#include "dart.h" +#include "vii.h" #include <vector> #include <string> -namespace dart +namespace vii { class DART_API Closures @@ -16,5 +16,5 @@ public: std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name) const; std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const; }; -} // namespace dart +} // namespace vii #endif //WClosures_H \ No newline at end of file diff --git a/vii/src/wDartCoupler.cpp b/vii/src/wDartCoupler.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cd40af47155f42ede0f8215ba6116422f3968025 --- /dev/null +++ b/vii/src/wDartCoupler.cpp @@ -0,0 +1,172 @@ +#include "wDartCoupler.h" +#include "../../dart/src/wNewton.h" +#include "../../dart/src/wBlowing.h" +#include "wViscSolver.h" + +using namespace vii; + +DartCoupler::DartCoupler(const Newton *_invSolver, ViscSolver *_viscSolver, const std::vector<*Blowing> _blw, unsigned long _nCouplIter = 150, double _tol = 1e-4) +{ + invSolver = _invSolver; + viscSolver = _viscSolver; + blw = _blw; + + maxCouplIter = _nCouplIter; + tol = _tol; + + Initialize(); +} +DartCoupler::~DartCoupler() + +void DartCoupler::Run(){ + double cdPrev = 0.; + double error = 1.; + + unsigned long nIters; + + while (nIters < maxCouplIter){ + + invSolver.run(); + + ImposeInviscidBC(); + + viscSolver.Run(); + + error = (viscSolver.Cdt - cdPrev) / viscSolver.Cdt; + + if (error < tol){ + break; + } + else{ + continue; + } + + cdPrev = viscSolver.Cdt; + + ImposeBlowing(); + + nIters++; + + } + + else{ + std::cout << "Maximum number of " << maxCouplIter << " reached." << std::endl; + } +} +void DartCoupler::ImposeInviscidBC(){ + + for (int iReg=0; iReg < blw.size(); ++iReg){ + + for (auto iNode : blw[iReg]->nodes){ + invSolver->U(iNode->row) + } + } +} + +void DartCoupler::ImposeBlowing() + +void DartCoupler::Initialize(){ + + /* Create connectivity list */ + + std::vector<std::vector<unsigned int>,2> connectListNodes; + std::vector<std::vector<unsigned int>,2> connectListElems; + + for (size_t iReg=0; iReg < blw.size(); ++iReg) + { + connectListNodes(iReg).resize(blw(iReg)->boundary->nodes.size(), 0); + connectListElems(iReg).resize(blw(iReg)->boudnary->tag->elems.size(), 0); + } + + + /* Find stagnation point */ + int idxStag = + + + 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 + else: + idxStag = self.idxStag + #idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2)) + 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/src/wDartCoupler.h b/vii/src/wDartCoupler.h new file mode 100644 index 0000000000000000000000000000000000000000..8910ecda39b957620d86509fe04c14d254ea6797 --- /dev/null +++ b/vii/src/wDartCoupler.h @@ -0,0 +1,35 @@ +#ifndef WDARTCOUPLER_H +#define WDARTCOUPLER_H + +#include "vii.h" +#include "../../dart/src/wNewton.h" +#include "../../dart/src/wBlowing.h" +#include "wViscSolver.h" + +namespace vii +{ + +class DART_API DartCoupler +{ + +private: + unsigned long maxCouplIter; + double tol; + + +public: + + Newton *invSolver; + ViscSolver *viscSolver; + std::vector<Blowing> *blw; + + + DartCoupler(); + ~DartCoupler(); + Run(); + ImposeBlowing(); + ImposeInviscidBC(); + +}; +} // namespace dart +#endif //WDARTCOUPLER_H diff --git a/dart/src/wInvBnd.cpp b/vii/src/wInvBnd.cpp similarity index 93% rename from dart/src/wInvBnd.cpp rename to vii/src/wInvBnd.cpp index d3ed14efc1670f0c92d7e709f7c6447ade862ccf..8d8722792b13f37e4afb4a033f2e8da8d05a60af 100644 --- a/dart/src/wInvBnd.cpp +++ b/vii/src/wInvBnd.cpp @@ -21,7 +21,7 @@ /* --- Namespaces -------------------------------------------------------------------- */ -using namespace dart; +using namespace vii; /* ----------------------------------------------------------------------------------- */ @@ -33,7 +33,7 @@ double InvBnd::GetMaxMach() const { if (Me.size() <= 0) { - throw std::runtime_error("dart::InvBnd Invalid Mach number vector."); + throw std::runtime_error("vii::InvBnd Invalid Mach number vector."); } double Mmax = 0.0; for (size_t i=0; i<Me.size(); ++i) diff --git a/dart/src/wInvBnd.h b/vii/src/wInvBnd.h similarity index 95% rename from dart/src/wInvBnd.h rename to vii/src/wInvBnd.h index d7f42990548a4a81baf65210df6ff6cbbb50e83b..b3f714cad2bc913f878b01546a2a04c430e26817 100644 --- a/dart/src/wInvBnd.h +++ b/vii/src/wInvBnd.h @@ -1,12 +1,12 @@ #ifndef WINVBND_H #define WINVBND_H -#include "dart.h" +#include "vii.h" #include "wBLRegion.h" #include <vector> #include <iostream> -namespace dart +namespace vii { class DART_API InvBnd @@ -40,5 +40,5 @@ class DART_API InvBnd double GetMaxMach() const; }; -} // namespace dart +} // namespace vii #endif //WINVBND_H \ No newline at end of file diff --git a/dart/src/wTimeSolver.cpp b/vii/src/wTimeSolver.cpp similarity index 99% rename from dart/src/wTimeSolver.cpp rename to vii/src/wTimeSolver.cpp index 3cae23f0cfc89cedcd83b8e01c83f8918129b8cb..dd420136a76b9353f2b23017068625233f1c60e5 100644 --- a/dart/src/wTimeSolver.cpp +++ b/vii/src/wTimeSolver.cpp @@ -9,7 +9,7 @@ using namespace Eigen; using namespace tbox; -using namespace dart; +using namespace vii; TimeSolver::TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac) { diff --git a/dart/src/wTimeSolver.h b/vii/src/wTimeSolver.h similarity index 97% rename from dart/src/wTimeSolver.h rename to vii/src/wTimeSolver.h index 6ab52df2eabf601d93afc2d45a0d1434dc23fe16..af485d4cdb6686206da9a04252711821ba256569 100644 --- a/dart/src/wTimeSolver.h +++ b/vii/src/wTimeSolver.h @@ -1,12 +1,12 @@ #ifndef WTIMESOLVER_H #define WTIMESOLVER_H -#include "dart.h" +#include "vii.h" #include "wViscFluxes.h" #include <vector> #include <string> -namespace dart +namespace vii { class DART_API TimeSolver diff --git a/dart/src/wViscFluxes.cpp b/vii/src/wViscFluxes.cpp similarity index 99% rename from dart/src/wViscFluxes.cpp rename to vii/src/wViscFluxes.cpp index c94070889a2582dca6a58cdb46425153161397c9..7fe7947fbfab3dd2f0b0ff734502c6c0b927dce6 100644 --- a/dart/src/wViscFluxes.cpp +++ b/vii/src/wViscFluxes.cpp @@ -9,7 +9,7 @@ #include <algorithm> #include <string> -using namespace dart; +using namespace vii; using namespace Eigen; ViscFluxes::ViscFluxes(double _Re) diff --git a/dart/src/wViscFluxes.h b/vii/src/wViscFluxes.h similarity index 91% rename from dart/src/wViscFluxes.h rename to vii/src/wViscFluxes.h index c295c9863e43f313175ed9ac03e135dea6ff0d49..a60bf35f34c0e2099b608e539694c2f4fd196a44 100644 --- a/dart/src/wViscFluxes.h +++ b/vii/src/wViscFluxes.h @@ -1,13 +1,13 @@ #ifndef WViscFluxes_H #define WViscFluxes_H -#include "dart.h" +#include "vii.h" #include "wBLRegion.h" #include <Eigen/Dense> #include <vector> #include <string> -namespace dart +namespace vii { class DART_API ViscFluxes @@ -25,5 +25,5 @@ private: Eigen::VectorXd BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u); double AmplificationRoutine(double hk, double ret, double theta) const; }; -} // namespace dart +} // namespace vii #endif //WViscFluxes_H \ No newline at end of file diff --git a/dart/src/wViscMesh.cpp b/vii/src/wViscMesh.cpp similarity index 95% rename from dart/src/wViscMesh.cpp rename to vii/src/wViscMesh.cpp index aafc15455d2111e3e7048df0553912fc2bf9f5e2..9cb054452bdf98b13c0f193d6493d7ab10d08d37 100644 --- a/dart/src/wViscMesh.cpp +++ b/vii/src/wViscMesh.cpp @@ -17,7 +17,7 @@ /* --- Namespaces -------------------------------------------------------------------- */ -using namespace dart; +using namespace vii; /* ----------------------------------------------------------------------------------- */ @@ -29,7 +29,7 @@ void ViscMesh::SetGlob(std::vector<double> _x, std::vector<double> _y, std::vect { if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size()) { - throw std::runtime_error("dart::ViscMesh Wrong mesh sizes."); + throw std::runtime_error("vii::ViscMesh Wrong mesh sizes."); } nMarkers = _x.size(); diff --git a/dart/src/wViscMesh.h b/vii/src/wViscMesh.h similarity index 97% rename from dart/src/wViscMesh.h rename to vii/src/wViscMesh.h index ff1a53e269d100461dcf173db863066fead0702b..0531097bc4bcbec7c300c2278860b0f9109134b7 100644 --- a/dart/src/wViscMesh.h +++ b/vii/src/wViscMesh.h @@ -1,12 +1,12 @@ #ifndef WVISCMESH_H #define WVISCMESH_H -#include "dart.h" +#include "vii.h" #include "wBLRegion.h" #include <vector> #include <iostream> -namespace dart +namespace vii { class DART_API ViscMesh diff --git a/dart/src/wViscSolver.cpp b/vii/src/wViscSolver.cpp similarity index 99% rename from dart/src/wViscSolver.cpp rename to vii/src/wViscSolver.cpp index 1dc6e70362b226f736b8e9b48c0277f3148a7336..e4270fed1af41494016b3686cacc5a1632cf4e54 100644 --- a/dart/src/wViscSolver.cpp +++ b/vii/src/wViscSolver.cpp @@ -13,7 +13,6 @@ #include "wViscSolver.h" #include "wBLRegion.h" #include "wTimeSolver.h" -#include "wInterpolator.h" #include <iostream> #include <iomanip> #include <vector> @@ -21,7 +20,7 @@ /* --- Namespaces -------------------------------------------------------------------- */ -using namespace dart; +using namespace vii; /* ----------------------------------------------------------------------------------- */ @@ -42,7 +41,6 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec verbose = _verbose; - /* Initialzing boundary layer */ Sections.resize(nSections); diff --git a/dart/src/wViscSolver.h b/vii/src/wViscSolver.h similarity index 98% rename from dart/src/wViscSolver.h rename to vii/src/wViscSolver.h index fdf87267c5a6333670e9989656dd9c3403cd0fb1..d5006506cfe021a6d90d52cc0098a577beb591ef 100644 --- a/dart/src/wViscSolver.h +++ b/vii/src/wViscSolver.h @@ -1,9 +1,9 @@ #ifndef WVISCSOLVER_H #define WVISCSOLVER_H -#include "dart.h" +#include "vii.h" #include "wBLRegion.h" -namespace dart +namespace vii { class DART_API ViscSolver diff --git a/dart/tests/bli_dualMesh.py b/vii/tests/bli2.py similarity index 62% rename from dart/tests/bli_dualMesh.py rename to vii/tests/bli2.py index ac96d71565c5808a2b2b82bee2a32e4419ed77b1..c46a10a8e9a8fc058a3c827bd315b07aff9bf728 100644 --- a/dart/tests/bli_dualMesh.py +++ b/vii/tests/bli2.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python3 +#!/usr/bin/env python3 # -*- coding: utf-8 -*- # Copyright 2020 University of Liège @@ -34,18 +34,18 @@ # This test is provided to ensure that the solver works properly. # Mesh refinement may have to be performed to obtain physical results. -import dart.utils as floU -import dart.default as floD -import dart.VII.viscUtils as viscU +import dart.utils as invUtils +import dart.default as invDefault -import dart.VII.vCoupler as floVC -import numpy as np +import vii.pyVII.vCoupler as viiCoupler +import vii.pyVII.vUtils as viscUtils import tbox import tbox.utils as tboxU import fwk from fwk.testing import * from fwk.coloring import ccolors +from matplotlib import pyplot as plt import math @@ -57,15 +57,11 @@ def main(): # define flow variables Re = 1e7 - alpha = 2*math.pi/180 - #alphaSeq = np.arange(-5, 5.5, 0.5) + alpha = 5.*math.pi/180 M_inf = 0. - meshFactor = 2 CFL0 = 1 - - plotVar = 'H' - + # ======================================================================================== c_ref = 1 @@ -74,13 +70,14 @@ def main(): # 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.001} - msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream']) + msh, gmshWriter = invDefault.mesh(dim, '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 = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) + pbl, _, _, bnd, blw = invDefault.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) tms['pre'].stop() # solve viscous problem @@ -88,53 +85,52 @@ def main(): print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) tms['solver'].start() - isolver = floD.newton(pbl) - vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor) + iSolver = invDefault.newton(pbl) + vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1) - coupler = floVC.Coupler(isolver, vsolver) - - coupler.CreateProblem(blw[0], blw[1]) - - #coupler.RunPolar(alphaSeq) - coupler.Run() - tms['solver'].stop() + Coupler = viiCoupler.Coupler(iSolver, vSolver) + Coupler.CreateProblem(blw[0], blw[1]) + Coupler.Run() # extract Cp - Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp) - tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '') + Cp = invUtils.extract(bnd.groups[0].tag.elems, iSolver.Cp) + tboxU.write(Cp, 'Cp.dat', '%1.5e', ', ', 'x, y, z, Cp', '') + + print('') # display results print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) print(' Re M alpha Cl Cd Cdp Cdf Cm') - print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm)) + print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, iSolver.Cl, vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolver.Cm)) + + # Write results to file + vSolution = viscUtils.GetSolution(0, vSolver) + viscUtils.WriteFile(vSolution, vSolver.GetRe()) + print('') # display timers tms['total'].stop() print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) print('CPU statistics') print(tms) - - xtr=vsolver.Getxtr() - - - """# visualize solution and plot results - #floD.initViewer(pbl) - tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True) - - blScalars, blVectors = viscU.ExtractVars(vsolver) - xtr=vsolver.Getxtr() - wData = viscU.Dictionary(blScalars, blVectors, xtr) - viscU.PlotVars(wData, plotVar)""" + print('SOLVERS statistics') + print(Coupler.tms) # check results print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) tests = CTests() if Re == 1e7 and M_inf == 0. and alpha == 2.*math.pi/180: - tests.add(CTest('Cl', isolver.Cl, 0.22113, 5e-3)) - tests.add(CTest('Cd', vsolver.Cdt, 0.00534, 1e-3, forceabs=True)) - tests.add(CTest('Cdp', vsolver.Cdp, 0.0009, 1e-3, forceabs=True)) - tests.add(CTest('xtr_top', xtr[0], 0.20, 0.03, forceabs=True)) - tests.add(CTest('xtr_bot', xtr[1], 0.50, 0.03, forceabs=True)) + tests.add(CTest('Cl', iSolver.Cl, 0.2208, 5e-3)) + tests.add(CTest('Cd', vSolution['Cdt'], 0.00531, 1e-3, forceabs=True)) + tests.add(CTest('Cdp', vSolution['Cdp'], 0.0009, 1e-3, forceabs=True)) + tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) + tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.50, 0.03, forceabs=True)) + if Re == 1e7 and M_inf == 0. and alpha == 5.*math.pi/180: + tests.add(CTest('Cl', iSolver.Cl, 0.554, 5e-3)) + tests.add(CTest('Cd', vSolution['Cdt'], 0.007, 1e-3, forceabs=True)) + tests.add(CTest('Cdp', vSolution['Cdp'],0.003, 1e-3, forceabs=True)) + tests.add(CTest('xtr_top', vSolution['xtrT'], 0.06, 0.03, forceabs=True)) + tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.74, 0.03, forceabs=True)) else: raise Exception('Test not defined for this flow') diff --git a/dart/tests/bli3.py b/vii/tests/bli3.py similarity index 100% rename from dart/tests/bli3.py rename to vii/tests/bli3.py