From e6441fcd8cbee70cef89e9a14930d8e4f2cb84c4 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paulzer@sentinelle.local>
Date: Mon, 16 May 2022 15:09:04 +0200
Subject: [PATCH] Improving blw interpolation
---
dart/pyVII/vInterpolator.py | 104 ++++++++----------------------------
dart/src/wBLRegion.cpp | 1 +
dart/src/wBLRegion.h | 2 +-
dart/src/wClosures.cpp | 4 +-
dart/src/wTimeSolver.cpp | 9 ----
dart/src/wViscFluxes.cpp | 4 +-
dart/src/wViscSolver.cpp | 2 -
dart/tests/bli2.py | 4 +-
dart/tests/bli3.py | 2 +-
9 files changed, 32 insertions(+), 100 deletions(-)
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index 7323c2f..66dee72 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -1,10 +1,9 @@
-from gc import isenabled
-from xmlrpc.client import INVALID_ENCODING_CHAR
import numpy as np
from matplotlib import pyplot as plt
from fwk.coloring import ccolors
import dart.pyVII.vUtils as blxU
+import sys
#import dart.pyVII.vStructures as vStruct
# 3D spline with Bézier curves
@@ -31,8 +30,9 @@ class Interpolator:
self.tms['GetWing'].start()
self.__GetWing(_blw, _cfg)
self.tms['GetWing'].stop()
- print(self.tms)
+
self.cfg = _cfg
+ print(self.tms)
self.stgPt=np.zeros(self.cfg['nSections'], dtype=int)
@@ -74,6 +74,7 @@ class Interpolator:
Uz = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], U[0][:,2], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], U[1][:,2], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
Me = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], M[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], M[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
Rhoe = [self.__RbfToSections(self.invStruct_tr[0][:,:nD], Rho[0], self.viscStruct_tr[iSec][0][:,:nD]), self.__RbfToSections(self.invStruct_tr[1][:, [0,2] if nD==3 else 0], Rho[1], self.viscStruct[iSec][1][:,[0,2] if nD==3 else 0])]
+
if couplIter == 0:
self.stgPt[iSec] = np.argmin(Ux[0]**2+Uy[0]**2)
@@ -99,6 +100,10 @@ class Interpolator:
UyUp, UyLw = self.__Remesh(Uy[0], self.stgPt[iSec])
UzUp, UzLw = self.__Remesh(Uz[0], self.stgPt[iSec])
+ plt.plot(xUp, UxUp)
+ plt.plot(xLw, 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])
@@ -137,32 +142,15 @@ class Interpolator:
x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
blwUp = vSolver.ExtractBlowingVel(iSec, 0)
blwLw = vSolver.ExtractBlowingVel(iSec, 1)
- plt.plot(blwLw)
- plt.show()
blw[iReg] = np.concatenate((blw[iReg], blwUp[::-1] , blwLw))
- plt.plot(self.viscCg_tr[iSec][iReg][:,0],np.concatenate((blwUp[::-1] , blwLw)))
- plt.title(iSec)
- plt.show()
elif iReg == 1:
x[iReg] = np.vstack((x[iReg], self.viscCg_tr[iSec][iReg][:,:3]))
blw[iReg] = np.concatenate((blw[iReg], vSolver.ExtractBlowingVel(iSec, 2)))
self.tms['Interpolation'].start()
nD = self.cfg['nDim']
- mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3], forceLin=1), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0], forceLin=1)]
+ mappedBlw = [self.__RbfToSections(x[0], blw[0], self.invCg_tr[0][:,:3]), self.__RbfToSections(x[1][:, [0,2] if nD==3 else 0], blw[1], self.invCg_tr[1][:, [0,2] if nD==3 else 0])]
-
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(self.invCg_tr[0][:,0], self.invCg_tr[0][:,2], mappedBlw[0], s=0.3)
- #ax.scatter(abs(invCg_tr[0][:,0]), invCg_tr[0][:,2], invCg_tr[0][:,1], s=0.3)
- for sec in range(len(self.viscStruct_tr)):
- ax.scatter(x[0][:,0], x[0][:,2], blw[0], color = 'red')
- #ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
self.tms['Interpolation'].stop()
for iElm in range(len(self.invCg_tr[0][:,0])):
self.blw[0].setU(iElm, mappedBlw[0][iElm])
@@ -199,11 +187,9 @@ class Interpolator:
xLw = vec[stgPt:]
return xUp, xLw
- def __RbfToSections(self, x, var, xs, forceLin=0):
+ 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)
- elif forceLin==1:
- varOut = RBFInterpolator(x, var, neighbors=50, kernel='linear', smoothing=0.1, degree=0)(xs)
else:
varOut = RBFInterpolator(x, var, neighbors=self.cfg['neighbors'], kernel=self.cfg['rbftype'], smoothing=self.cfg['smoothing'], degree=self.cfg['degree'])(xs)
return varOut
@@ -254,7 +240,7 @@ class Interpolator:
invStruct[0] = np.vstack((upperNodes, lowerNodes))
# Remove duplicates
- invStruct_tr[0] = np.unique(invStruct_tr[0], axis=0)
+ #invStruct_tr[0] = np.unique(invStruct_tr[0], axis=0)
# Wake part
for iNode, n in enumerate(blw[1].nodes):
@@ -278,24 +264,12 @@ class Interpolator:
spanDistri[0] +=0.1
spanDistri[-1] -=0.1
- print('spanDistri', spanDistri)
-
# Spline wing
# Point in the wing tip region are not accounted for during Bezier curve spline
- #xUp = upperNodes[upperNodes[:,0]<config['span'], :]
- #xLw = lowerNodes[lowerNodes[:,0]<config['span'], :]
-
- #xUp[:,[1,2]] = xUp[:,[2,1]]
- #xLw[:,[1,2]] = xLw[:,[2,1]]
-
- #sUp = (len(xUp)-np.sqrt(2*len(xUp)) + len(xUp)+np.sqrt(2*len(xUp)))/2
- #sLw = (len(xLw)-np.sqrt(2*len(xLw)) + len(xLw)+np.sqrt(2*len(xLw)))/2
- #tck_wingUp = bisplrep(xUp[:,0], xUp[:,2], xUp[:,1], kx=3, ky=3, s=sUp, quiet=1)
- #tck_wingLw = bisplrep(xLw[:,0], xLw[:,2], xLw[:,1], kx=3, ky=3, s=sLw, quiet=1)
m = len(invStruct_tr[0][:,0])
_s = (m-np.sqrt(2*m) + m+np.sqrt(2*m))/2
- tck_wing = bisplrep(invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],0], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],2], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],1], kx=5, ky=5, s=1e-8)
+ tck_wing = bisplrep(invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],0], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],2], invStruct_tr[0][invStruct_tr[0][:,2]<config['span'],1], kx=5, ky=5, s=1e-6, quiet=1)
# Create sections
self.tms['iSecLoop'].start()
@@ -322,22 +296,22 @@ class Interpolator:
elif config['nDim'] == 2:
viscStruct_tr[iSec][0][:,1] = fInterp(viscStruct_tr[iSec][0][:,0])
- viscStruct_tr[iSec][0][0,1] = 0.
- viscStruct_tr[iSec][0][-1,1] = 0.
+ #viscStruct_tr[iSec][0][0,1] = 0.
+ #viscStruct_tr[iSec][0][-1,1] = 0.
viscStruct[iSec][0][:,1] = viscStruct_tr[iSec][0][:,1].copy()
# Wake
wakeLength = np.max(invStruct[1][:,0]) - np.max(viscStruct[iSec][0][:,0])
viscStruct_tr[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_tr[iSec][1][:,2] = z*np.ones(len(viscStruct[iSec][1][:,0]))
- viscStruct_tr[iSec][1][:,1] = interp1d(invStruct_tr[1][:,0], invStruct_tr[1][:,1], bounds_error=False, fill_value='extrapolate')(viscStruct_tr[iSec][1][:,0])
+ viscStruct_tr[iSec][1][:,1] = interp1d(invStruct_tr[1][:,0], invStruct_tr[1][:,1])(viscStruct_tr[iSec][1][:,0])
viscStruct[iSec][1] = viscStruct_tr[iSec][1].copy()
# Cg positions
self.tms['iSecLoop'].stop()
- fig = plt.figure()
+ """fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct_tr[0][:,1], s=0.3)
ax.scatter(invStruct_tr[1][:,0], invStruct_tr[1][:,2], invStruct_tr[1][:,1], s=0.3)
@@ -348,17 +322,6 @@ class Interpolator:
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
- plt.show()
- """fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(invStruct[0][:,0], invStruct[0][:,2], invStruct[0][:,1], s=0.3)
- for sec in viscStruct:
- ax.scatter(sec[0][:,0], sec[0][:,2], sec[0][:,1], color = 'red')
- #ax.scatter(sec[1][:,0], sec[1][:,2], sec[1][:,1], color = 'red')
- ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
plt.show()"""
viscCg_tr = [[np.zeros((len(viscStruct[iSec][iReg][:,0]) - 1, 3)) for iReg in range(2)] for iSec in range(len(viscStruct))]
@@ -390,38 +353,15 @@ class Interpolator:
invCg_tr[iReg][:,[1,2]] = invCg_tr[iReg][:,[2,1]]
self.tms['CgLoop'].stop()
- """if config['nDim'] == 3:
- fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(invStruct_tr[0][:,0], invStruct_tr[0][:,2], invStruct[0][:,1], s=0.3)
- #ax.scatter(abs(invCg_tr[0][:,0]), invCg_tr[0][:,2], invCg_tr[0][:,1], s=0.3)
- for sec in viscStruct_tr:
- ax.scatter(abs(sec[0][:,0]), sec[0][:,2], sec[0][:,1], color = 'red')
- ax.scatter(abs(sec[1][:,0]), sec[1][:,2], sec[1][:,1], color = 'red')
- #ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
- elif config['nDim'] == 2:
- plt.plot(abs(invStruct_tr[0][:,0]), invStruct[0][:,1], 'x', color='blue')
- plt.plot(abs(invStruct_tr[1][:,0]), invStruct[1][:,1], 'x', color='blue')
- plt.plot(abs(invCg_tr[0][:,0]), invCg_tr[0][:,1], 'o', color='red')
- plt.plot(abs(invCg_tr[1][:,0]), invCg_tr[1][:,1], 'o', color='red')
- plt.title('Inviscid')
- plt.show()
- plt.plot(viscStruct_tr[0][0][:,0], viscStruct_tr[0][0][:,1])
- plt.plot(viscStruct_tr[0][1][:,0], viscStruct_tr[0][1][:,1])
- plt.plot(invStruct_tr[0][:,0], invStruct[0][:,1], 'o', color='blue')
- plt.plot(invStruct_tr[1][:,0], invStruct[1][:,1], 'x', color='blue')
- plt.title('Viscous vs Inviscid')
- plt.show()"""
-
-
self.invStruct = invStruct
self.invStruct_tr = invStruct_tr
self.viscStruct = viscStruct
self.viscStruct_tr = viscStruct_tr
self.invCg_tr = invCg_tr
- self.viscCg_tr = viscCg_tr
\ No newline at end of file
+ self.viscCg_tr = viscCg_tr
+ np.set_printoptions(threshold=sys.maxsize)
+ for sec in self.viscStruct:
+ print('next')
+ sec[0][:,[1,2]] = sec[0][:,[2,1]]
+ print(sec[0])
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 0998854..f87aadb 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -161,6 +161,7 @@ void BLRegion::ComputeBlwVel()
void BLRegion::PrintSolution(size_t iPoint) const
{
std::cout << "Pt " << iPoint << "Reg " << Regime[iPoint] << std::endl;
+ std::cout << "x " << mesh->Getx(iPoint) << "xx " << mesh->GetLoc(iPoint) << "xxExt " << mesh->GetExt(iPoint) << std::endl;
std::cout << std::scientific << std::setprecision(15);
std::cout << std::setw(10) << "T " << U[iPoint*nVar + 0] << "\n"
<< std::setw(10) << "H " << U[iPoint*nVar + 1] << "\n"
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index 00e111c..afd73d5 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -76,4 +76,4 @@ public:
};
} // namespace dart
-#endif //WBLREGION_H
\ No newline at end of file
+#endif //WBLREGION_H
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 7f867be..3e29976 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -108,10 +108,10 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const
{
- if (std::isnan(Ret))
+ /*if (std::isnan(Ret))
{
std::cout << "Ret is nan. H:" << H << std::endl;
- }
+ }*/
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
//Ct = std::max(std::min(Ct, 0.30), 0.0000001);
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 2b4d6f5..3cae23f 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -99,17 +99,8 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
unsigned int nErrors = 0; // Number of errors encountered
unsigned int innerIters = 0; // Inner (non-linear) iterations
- /* for (int i=0; i<1e10; ++i)
- {
- double a = 1+1;
- } */
while (innerIters < maxIt)
{
- /* if (bl->name=="upper" && iPoint == 1 && innerIters == 0)
- {
- std::cout << "it " << innerIters << " res " << std::log10(normSysRes/normSysRes0) << std::endl;
- bl->PrintSolution(0);
- } */
if (innerIters % itFrozenJac == 0)
{
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 0346a58..c940708 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -75,10 +75,10 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
dissipRatio = 1.;
} // End else
- if (std::isnan(u[1]))
+ /*if (std::isnan(u[1]))
{
std::cout << "H is nan. Point " << iPoint << bl->name << std::endl;
- }
+ }*/
bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); /* Reynolds number based on theta */
bl->DeltaStar[iPoint] = u[0] * u[1]; /* Displacement thickness */
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 2a06c2c..1dc6e70 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -223,8 +223,6 @@ int ViscSolver::Run(unsigned int const couplIter)
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
- // Sections[iSec][iRegion]->U[iPoint * Sections[iSec][iRegion]->GetnVar() + 2] = 10.0;
-
/* Unsucessfull convergence output */
if (pointExitCode != 0)
diff --git a/dart/tests/bli2.py b/dart/tests/bli2.py
index 26c8a75..e16ea38 100644
--- a/dart/tests/bli2.py
+++ b/dart/tests/bli2.py
@@ -60,7 +60,7 @@ def main():
# define flow variables
Re = 1e7
alpha = 2*math.pi/180
- M_inf = 0.715
+ M_inf = 0.73
CFL0 = 1
@@ -107,6 +107,8 @@ def main():
# 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
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 792d16b..718ba48 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -77,7 +77,7 @@ def main():
tms['pre'].stop()
# solve problem
- config = {'nDim': dim, 'nSections':nSections, 'span':spn, 'nPoints':[200, 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.001, 'degree': 1, 'neighbors': 5}
+ config = {'nDim': dim, 'nSections':nSections, 'span':spn, 'nPoints':[100, 25], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 1e-8, 'degree': 1, 'neighbors': 50}
iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw, config)
vSolver = viscU.initBL(Re, M_inf, CFL0, nSections, 2)
#iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
--
GitLab