From a907d1298bc85c9f7c303fa80e452d5dec2fe183 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paulzer@Pauls-MacBook-Pro.local>
Date: Fri, 29 Apr 2022 14:08:26 +0200
Subject: [PATCH] Petite blague
---
dart/default.py | 5 -
dart/models/n0012_3.geo | 443 +++++++++++++++++++++++++++++++++++
dart/models/old/n0012-3_.geo | 5 +-
dart/pyVII/vCoupler.py | 2 +-
dart/pyVII/vCoupler2.py | 20 +-
dart/pyVII/vInterpolator.py | 349 +++++++++++++--------------
dart/pyVII/vStructures.py | 46 ++++
dart/pyVII/vUtils.py | 22 +-
dart/src/wTimeSolver.cpp | 1 +
dart/src/wViscFluxes.cpp | 72 +++---
dart/src/wViscSolver.cpp | 8 +-
dart/src/wViscSolver.h | 4 +-
dart/tests/bli3.py | 6 +-
dart/tests/bli_lowLift.py | 17 +-
dart/tests/bli_lowLift2.py | 91 ++++---
15 files changed, 799 insertions(+), 292 deletions(-)
create mode 100644 dart/models/n0012_3.geo
create mode 100644 dart/pyVII/vStructures.py
diff --git a/dart/default.py b/dart/default.py
index f7909a8..caf91dc 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,11 +23,6 @@ from fwk.wutils import parseargs
import dart
import tbox
-def initBL(Re, Minf, CFL0, meshFactor):
- solver = dart.ViscSolver(Re, Minf, CFL0, meshFactor)
- return solver
-
-
def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
"""Initialize mesh and mesh writer
diff --git a/dart/models/n0012_3.geo b/dart/models/n0012_3.geo
new file mode 100644
index 0000000..8ca6385
--- /dev/null
+++ b/dart/models/n0012_3.geo
@@ -0,0 +1,443 @@
+/* Rectangular NACA0012 wing */
+// Initially generated by unsRgridWingGen.m
+// For gmsh 4, use line 334 instead of line 335 (Bezier <- BSpline)
+
+// Parameters
+// domain and mesh
+DefineConstant[ spn = { 2.0, Name "Wing span" } ];
+DefineConstant[ lgt = { 3.0, Name "Channel length" } ];
+DefineConstant[ wdt = { 3.0, Name "Channel width" } ];
+DefineConstant[ hgt = { 3.0, Name "Channel height" } ];
+DefineConstant[ msLe = { 0.05, Name "Leading edge mesh size" } ];
+DefineConstant[ msTe = { 0.05, Name "Trailing edge mesh size" } ];
+DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ];
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1: naca0012, 101 points
+Point(1) = {1.000000,0.000000,0.000000,msTe};
+Point(2) = {0.999013,0.000000,0.000143};
+Point(3) = {0.996057,0.000000,0.000572};
+Point(4) = {0.991144,0.000000,0.001280};
+Point(5) = {0.984292,0.000000,0.002260};
+Point(6) = {0.975528,0.000000,0.003501};
+Point(7) = {0.964888,0.000000,0.004990};
+Point(8) = {0.952414,0.000000,0.006710};
+Point(9) = {0.938153,0.000000,0.008643};
+Point(10) = {0.922164,0.000000,0.010770};
+Point(11) = {0.904508,0.000000,0.013071,msTe};
+Point(12) = {0.885257,0.000000,0.015523};
+Point(13) = {0.864484,0.000000,0.018106};
+Point(14) = {0.842274,0.000000,0.020795};
+Point(15) = {0.818712,0.000000,0.023569};
+Point(16) = {0.793893,0.000000,0.026405};
+Point(17) = {0.767913,0.000000,0.029279};
+Point(18) = {0.740877,0.000000,0.032168};
+Point(19) = {0.712890,0.000000,0.035048};
+Point(20) = {0.684062,0.000000,0.037896};
+Point(21) = {0.654508,0.000000,0.040686};
+Point(22) = {0.624345,0.000000,0.043394};
+Point(23) = {0.593691,0.000000,0.045992};
+Point(24) = {0.562667,0.000000,0.048455};
+Point(25) = {0.531395,0.000000,0.050754};
+Point(26) = {0.500000,0.000000,0.052862};
+Point(27) = {0.468605,0.000000,0.054749};
+Point(28) = {0.437333,0.000000,0.056390};
+Point(29) = {0.406309,0.000000,0.057755};
+Point(30) = {0.375655,0.000000,0.058819};
+Point(31) = {0.345492,0.000000,0.059557};
+Point(32) = {0.315938,0.000000,0.059947};
+Point(33) = {0.287110,0.000000,0.059971,msTe};
+Point(34) = {0.259123,0.000000,0.059614};
+Point(35) = {0.232087,0.000000,0.058863};
+Point(36) = {0.206107,0.000000,0.057712};
+Point(37) = {0.181288,0.000000,0.056159};
+Point(38) = {0.157726,0.000000,0.054206};
+Point(39) = {0.135516,0.000000,0.051862};
+Point(40) = {0.114743,0.000000,0.049138};
+Point(41) = {0.095492,0.000000,0.046049};
+Point(42) = {0.077836,0.000000,0.042615};
+Point(43) = {0.061847,0.000000,0.038859};
+Point(44) = {0.047586,0.000000,0.034803};
+Point(45) = {0.035112,0.000000,0.030473};
+Point(46) = {0.024472,0.000000,0.025893};
+Point(47) = {0.015708,0.000000,0.021088};
+Point(48) = {0.008856,0.000000,0.016078};
+Point(49) = {0.003943,0.000000,0.010884};
+Point(50) = {0.000987,0.000000,0.005521};
+Point(51) = {0.000000,0.000000,0.000000,msLe};
+Point(52) = {0.000987,0.000000,-0.005521};
+Point(53) = {0.003943,0.000000,-0.010884};
+Point(54) = {0.008856,0.000000,-0.016078};
+Point(55) = {0.015708,0.000000,-0.021088};
+Point(56) = {0.024472,0.000000,-0.025893};
+Point(57) = {0.035112,0.000000,-0.030473};
+Point(58) = {0.047586,0.000000,-0.034803};
+Point(59) = {0.061847,0.000000,-0.038859};
+Point(60) = {0.077836,0.000000,-0.042615};
+Point(61) = {0.095492,0.000000,-0.046049};
+Point(62) = {0.114743,0.000000,-0.049138};
+Point(63) = {0.135516,0.000000,-0.051862};
+Point(64) = {0.157726,0.000000,-0.054206};
+Point(65) = {0.181288,0.000000,-0.056159};
+Point(66) = {0.206107,0.000000,-0.057712};
+Point(67) = {0.232087,0.000000,-0.058863};
+Point(68) = {0.259123,0.000000,-0.059614};
+Point(69) = {0.287110,0.000000,-0.059971,msTe};
+Point(70) = {0.315938,0.000000,-0.059947};
+Point(71) = {0.345492,0.000000,-0.059557};
+Point(72) = {0.375655,0.000000,-0.058819};
+Point(73) = {0.406309,0.000000,-0.057755};
+Point(74) = {0.437333,0.000000,-0.056390};
+Point(75) = {0.468605,0.000000,-0.054749};
+Point(76) = {0.500000,0.000000,-0.052862};
+Point(77) = {0.531395,0.000000,-0.050754};
+Point(78) = {0.562667,0.000000,-0.048455};
+Point(79) = {0.593691,0.000000,-0.045992};
+Point(80) = {0.624345,0.000000,-0.043394};
+Point(81) = {0.654508,0.000000,-0.040686};
+Point(82) = {0.684062,0.000000,-0.037896};
+Point(83) = {0.712890,0.000000,-0.035048};
+Point(84) = {0.740877,0.000000,-0.032168};
+Point(85) = {0.767913,0.000000,-0.029279};
+Point(86) = {0.793893,0.000000,-0.026405};
+Point(87) = {0.818712,0.000000,-0.023569};
+Point(88) = {0.842274,0.000000,-0.020795};
+Point(89) = {0.864484,0.000000,-0.018106};
+Point(90) = {0.885257,0.000000,-0.015523};
+Point(91) = {0.904508,0.000000,-0.013071,msTe};
+Point(92) = {0.922164,0.000000,-0.010770};
+Point(93) = {0.938153,0.000000,-0.008643};
+Point(94) = {0.952414,0.000000,-0.006710};
+Point(95) = {0.964888,0.000000,-0.004990};
+Point(96) = {0.975528,0.000000,-0.003501};
+Point(97) = {0.984292,0.000000,-0.002260};
+Point(98) = {0.991144,0.000000,-0.001280};
+Point(99) = {0.996057,0.000000,-0.000572};
+Point(100) = {0.999013,0.000000,-0.000143,msTe};
+
+// Airfoil 2: naca0012, 101 points
+Point(102) = {1.000000,spn,0.000000,msTe};
+Point(103) = {0.999013,spn,0.000143};
+Point(104) = {0.996057,spn,0.000572};
+Point(105) = {0.991144,spn,0.001280};
+Point(106) = {0.984292,spn,0.002260};
+Point(107) = {0.975528,spn,0.003501};
+Point(108) = {0.964888,spn,0.004990};
+Point(109) = {0.952414,spn,0.006710};
+Point(110) = {0.938153,spn,0.008643};
+Point(111) = {0.922164,spn,0.010770};
+Point(112) = {0.904508,spn,0.013071,msTe};
+Point(113) = {0.885257,spn,0.015523};
+Point(114) = {0.864484,spn,0.018106};
+Point(115) = {0.842274,spn,0.020795};
+Point(116) = {0.818712,spn,0.023569};
+Point(117) = {0.793893,spn,0.026405};
+Point(118) = {0.767913,spn,0.029279};
+Point(119) = {0.740877,spn,0.032168};
+Point(120) = {0.712890,spn,0.035048};
+Point(121) = {0.684062,spn,0.037896};
+Point(122) = {0.654508,spn,0.040686};
+Point(123) = {0.624345,spn,0.043394};
+Point(124) = {0.593691,spn,0.045992};
+Point(125) = {0.562667,spn,0.048455};
+Point(126) = {0.531395,spn,0.050754};
+Point(127) = {0.500000,spn,0.052862};
+Point(128) = {0.468605,spn,0.054749};
+Point(129) = {0.437333,spn,0.056390};
+Point(130) = {0.406309,spn,0.057755};
+Point(131) = {0.375655,spn,0.058819};
+Point(132) = {0.345492,spn,0.059557};
+Point(133) = {0.315938,spn,0.059947};
+Point(134) = {0.287110,spn,0.059971,msTe};
+Point(135) = {0.259123,spn,0.059614};
+Point(136) = {0.232087,spn,0.058863};
+Point(137) = {0.206107,spn,0.057712};
+Point(138) = {0.181288,spn,0.056159};
+Point(139) = {0.157726,spn,0.054206};
+Point(140) = {0.135516,spn,0.051862};
+Point(141) = {0.114743,spn,0.049138};
+Point(142) = {0.095492,spn,0.046049};
+Point(143) = {0.077836,spn,0.042615};
+Point(144) = {0.061847,spn,0.038859};
+Point(145) = {0.047586,spn,0.034803};
+Point(146) = {0.035112,spn,0.030473};
+Point(147) = {0.024472,spn,0.025893};
+Point(148) = {0.015708,spn,0.021088};
+Point(149) = {0.008856,spn,0.016078};
+Point(150) = {0.003943,spn,0.010884};
+Point(151) = {0.000987,spn,0.005521};
+Point(152) = {0.000000,spn,0.000000,msLe};
+Point(153) = {0.000987,spn,-0.005521};
+Point(154) = {0.003943,spn,-0.010884};
+Point(155) = {0.008856,spn,-0.016078};
+Point(156) = {0.015708,spn,-0.021088};
+Point(157) = {0.024472,spn,-0.025893};
+Point(158) = {0.035112,spn,-0.030473};
+Point(159) = {0.047586,spn,-0.034803};
+Point(160) = {0.061847,spn,-0.038859};
+Point(161) = {0.077836,spn,-0.042615};
+Point(162) = {0.095492,spn,-0.046049};
+Point(163) = {0.114743,spn,-0.049138};
+Point(164) = {0.135516,spn,-0.051862};
+Point(165) = {0.157726,spn,-0.054206};
+Point(166) = {0.181288,spn,-0.056159};
+Point(167) = {0.206107,spn,-0.057712};
+Point(168) = {0.232087,spn,-0.058863};
+Point(169) = {0.259123,spn,-0.059614};
+Point(170) = {0.287110,spn,-0.059971,msTe};
+Point(171) = {0.315938,spn,-0.059947};
+Point(172) = {0.345492,spn,-0.059557};
+Point(173) = {0.375655,spn,-0.058819};
+Point(174) = {0.406309,spn,-0.057755};
+Point(175) = {0.437333,spn,-0.056390};
+Point(176) = {0.468605,spn,-0.054749};
+Point(177) = {0.500000,spn,-0.052862};
+Point(178) = {0.531395,spn,-0.050754};
+Point(179) = {0.562667,spn,-0.048455};
+Point(180) = {0.593691,spn,-0.045992};
+Point(181) = {0.624345,spn,-0.043394};
+Point(182) = {0.654508,spn,-0.040686};
+Point(183) = {0.684062,spn,-0.037896};
+Point(184) = {0.712890,spn,-0.035048};
+Point(185) = {0.740877,spn,-0.032168};
+Point(186) = {0.767913,spn,-0.029279};
+Point(187) = {0.793893,spn,-0.026405};
+Point(188) = {0.818712,spn,-0.023569};
+Point(189) = {0.842274,spn,-0.020795};
+Point(190) = {0.864484,spn,-0.018106};
+Point(191) = {0.885257,spn,-0.015523};
+Point(192) = {0.904508,spn,-0.013071,msTe};
+Point(193) = {0.922164,spn,-0.010770};
+Point(194) = {0.938153,spn,-0.008643};
+Point(195) = {0.952414,spn,-0.006710};
+Point(196) = {0.964888,spn,-0.004990};
+Point(197) = {0.975528,spn,-0.003501};
+Point(198) = {0.984292,spn,-0.002260};
+Point(199) = {0.991144,spn,-0.001280};
+Point(200) = {0.996057,spn,-0.000572};
+Point(201) = {0.999013,spn,-0.000143,msTe};
+
+// Tip:
+Point(5101) = {0.999013,spn,0.000000};
+Point(5102) = {0.996057,spn+0.000125,0.000000};
+Point(5103) = {0.991144,spn+0.000331,0.000000};
+Point(5104) = {0.984292,spn+0.000620,0.000000};
+Point(5105) = {0.975528,spn+0.000989,0.000000};
+Point(5106) = {0.964888,spn+0.001437,0.000000};
+Point(5107) = {0.952414,spn+0.001963,0.000000};
+Point(5108) = {0.938153,spn+0.002563,0.000000};
+Point(5109) = {0.922164,spn+0.003237,0.000000};
+Point(5110) = {0.904508,spn+0.003981,0.000000,msTe};
+Point(5111) = {0.885257,spn+0.004791,0.000000};
+Point(5112) = {0.864484,spn+0.005666,0.000000};
+Point(5113) = {0.842274,spn+0.006602,0.000000};
+Point(5114) = {0.818712,spn+0.007594,0.000000};
+Point(5115) = {0.793893,spn+0.008640,0.000000};
+Point(5116) = {0.767913,spn+0.009734,0.000000};
+Point(5117) = {0.740877,spn+0.010873,0.000000};
+Point(5118) = {0.712890,spn+0.012052,0.000000};
+Point(5119) = {0.684062,spn+0.013266,0.000000};
+Point(5120) = {0.654508,spn+0.014511,0.000000};
+Point(5121) = {0.624345,spn+0.015781,0.000000};
+Point(5122) = {0.593691,spn+0.017072,0.000000};
+Point(5123) = {0.562667,spn+0.018379,0.000000};
+Point(5124) = {0.531395,spn+0.019696,0.000000};
+Point(5125) = {0.500000,spn+0.021019,0.000000};
+Point(5126) = {0.468605,spn+0.022341,0.000000};
+Point(5127) = {0.437333,spn+0.023658,0.000000};
+Point(5128) = {0.406309,spn+0.024965,0.000000};
+Point(5129) = {0.375655,spn+0.026256,0.000000};
+Point(5130) = {0.345492,spn+0.027526,0.000000};
+Point(5131) = {0.315938,spn+0.028771,0.000000};
+Point(5132) = {0.287110,spn+0.029985,0.000000,msTe};
+Point(5133) = {0.019492,spn+0.041801,0.000000};
+
+// Dummy tip center:
+Point(5349) = {0.904508,spn,0.000000};
+Point(5350) = {0.287110,spn,0.000000};
+
+// Box
+Point(10001) = {1+lgt, 0, 0, msF};
+Point(10002) = {1+lgt, 0, hgt, msF};
+Point(10003) = {-lgt, 0, hgt, msF};
+Point(10004) = {-lgt, 0, -hgt, msF};
+Point(10005) = {1+lgt, 0, -hgt, msF};
+Point(10011) = {1+lgt, wdt, hgt, msF};
+Point(10012) = {-lgt, wdt, hgt, msF};
+Point(10013) = {-lgt, wdt, -hgt, msF};
+Point(10014) = {1+lgt, wdt, -hgt, msF};
+
+// Wake
+Point(10021) = {1+lgt, spn, 0, msF};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1,2,3,4,5,6,7,8,9,10,11};
+Spline(2) = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33};
+Spline(3) = {33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51};
+Spline(4) = {51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69};
+Spline(5) = {69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91};
+Spline(6) = {91,92,93,94,95,96,97,98,99,100,1};
+
+// Airfoil 2:
+Spline(7) = {102,103,104,105,106,107,108,109,110,111,112};
+Spline(8) = {112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134};
+Spline(9) = {134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152};
+Spline(10) = {152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170};
+Spline(11) = {170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192};
+Spline(12) = {192,193,194,195,196,197,198,199,200,201,102};
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,102};
+Line(62) = {11,112};
+Line(63) = {33,134};
+Line(64) = {51,152};
+Line(65) = {69,170};
+Line(66) = {91,192};
+
+// Tip:
+Spline(121) = {102,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110};
+Spline(122) = {5110,5111,5112,5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132};
+If(GMSH_MAJOR_VERSION >= 4)
+ Bezier(123) = {5132,5133,152};
+Else
+ BSpline(123) = {5132,5133,152};
+EndIf
+Ellipse(124) = {112,5349,112,5110};
+Ellipse(125) = {134,5350,134,5132};
+Ellipse(126) = {170,5350,170,5132};
+Ellipse(127) = {192,5349,192,5110};
+
+// Box
+Line(201) = {10001,10002};
+Line(202) = {10002,10003};
+Line(203) = {10003,10004};
+Line(204) = {10004,10005};
+Line(205) = {10005,10001};
+Line(211) = {10011,10012};
+Line(212) = {10012,10013};
+Line(213) = {10013,10014};
+Line(214) = {10014,10011};
+Line(221) = {1, 10001};
+Line(222) = {102, 10021};
+Line(231) = {10002,10011};
+Line(232) = {10003,10012};
+Line(233) = {10004,10013};
+Line(234) = {10005,10014};
+
+// Wake
+Line(241) = {10001, 10021};
+
+/// Line loops & Surfaces
+
+// Wing 1:
+Line Loop(1) = {1,62,-7,-61};
+Line Loop(2) = {2,63,-8,-62};
+Line Loop(3) = {3,64,-9,-63};
+Line Loop(4) = {4,65,-10,-64};
+Line Loop(5) = {5,66,-11,-65};
+Line Loop(6) = {6,61,-12,-66};
+Surface(1) = {-1};
+Surface(2) = {-2};
+Surface(3) = {-3};
+Surface(4) = {-4};
+Surface(5) = {-5};
+Surface(6) = {-6};
+
+// Wingtip:
+Line Loop(11) = {7,124,-121};
+Line Loop(12) = {8,125,-122,-124};
+Line Loop(13) = {9,-123,-125};
+Line Loop(14) = {10,126,123};
+Line Loop(15) = {11,127,122,-126};
+Line Loop(16) = {12,121,-127};
+Surface(11) = {-11};
+Surface(12) = {-12};
+Surface(13) = {-13};
+Surface(14) = {-14};
+Surface(15) = {-15};
+Surface(16) = {-16};
+
+// Symmetry
+Line Loop(21) = {201, 202, 203, 204, 205};
+Line Loop(22) = {1, 2, 3, 4, 5, 6};
+Plane Surface(21) = {-21, 22};
+
+// Downsteam
+Line Loop(31) = {201, 231, -214, -234, 205};
+Plane Surface(31) = {31};
+
+// Farfield
+Line Loop(41) = {233, -212, -232, 203};
+Plane Surface(41) = {41};
+Line Loop(42) = {204, 234, -213, -233};
+Plane Surface(42) = {42};
+Line Loop(43) = {202, 232, -211, -231};
+Plane Surface(43) = {43};
+Line Loop(44) = {213, 214, 211, 212};
+Plane Surface(44) = {44};
+
+// Wake
+Line Loop(51) = {221, 241, -222, -61};
+Surface(51) = {51};
+
+/// Volume
+Surface Loop(1) = {1,2,3,4,5,6,11,12,13,14,15,16,21,31,41,42,43,44};
+Volume(1) = {1};
+
+/// Embbeded
+Line{221} In Surface{21};
+Line{241} In Surface{31};
+Surface{51} In Volume{1};
+
+//// MESHING ALGORITHM
+
+/// 2D:
+
+///Wing, farfield and symmetry plane:
+MeshAlgorithm Surface{1,2,3,4,5,6,21,31,41,42,43,44,51} = 5;
+
+///Tip:
+MeshAlgorithm Surface{11,12,13,14,15,16} = 1;
+
+/// 3D:
+
+Mesh.Algorithm3D = 2;
+//Mesh.OptimizeNetgen = 1;
+Mesh.Optimize = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {222};
+Physical Line("te") = {61};
+Physical Line("teTip") = {61, 222};
+
+// Internal Field:
+Physical Volume("field") = {1};
+
+// Wing:
+Physical Surface("wing") = {1,2,3,11,12,13};
+Physical Surface("wing_") = {4,5,6,14,15,16};
+
+// Symmetry:
+Physical Surface("symmetry") = {21};
+
+// Downstream:
+Physical Surface("downstream") = {31};
+
+// Farfield:
+Physical Surface("upstream") = {41};
+Physical Surface("farfield") = {42,43,44};
+
+// Wake:
+Physical Surface("wake") = {51};
diff --git a/dart/models/old/n0012-3_.geo b/dart/models/old/n0012-3_.geo
index 0bba944..1c31dbe 100644
--- a/dart/models/old/n0012-3_.geo
+++ b/dart/models/old/n0012-3_.geo
@@ -326,8 +326,9 @@ Mesh.SmoothNormals = 1;
// Physical
Physical Line("te") = {5};
-Physical Surface("cylinder") = {1};
-Physical Surface("cylinder_") = {2};
+//Physical Surface("cylinder") = {1};
+//Physical Surface("cylinder_") = {2};
+Physical Surface("wing") = {1,2};
Physical Surface("symmetry") = {3, 9};
Physical Surface("symmetry_") = {4, 10};
Physical Surface("upstream") = {5, 6};
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index fb95ffc..cde32fe 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -39,7 +39,7 @@ class Coupler:
self.vSolver.verbose = self.iSolver.verbose
#self.vSolver.verbose = 0
- self.maxCouplIter = 150
+ self.maxCouplIter = 1
self.couplTol = 1e-4
self.tms=fwk.Timers()
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
index fe8baee..1d86205 100644
--- a/dart/pyVII/vCoupler2.py
+++ b/dart/pyVII/vCoupler2.py
@@ -29,14 +29,14 @@ class Coupler:
def __init__(self, _iSolverAPI, vSolver) -> None:
self.iSolverAPI = _iSolverAPI
self.vSolver = vSolver
-
- self.vSolver.verbose = self.iSolverAPI.iSolver.verbose
self.maxCouplIter = 150
self.couplTol = 1e-4
self.tms=fwk.Timers()
self.resetInv = False
+
+ self.iSolverAPI.iSolver.printOn=False
pass
def Run(self):
@@ -44,7 +44,7 @@ class Coupler:
# Initilize meshes in c++ objects
- dragIter = np.zeros((0,2))
+ coeffConv = np.zeros((0,2))
self.nElmUpperPrev = 0
couplIter = 0
@@ -71,13 +71,12 @@ class Coupler:
self.tms['viscous'].start()
exitCode = self.vSolver.Run(couplIter)
self.tms['viscous'].stop()
- # print('Viscous ops terminated with exit code ', exitCode)
self.tms['processing'].start()
- self.iSolverAPI.SetBlowingVel(self.vSolver)
+ self.iSolverAPI.SetBlowingVel(self.vSolver, couplIter)
self.tms['processing'].stop()
- dragIter = np.vstack((dragIter, [self.iSolverAPI.GetCl(), self.vSolver.Cdt]))
+ coeffConv = np.vstack((coeffConv, [self.iSolverAPI.GetCl(), self.vSolver.Cdt]))
error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
if error <= self.couplTol:
@@ -86,7 +85,7 @@ class Coupler:
self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
self.tms['tot'].stop()
- return dragIter
+ return coeffConv
cdPrev = self.vSolver.Cdt
@@ -95,20 +94,15 @@ class Coupler:
self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
couplIter += 1
- """if couplIter == 5:
- quit()"""
- self.tms['processing'].stop()
else:
print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
self.tms['tot'].stop()
- return dragIter
+ return coeffConv
def RunPolar(self, alphaSeq):
- import math
-
self.iSolver.printOn = False
self.vSolver.printOn = False
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index f4b3948..1304404 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -3,10 +3,13 @@ from matplotlib import pyplot as plt
from fwk.coloring import ccolors
import dart.pyVII.vUtils as blxU
+import dart.pyVII.vStructures as vStruct
# 3D spline with Bézier curves
from scipy.interpolate import bisplrep
from scipy.interpolate import bisplev
+from scipy.interpolate import splrep
+from scipy.interpolate import splev
from scipy.interpolate import interp1d
@@ -15,21 +18,26 @@ from scipy.interpolate import RBFInterpolator
import fwk
-
class Interpolator:
def __init__(self, _dartSolver, bnd_wing, bnd_wake, _cfg):
self.iSolver = _dartSolver
+ self.invStruct = vStruct.invStruct()
+ self.viscStruct = vStruct.viscStruct()
+
self.bndWing = bnd_wing
self.bndWake = bnd_wake
- self.nodesList, self.transformedCoord, self.wakeCoord, self.viscStruct, self.viscStruct_tr = self.__GetWing(bnd_wing, bnd_wake, _cfg)
- print('Sections created')
self.config = _cfg
+ if self.config['nDim'] == 2:
+ self.nodesList, self.transformedCoord, self.wakeCoord, self.viscStruct, self.viscStruct_tr = self.__GetAirfoil(bnd_wing, bnd_wake, _cfg)
+ elif self.config['nDim'] == 3:
+
+
+ print('Sections created')
- self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(1)]
+ self.xxPrev = [[None for iReg in range(3)] for iSec in range(1)]
#self.uePrev = [[np.zeros(len(self.viscStruct[iSec][iReg][:,0])) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
- self.DeltaStarPrev = [[np.zeros(0
- ) for iReg in range(3)] for iSec in range(1)]
+ self.DeltaStarPrev = [[None for iReg in range(3)] for iSec in range(1)]
self.tms = fwk.Timers()
def GetInvBC(self, vSolver, couplIter):
@@ -41,8 +49,9 @@ class Interpolator:
# Wing
- nElm = len(self.nodesList)
- U_wg, M_wg, Rho_wg = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+ U_wg = np.zeros((len(self.nodesList), 3))
+ M_wg = np.zeros(len(self.nodesList))
+ Rho_wg = np.zeros(len(self.nodesList))
for iNode, n in enumerate(self.nodesList):
U_wg[iNode,0] = self.iSolver.U[n.row][0]
U_wg[iNode,1] = self.iSolver.U[n.row][1]
@@ -52,10 +61,9 @@ class Interpolator:
# Wake
- nElm = len(self.bndWake.nodes)
- U_wk =np.zeros((nElm, 3))
- M_wk = np.zeros(nElm)
- Rho_wk = np.zeros(nElm)
+ U_wk =np.zeros((len(self.bndWake.nodes), 3))
+ M_wk = np.zeros(len(self.bndWake.nodes))
+ Rho_wk = np.zeros(len(self.bndWake.nodes))
for iNode, n in enumerate(self.bndWake.nodes):
U_wk[iNode,0] = self.iSolver.U[n.row][0]
@@ -67,13 +75,14 @@ class Interpolator:
if self.config['nDim'] == 3: # x and y inverted in 3d
U_wg[:, [1,2]] = U_wg[:,[2,1]]
U_wk[:, [1,2]] = U_wk[:,[2,1]]
+
+ Ux = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,0], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,0], self.viscStruct[0][1][:,:1])]
+ Uy = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,1], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,1], self.viscStruct[0][1][:,:1])]
+ Uz = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,2], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,2], self.viscStruct[0][1][:,:1])]
- Ux = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,0], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,0], self.viscStruct[0][1][:,:1], 1)]
- Uy = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,1], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,1], self.viscStruct[0][1][:,:1], 1)]
- Uz = [self.__RbfToSections(self.transformedCoord[:,:2], U_wg[:,2], self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], U_wk[:,2], self.viscStruct[0][1][:,:1], 1)]
+ Me = [self.__RbfToSections(self.transformedCoord[:,:2], M_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], M_wk, self.viscStruct[0][1][:,:1])]
- Me = [self.__RbfToSections(self.transformedCoord[:,:2], M_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], M_wk, self.viscStruct[0][1][:,:1], 1)]
- Rho = [self.__RbfToSections(self.transformedCoord[:,:2], Rho_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], Rho_wk, self.viscStruct[0][1][:,:1], 1)]
+ Rho = [self.__RbfToSections(self.transformedCoord[:,:2], Rho_wg, self.viscStruct_tr[0][0][:,:2]), self.__RbfToSections(self.wakeCoord[:, :1], Rho_wk, self.viscStruct[0][1][:,:1])]
# Find stagnation point
for iSec in range(len(self.viscStruct)):
@@ -103,6 +112,11 @@ class Interpolator:
UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt)
UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt)
+ if couplIter == 0:
+ plt.plot(self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)[0], UxUp)
+ plt.plot(self.__Remesh(self.viscStruct[iSec][0][:,0], self.stgPt)[1], UxLw)
+ plt.show()
+
vSolver.SetVelocities(iSec, 0, UxUp, UyUp, UzUp)
vSolver.SetVelocities(iSec, 1, UxLw, UyLw, UzLw)
vSolver.SetVelocities(iSec, 2, Ux[1], Uy[1], Uz[1])
@@ -119,12 +133,10 @@ class Interpolator:
vSolver.SetRhoe(iSec, 1, RhoLw)
vSolver.SetRhoe(iSec, 2, Rho[1])
- #DeltaStarPrevUp, DeltaStarPrevLw = self.__Remesh(self.DeltaStarPrev[iSec][0], self.stgPt)
vSolver.SetDeltaStarExt(iSec, 0, self.DeltaStarPrev[0][0])
vSolver.SetDeltaStarExt(iSec, 1, self.DeltaStarPrev[0][1])
vSolver.SetDeltaStarExt(iSec, 2, self.DeltaStarPrev[0][2])
- #xxPrevUp, xxPrevLw = self.__Remesh(self.xxPrev[iSec][0], self.stgPt)
vSolver.SetxxExt(iSec, 0, self.xxPrev[0][0])
vSolver.SetxxExt(iSec, 1, self.xxPrev[0][1])
vSolver.SetxxExt(iSec, 2, self.xxPrev[0][2])
@@ -136,7 +148,7 @@ class Interpolator:
- def SetBlowingVel(self, vSolver):
+ def SetBlowingVel(self, vSolver, couplIter):
BlowingVel = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
for iSec in range(len(self.viscStruct)):
for iRegion in range(3):
@@ -147,12 +159,6 @@ class Interpolator:
if self.config['nDim'] == 2:
BlwWing = np.concatenate((BlowingVel[0][0][::-1], BlowingVel[0][1]))
BlwWake = BlowingVel[0][2]
- """self.DeltaStarPrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractDeltaStar(iSec, 0))[::-1], np.asarray(vSolver.ExtractDeltaStar(iSec, 1))[1:]), axis=None)
- self.DeltaStarPrev[iSec][1] = np.asarray(vSolver.ExtractDeltaStar(iSec, 2))
- self.xxPrev[iSec][0] = np.concatenate((np.asarray(vSolver.Extractxx(iSec, 0))[::-1], np.asarray(vSolver.Extractxx(iSec, 1))[1:]), axis=None)
- self.xxPrev[iSec][1] = np.asarray(vSolver.Extractxx(iSec, 2))
- self.uePrev[iSec][0] = np.concatenate((np.asarray(vSolver.ExtractUe(iSec, 0))[::-1], np.asarray(vSolver.ExtractUe(iSec, 1))[1:]), axis=None)
- self.uePrev[iSec][1] = np.asarray(vSolver.ExtractUe(iSec, 2))"""
meany = (self.viscStruct[0][0][np.argmin(self.viscStruct[0][0][:,0]), 1] - self.viscStruct[0][0][np.argmax(self.viscStruct[0][0][:,0]), 1]) / 2
@@ -200,9 +206,6 @@ class Interpolator:
mappedBlw = [self.__RbfToSections(viscCgWing[:, :2], BlwWing, invCgWing[:,:2]), self.__RbfToSections(viscCgWake[:, :1], BlwWake, invCgWake[:, :1])]
- #mappedBlwWing = RBFInterpolator(viscCgWing[:,:2], BlwWing, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(invCgWing[:,:2])
- #mappedBlwWake = RBFInterpolator(viscCgWake[:, :2], BlwWake, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(invCgWake[:, :2])
-
sln = blxU.GetSolution(0, vSolver)
self.tms['Interpolation'].stop()
@@ -222,7 +225,6 @@ class Interpolator:
def GetCp(self,bnd):
import dart.utils as floU
-
Cp = floU.extract(bnd, self.iSolver.Cp)
return Cp
@@ -245,20 +247,20 @@ class Interpolator:
xLw = vec[idx_stag:]
return xUp, xLw
- def __RbfToSections(self, x, var, xs, wake=0):
-
- if wake==0:
- return RBFInterpolator(x, var, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(xs)
+ def __RbfToSections(self, x, var, xs):
+ if np.all(var == var[0]):
+ varOut = RBFInterpolator(x, var, neighbors=1, kernel='linear', smoothing=0.0, degree=0)(xs)
else:
- return RBFInterpolator(x, var, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(xs)
-
+ varOut = RBFInterpolator(x, var, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(xs)
+ return varOut
+
def __GetWing(self, bnd_wing, bnd_wake, config):
if config['nDim'] == 2 and config['nSections'] != 1:
raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
- viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+ viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(1)] for _ in range(config['nSections'])]
# Create node list
nodesCoord = np.zeros((0, 3))
@@ -268,18 +270,135 @@ class Interpolator:
if config['nDim'] == 2:
spanDistri = np.zeros(1)
spanDistri[0] = nodesCoord[0,2] # Span direction is the z direction
+
+
elif config['nDim'] == 3:
nodesCoord[:, [1,2]] = nodesCoord[:, [2, 1]]
spanDistri = np.linspace(min(nodesCoord[:,2]), max(nodesCoord[:,2]), config['nSections'])
spanDistri[0] +=0.1
spanDistri[-1] -=0.1
+ print('spanDistri', spanDistri)
else:
raise RuntimeError("Number of dimensions not recognized", config['nDim'])
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(nodesCoord[:,0], nodesCoord[:,1], nodesCoord[:,2], s=0.3)
+ #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
+ #ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+
wingSpan = max(nodesCoord[:,2]) - min(nodesCoord[:,2])
+ upperNodes, lowerNodes = self.__SeparateWing(nodesCoord)
+
+ lowerNodes[:,0] *= -1
+ transformedCoord = np.vstack((upperNodes, lowerNodes))
+
+ tck_wing = bisplrep(transformedCoord[:,0], transformedCoord[:,2], transformedCoord[:,1], kx=3, ky=3, s=1e-3, quiet=1)
+
+ wakeCoord = np.zeros((len(bnd_wake.nodes), 3))
+ for iNode, n in enumerate(bnd_wake.nodes):
+ wakeCoord[iNode, 0] = n.pos[0]
+ wakeCoord[iNode, 1] = n.pos[1]
+ wakeCoord[iNode, 2] = n.pos[2]
+
+ if config['nDim'] == 3:
+ wakeCoord[:, [1,2]] = wakeCoord[:, [1,2]]
+
+ #tck_wake = splrep(wakeCoord)
+
+ for iSec, z in enumerate(spanDistri):
+
+ portion = nodesCoord[abs(nodesCoord[:,2] - z) <= 0.01 * wingSpan,:]
+
+ chordLength = max(portion[:, 0]) - min(portion[:, 0])
+ xDistri = min(portion[:,0]) + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
+
+ viscStruct[iSec][0][:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
+ viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ viscStruct_tr[iSec][0][:,0] = np.concatenate((xDistri[:-1], -xDistri[::-1]), axis=None).copy()
+ viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
+
+ # Interpolation in the transformed space.
+ for iPoint in range(len(viscStruct_tr[iSec][0][:,0])):
+ viscStruct_tr[iSec][0][iPoint,1] = bisplev(viscStruct_tr[iSec][0][iPoint,0], viscStruct_tr[iSec][0][iPoint,2], tck_wing)
+
+ viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(transformedCoord[:,0], transformedCoord[:,2], transformedCoord[:,1], s=0.3)
+ #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)``
+ for iSec in range(len(viscStruct)):
+ ax.scatter(viscStruct_tr[iSec][0][:,0], viscStruct_tr[iSec][0][:,2], viscStruct[iSec][0][:,1], color='red')
+ ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+
+ #viscStruct_tr[iSec][0][:,1] = splev(viscStruct_tr[iSec][0][:,0], spl)
+ #viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
+ #viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
+
+ """# Wake
+ wakeLength = np.max(wakeCoord[:,0]) - np.max(viscStruct[iSec][0][:,0])
+ viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
+ viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
+ viscStruct[iSec][1][:,1] = interp1d(wakeCoord[:,0], wakeCoord[:,1])(viscStruct[iSec][1][:,0])"""
+
+ fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(nodesCoord[:,0], nodesCoord[:,2], nodesCoord[:,1], s=0.3)
+ #ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
+ for iSec in range(len(viscStruct)):
+ ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,2], viscStruct[iSec][0][:,1], color='red')
+ ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+ return nodesList, transformedCoord, wakeCoord, viscStruct, viscStruct_tr
+
+
+ def __GetAirfoil(self, bnd_wing, bnd_wake, config):
+
+ if config['nSections'] != 1:
+ raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', config['nSections'])
+
+ viscStruct = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(2)] for _ in range(config['nSections'])]
+ viscStruct_tr = [[np.zeros((config['nPoints'][1] if i==1 else 2*config['nPoints'][0]-1, 3)) for i in range(1)] for _ in range(config['nSections'])]
+
+ # Create node list
+ nodesCoord = np.zeros((0, 3))
+ for n in bnd_wing.nodes:
+ nodesCoord = np.vstack((nodesCoord, [n.pos[0], n.pos[1], n.pos[2]]))
+
+ if config['nDim'] == 2:
+ spanDistri = np.zeros(1)
+ spanDistri[0] = nodesCoord[0,2] # Span direction is the z direction
+
+ elif config['nDim'] == 3:
+ nodesCoord[:, [1,2]] = nodesCoord[:, [2, 1]]
+ spanDistri = np.linspace(min(nodesCoord[:,2]), max(nodesCoord[:,2]), config['nSections'])
+ spanDistri[0] +=0.1
+ spanDistri[-1] -=0.1
+ print('spanDistri', spanDistri)
+
+ else:
+ raise RuntimeError("Number of dimensions not recognized", config['nDim'])
+
+ wingSpan = max(nodesCoord[:,2]) - min(nodesCoord[:,2])
+
+ upperNodes, lowerNodes = self.__SeparateWing(nodesCoord)
+
for iSec, z in enumerate(spanDistri):
portion = nodesCoord[abs(nodesCoord[:,2] - z) <= 0.01 * wingSpan,:]
@@ -362,7 +481,11 @@ class Interpolator:
viscStruct_tr[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
# Interpolation in the transformed space.
- viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
+ splineTr = transformedCoord.copy()
+ splineTr = splineTr[splineTr[:, 0].argsort()]
+ spl = splrep(splineTr[:,0], splineTr[:,1], k=5)
+ viscStruct_tr[iSec][0][:,1] = splev(viscStruct_tr[iSec][0][:,0], spl)
+ #viscStruct_tr[iSec][0][:,1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(viscStruct_tr[iSec][0][:,0])
viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
# Wake
@@ -371,151 +494,13 @@ class Interpolator:
viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
viscStruct[iSec][1][:,1] = interp1d(wakeCoord[:,0], wakeCoord[:,1])(viscStruct[iSec][1][:,0])
-
-
return nodesList, transformedCoord, wakeCoord, viscStruct, viscStruct_tr
- """for iSec, z in enumerate(spanDistri):
-
- # Create x distribution at the corresponding span position
- minX = min(upperCoord[abs(upperCoord[:,2]-z)<=config['eps'], 0])
- chordLength = max(upperWing[abs(upperWing[:,2]-z)<=config['eps'], 0]) - minX
- xDistri = minX + 1/2*(1 + np.cos(np.linspace(0, np.pi, config['nPoints'][0])))*chordLength
- nElmSide = len(xDistri)
- viscStruct[iSec][0][:,0] = np.concatenate((xDistri, xDistri[::-1]), axis=None)
- viscStruct[iSec][0][:,2] = z*np.ones(len(viscStruct[iSec][0][:,0]))
-
- wakeLength = np.max(wakeCoord[:,0]) - np.max(viscStruct[iSec][0][:,0])
- viscStruct[iSec][1][:,0] = np.max(viscStruct[iSec][0][:,0]) + (np.max(viscStruct[iSec][0][:,0]) + np.cos(np.linspace(np.pi, np.pi/2, config['nPoints'][1]))) * wakeLength
- viscStruct[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
-
- if config['nDim'] == 2:
- transformedViscCoord = np.zeros((np.shape(viscStruct[iSec][0])))
- transformedCoord
- transformedViscCoord[:, 1] = interp1d(transformedCoord[:,0], transformedCoord[:,1])(transformedCoord[:,0])
- viscStruct[iSec][0][:nElmSide,1] = interp1d(upper[:,0], upper[:,1])(viscStruct[iSec][0][:nElmSide,0])
- viscStruct[iSec][0][nElmSide:,1] = interp1d(lower[:,0], lower[:,1])(viscStruct[iSec][0][nElmSide:,0])
- viscStruct[iSec][0] = np.delete(viscStruct[iSec][0], nElmSide, axis=0)
-
- viscStruct[iSec][1][:,1] = interp1d(invStruct[1][:,0], invStruct[1][:,1])(viscStruct[iSec][1][:,0])
-
- plt.plot(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], '-')
- plt.plot(viscStruct[iSec][0][:nElmSide,0], viscStruct[iSec][0][:nElmSide,1], 'x')
- plt.plot(viscStruct[iSec][0][nElmSide:,0], viscStruct[iSec][0][nElmSide:,1], 'o')
- plt.plot(viscStruct[iSec][1][:,0], np.zeros(len(viscStruct[iSec][1][:,0])), 'x-')
- plt.xlim([-0.5,2])
- plt.ylim([-0.5, 0.5])
- plt.show()
-
- elif config['nDim'] == 3:
- nPointsMin = (config['k'][0]+1) * (config['k'][1]+1) # See Scipy bisplrep doc
-
- wingSpan = max(invStruct[0][:,2]) - min(invStruct[0][:,2])
- portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.01 * wingSpan,:]
-
- portionUpper, portionLower = self.__SeparateArf(portion[:,0], portion[:,1], portion[:,2])
- if len(portionUpper)<=nPointsMin or len(portionLower)<=nPointsMin:
- print(ccolors.ANSI_YELLOW,'dart::Interpolator Warning: Found {:.0f} points within 1% of spanwise position z={:.2f}. Trying 10%. '.format(
- len(portion), z),ccolors.ANSI_RESET)
- portion = invStruct[0][abs(invStruct[0][:,2] - z) < 0.1 * wingSpan,:]
- if len(portion)<=nPointsMin:
- print(ccolors.ANSI_RED, 'Failed. Continuing\n', ccolors.ANSI_RESET)
- else:
- print(ccolors.ANSI_GREEN, 'Sucess.\n', ccolors.ANSI_RESET)
-
- # Upper
- sWing = (len(portionUpper[:,0])-np.sqrt(2*len(portionUpper[:,0])) + len(portionUpper[:,0])+np.sqrt(2*len(portionUpper[:,0])))/2
- tckWingUp = bisplrep(portionUpper[:,0], portionUpper[:,2], portionUpper[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
- for iPoint in range(nElmSide):
- viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
-
- # Lower
- sWing = (len(portionLower[:,0])-np.sqrt(2*len(portionLower[:,0])) + len(portionLower[:,0])+np.sqrt(2*len(portionLower[:,0])))/2
- tckWingUp = bisplrep(portionLower[:,0], portionLower[:,2], portionLower[:,1], kx=config['k'][0], ky=config['k'][1], s=sWing)
- for iPoint in range(nElmSide, len(viscStruct[iSec][0][:,0])):
- viscStruct[iSec][0][iPoint,1] = bisplev(viscStruct[iSec][0][iPoint,0], viscStruct[iSec][0][iPoint,2], tckWingUp)
-
- invStruct[1] = invStruct[1][invStruct[1][:, 0].argsort()]
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(invStruct[1][:,0], invStruct[1][:,1], invStruct[1][:,2], s=0.2)
- ax.set_ylim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
- # Wake
- sWake = (len(invStruct[1][:,0])-np.sqrt(2*len(invStruct[1][:,0])) + len(invStruct[1][:,0])+np.sqrt(2*len(invStruct[1][:,0])))/2
- tckWake = bisplrep(invStruct[1][:,0], invStruct[1][:,2], invStruct[1][:,1], kx=config['k'][0], ky=config['k'][1], s=sWake)
- for iPoint in range(len(viscStruct[iSec][1][:,0])):
- #viscStruct[iSec][1][iPoint,1] = bisplev(viscStruct[iSec][1][iPoint,0], viscStruct[iSec][1][iPoint,2], tckWake)
- viscStruct[iSec][1][iPoint,1] = 0
-
- if config['nDim'] ==3:
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(invStruct[0][:,0], invStruct[0][:,1], invStruct[0][:,2], s=0.2)
- for iSec in range(len(viscStruct)):
- ax.scatter(viscStruct[iSec][0][:,0], viscStruct[iSec][0][:,1], viscStruct[iSec][0][:,2], s=0.9, color='red')
- ax.set_ylim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
- return invStruct, viscStruct"""
-
- def __SeparateArf(self, xInp, yInp, zInp=None):
- if zInp is None:
- zInp = np.zeros(len(xInp))
- if len(xInp) != len(yInp) or len(xInp) != len(zInp):
- raise RuntimeError ('dart::Interpolator::__SeparateArf Missmatch between intput vector sizes.')
-
- Q=np.zeros((len(xInp), 3))
- Q[:,0] = xInp.copy()
- Q[:,1] = yInp.copy()
- Q[:,2] = zInp.copy()
-
- Q = Q[Q[:, 0].argsort()]
-
- le = Q[0,:]
- te = Q[-1,:]
- ymean = (le[1] + te[1]) / 2
-
- upper, lower, save = le, le, le
-
- flag1 = 0
- iPoint = 0
- for i in range(len(Q[:,0])):
- if Q[iPoint,1] == 0.0:
- if flag1 == 1:
- Q=np.delete(Q, iPoint, axis = 0)
- iPoint-=1
- else:
- flag1 = 1
- iPoint +=1
-
- for iPoint in range(1, len(Q[:,0])):
- if Q[iPoint,1] > ymean:
- upper = np.vstack((upper, Q[iPoint,:]))
- elif Q[iPoint, 1] < ymean:
- lower = np.vstack((lower, Q[iPoint,:]))
- else:
- save = np.vstack((save, Q[iPoint, :]))
-
- _, idx = np.unique(save, return_index=True, axis=0)
- save=save[np.sort(idx)]
-
- #save = np.delete(save, 0, axis=0)
-
+ def __SeparateWing(self, nodesCoord):
- try:
- if len(save[:,0])>0:
- upper = np.vstack((upper, save))
- lower = np.vstack((lower, save))
- except:
- pass
-
- upper = np.vstack((upper, te))
- lower = np.vstack((lower, te))
+ print(np.max(nodesCoord[:,0]))
+ meany = nodesCoord[np.argmax(nodesCoord[:,0]), 1] - nodesCoord[np.argmin(nodesCoord[:,0]), 1]
- return upper, lower
+ upperNodes = nodesCoord[nodesCoord[:,1]>=meany].copy()
+ lowerNodes = nodesCoord[nodesCoord[:,1]<=meany].copy()
+ return upperNodes, lowerNodes
\ No newline at end of file
diff --git a/dart/pyVII/vStructures.py b/dart/pyVII/vStructures.py
new file mode 100644
index 0000000..fd6f7f5
--- /dev/null
+++ b/dart/pyVII/vStructures.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+class viscStruct:
+ def __init__(self, nPoints_surf, nPoints_wake):
+ self.nSurfNodes = nPoints_surf
+ self.nWakeNodes = nPoints_wake
+
+ self.nodesReal = np.zeros((self.nSurfNodes, 3)) # Coordinates of surface nodes in the real space
+ self.nodesTrans = np.zeros((self.nWakeNodes, 3)) # Coordinates of surface nodes in the transformed space
+
+ # Usefull variables for section representation
+ self.chordLength = 0.
+ self.leCoord = 0.
+
+
+ self.nodesWake = np.zeros((nPoints_wake, 3))
+
+ def GenerateDistri(self):
+ xDistri = self.le_coord + 1/2*(1 + np.cos(np.linspace(0, np.pi, len(self.nSurfNodes))))*self.chordLength
+ self.nodesReal[:,0] = np.concatenate((xDistri[:-1], xDistri[::-1]), axis=None)
+
+class invStruct:
+ def __init__(self, _bndSurf, _bndWake, _nDim):
+ self.nDim = _nDim
+
+ self.bndSurf = _bndSurf
+ self.bndWake = _bndWake
+
+ self.nSurfNodes = self.bndSurf.nodes.size()
+
+ # np matrix of nodes coordinate in the transformed space associated to the c++ structures 'node' in nodesList
+ # (coordinates at indice i in 'nodesTrans' correspond to the node at position i in 'nodesList')
+ self.nodesTrans = np.zeros(self.nSurfNodes,3)
+ self.nodesList = [None for _ in range(self.nSurfNodes)]
+
+ def CreateList(self):
+
+ nodesReal = np.zeros((self.nSurfNodes, 3))
+ nodeList = [None for _ in range(self.nSurfNodes)]
+ for iNode, n in enumerate(self.bndSurf.nodes):
+ nodesReal[iNode, :] = [n.pos[0], n.pos[1], n.pos[2]]
+ nodeList[iNode] = n
+
+
+
+
diff --git a/dart/pyVII/vUtils.py b/dart/pyVII/vUtils.py
index 32410ed..a8423fd 100644
--- a/dart/pyVII/vUtils.py
+++ b/dart/pyVII/vUtils.py
@@ -1,7 +1,23 @@
import dart
-
-def initBL(Re, Minf, CFL0, nSections):
- solver = dart.ViscSolver(Re, Minf, CFL0, nSections)
+from fwk.coloring import ccolors
+
+def initBL(Re, Minf, CFL0, nSections, verbose=1):
+ if Re<=0.:
+ print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Reynolds number.", Re, ccolors.ANSI_RESET)
+ raise RuntimeError("Invalid parameter")
+ if Minf<0.:
+ print(ccolors.ANSI_RED, "dart::vUtils Error : Negative Mach number.", Minf, ccolors.ANSI_RESET)
+ raise RuntimeError("Invalid parameter")
+ elif Minf>=1.:
+ print(ccolors.ANSI_YELLOW, "dart::vUtils Warning : (Super)sonic freestream Mach number.", Minf, ccolors.ANSI_RESET)
+ if nSections < 0:
+ print(ccolors.ANSI_RED, "dart::vUtils Fatal error : Negative number of sections.", nSections, ccolors.ANSI_RESET)
+ raise RuntimeError("Invalid parameter")
+ if 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)
return solver
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 1c3534f..2b4d6f5 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -155,6 +155,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
if (std::isnan(normSysRes) || normSysRes/normSysRes0 > 1e-3)
{
ResetSolution(iPoint, bl, Sln0);
+ return -1;
}
return innerIters;
} // End Integration
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 9ffd72a..15daf4a 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -230,30 +230,19 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
timeMatrix(4, 4) = 1.0;
}
- /*Vector<double, 5> result = timeMatrix.partialPivLu().solve(spaceVector);
- if (std::isnan(result.norm()))
- {
- std::cout << "Point " << iPoint << " det " << timeMatrix.determinant() << std::endl;
- std::cout << timeMatrix << std::endl;
- std::cout << spaceVector << std::endl;
- std::cout << dT_dx << " " << (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx << " " << " cf " << - Cfa / 2 << std::endl;
- std::cout << " u[0] " << u[0] << " u[1] " << u[1] << " u[2] " << u[2] << " u[3] " << u[3] << " u[4] " << u[4] << " Usa " << Usa << std::endl;
- throw;
- }*/
-
return timeMatrix.partialPivLu().solve(spaceVector);
}
+
double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) const
{
- double Ax = 0.;
-
double Dgr = 0.08;
double Hk1 = 3.5;
double Hk2 = 4.;
- double Hmi = (1. / (Hk - 1.));
- double logRet_crit = 2.492 * pow(1 / (Hk - 1), 0.43) + 0.7 * (tanh(14 * Hmi - 9.24) + 1.0); // Value at which the amplification starts to grow
- if (Ret <= 0.)
+ double Hmi = (1/(Hk-1));
+ double logRet_crit = 2.492*std::pow(1/(Hk-1.),0.43) + 0.7*(std::tanh(14*Hmi-9.24)+1.0); // value at which the amplification starts to grow
+ double Ax = 0.;
+ if (Ret <=0.)
{
Ret = 1.;
}
@@ -263,49 +252,52 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
}
else
{
- double Rnorm = (std::log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
+ double Rnorm = (std::log10(Ret) - (logRet_crit-Dgr))/(2*Dgr);
double Rfac = 0.;
- if (Rnorm <= 0.)
+ if (Rnorm <=0)
{
- Rfac = 0;
+ Rfac = 0.;
}
- else if (Rnorm >= 1)
+ if (Rnorm >= 1)
{
- Rfac = 1;
+ Rfac = 1.;
}
else
{
- Rfac = 3 * Rnorm * Rnorm - 2 * Rnorm * Rnorm * Rnorm;
+ Rfac = 3*Rnorm*Rnorm - 2*Rnorm*Rnorm*Rnorm;
}
- // Normal envelope amplification rate
- double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi + 0.1 * std::exp(-20. * Hmi);
- double arg = 3.87 * Hmi - 2.52;
- double dndRet = 0.028 * (Hk - 1) - 0.0345 * std::exp(-(arg * arg));
- Ax = (m * dndRet / theta) * Rfac;
- // Set correction for separated profiles
+ // normal envelope amplification rate
+ double m = -0.05+2.7*Hmi-5.5*Hmi*Hmi+3*Hmi*Hmi*Hmi+0.1*std::exp(-20*Hmi);
+ double arg = 3.87*Hmi-2.52;
+ double dndRet = 0.028*(Hk-1)-0.0345*std::exp(-(arg*arg));
+ Ax = (m*dndRet/theta)*Rfac;
+ // set correction for separated profiles
if (Hk > Hk1)
{
- double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
+ double Hnorm = (Hk - Hk1)/(Hk2-Hk1);
double Hfac = 0.;
- if (Hnorm >= 1.)
+ if (Hnorm >=1)
{
- Hfac = 1;
+ Hfac = 1.;
}
else
{
- Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm;
+ Hfac = 3*Hnorm*Hnorm - 2*Hnorm*Hnorm*Hnorm;
}
double Ax1 = Ax;
- double Gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5));
- double Tnr = std::tanh(1.2 * (std::log10(Ret) - Gro));
- double Ax2 = (0.086 * Tnr - 0.25 / std::pow(Hk - 1., 1.5)) / theta;
- if (Ax2 < 0)
+ double Gro = 0.3+0.35*std::exp(-0.15*(Hk-5));
+ double Tnr = std::tanh(1.2*(std::log10(Ret)-Gro));
+ double Ax2 = (0.086*Tnr - 0.25/std::pow((Hk-1),1.5))/theta;
+ if (Ax2 < 0.)
{
- Ax2 = 0;
+ Ax2 = 0.;
+ }
+ Ax = Hfac*Ax2 + (1-Hfac)*Ax1;
+ if (Ax < 0.)
+ {
+ Ax = 0.;
}
- Ax = Hfac * Ax2 + (1 - Hfac) * Ax1;
- Ax = std::max(Ax, 0.);
}
}
return Ax;
-}
\ No newline at end of file
+}
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 65cbe38..a3d0f9b 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -28,7 +28,7 @@ using namespace dart;
/**
* @brief Viscous solver and boundary layer(s) initialization
*/
-ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections)
+ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections, unsigned int _verbose)
{
PrintBanner();
@@ -40,7 +40,7 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
Cdt = 0.; Cdf = 0.; Cdp = 0.;
- verbose = 1;
+ verbose = _verbose;
/* Initialzing boundary layer */
@@ -100,7 +100,7 @@ ViscSolver::~ViscSolver()
* @brief Runs viscous operations. Temporal solver parametrisation is first adapted,
* Solutions are initialized than time marched towards steady state successively for each points.
*/
-int ViscSolver::Run(unsigned int couplIter)
+int ViscSolver::Run(unsigned int const couplIter)
{
bool lockTrans=false;
int pointExitCode = 0;
@@ -223,7 +223,7 @@ int ViscSolver::Run(unsigned int couplIter)
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
- //Sections[iSec][iRegion]->U[iPoint * Sections[iSec][iRegion]->GetnVar() + 2] = 10.0;
+ // Sections[iSec][iRegion]->U[iPoint * Sections[iSec][iRegion]->GetnVar() + 2] = 10.0;
/* Unsucessfull convergence output */
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index d56aabe..fdf8726 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -28,9 +28,9 @@ private:
TimeSolver *tSolver;
public:
- ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections);
+ ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections, unsigned int verbose=1);
~ViscSolver();
- int Run(unsigned int couplIter);
+ int Run(unsigned int const couplIter);
double GetRe() const {return Re;};
void SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 2f0a77e..6ef0ad4 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -59,13 +59,13 @@ def main():
S_ref = 1.*spn # reference area
c_ref = 1 # reference chord
fms = 1.0 # farfield mesh size
- nms = 0.01 # nearfield mesh size
+ nms = 0.02 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : 0.001, 'msTe' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/oneraM6.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 9272fca..544b06a 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -39,6 +39,7 @@ 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
@@ -81,7 +82,7 @@ def main():
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ 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()
@@ -96,7 +97,7 @@ def main():
tms['solver'].start()
isolver = floD.newton(pbl)
- vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
+ 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)
@@ -106,8 +107,18 @@ def main():
#isolver.reset()
nomDeVariableDePorc_x[i] = coupler.sln[15].copy()
- nomDeVariableDePorc_var[i] = coupler.sln[2].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-')
diff --git a/dart/tests/bli_lowLift2.py b/dart/tests/bli_lowLift2.py
index 7e738fa..7704020 100644
--- a/dart/tests/bli_lowLift2.py
+++ b/dart/tests/bli_lowLift2.py
@@ -34,12 +34,14 @@
# 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
@@ -59,7 +61,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 0*math.pi/180
+ alpha = 2*math.pi/180
M_inf = 0.
CFL0 = 1
@@ -71,32 +73,57 @@ def main():
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.01}
- 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, nSections)
- config = {'nDim': dim, 'nSections':nSections, 'nPoints':[103, 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 20}
-
- iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
- coupler = viscC.Coupler(iSolverAPI, vsolver)
- aeroCoeff = coupler.Run()
- #coupler.RunPolar(alphaSeq)
+ nPv = [100, 100, 100]
+ aeroCoeff = [np.zeros(2) for i in range(len(nPv))]
+
+ fig, axs = plt.subplots(2,1)
+
+ 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, 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, 0)
+ config = {'nDim': dim, 'nSections':nSections, 'nPoints':[nPv[i], 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 3}
+
+ iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
+
+ coupler = viscC.Coupler(iSolverAPI, vsolver)
+ aeroCoeff[i] = coupler.Run()
+ axs[0].plot(aeroCoeff[i][:,0], 'x-')
+ axs[1].plot(aeroCoeff[i][:,1], 'x-')
+
+ del isolver
+ del vsolver
+ del iSolverAPI
+ del pbl
+ del bnd
+ del blw
+ del msh
+ del gmshWriter
+
+ axs[0].set_xlabel('x/c')
+ axs[0].set_ylabel('Cl')
+ axs[0].set_ylim([-.1, .1])
+ axs[1].set_xlabel('x/c')
+ axs[1].set_ylabel('Cd')
+ plt.show()
tms['solver'].stop()
# extract Cp
@@ -115,15 +142,11 @@ def main():
print('SOLVERS statistics')
print(coupler.tms)
- xtr=vsolver.Getxtr()
+ xtr=[vsolver.Getxtr(0,0),vsolver.Getxtr(0,1)]
- """# visualize solution and plot results
+ # 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)"""
+ #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
--
GitLab