diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 887e692cf6f87f5435a20da24803f444b6da4683..05c672533a6f1455e5cdea974fcfe5fdf38e11db 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -55,6 +55,8 @@ class Coupler:
# Initilize meshes in c++ objects
+ dragIter = np.zeros((0,2))
+
self.nElmUpperPrev = 0
couplIter = 0
cdPrev = 0
@@ -87,16 +89,17 @@ class Coupler:
self.tms['processing'].stop()
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]))
if error <= self.couplTol and exitCode==0:
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 0
+ return dragIter
cdPrev = self.vSolver.Cdt
-
if couplIter%1==0:
xtr = [1, 1]
xtr[0]=self.vSolver.Getxtr(0,0)
@@ -108,7 +111,7 @@ class Coupler:
self.tms['processing'].stop()
else:
print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
- return couplIter
+ return dragIter
def RunPolar(self, alphaSeq):
@@ -179,7 +182,7 @@ class Coupler:
dataIn=[None,None]
for n in range(0,len(self.group)):
- self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList()
+ self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList(couplIter)
fMeshDict, _, _ = GroupConstructor().mergeVectors(dataIn, 0, 1)
if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
index d0f24b4efa984dc0dabbcb4789f8f8d0bfc44925..355b106d857fa3a5025ecceb87959d429dd3e2e7 100644
--- a/dart/pyVII/vCoupler2.py
+++ b/dart/pyVII/vCoupler2.py
@@ -34,7 +34,7 @@ class Coupler:
self.iSolverAPI = _iSolverAPI
self.vSolver = vSolver
- self.maxCouplIter = 150
+ self.maxCouplIter = 25
self.couplTol = 1e-4
self.tms=fwk.Timers()
@@ -46,7 +46,7 @@ class Coupler:
# Initilize meshes in c++ objects
- dragIter = []
+ dragIter = np.zeros((0,2))
self.nElmUpperPrev = 0
couplIter = 0
@@ -79,7 +79,7 @@ class Coupler:
self.iSolverAPI.SetBlowingVel(self.vSolver)
self.tms['processing'].stop()
- dragIter.append(self.vSolver.Cdt)
+ dragIter = np.vstack((dragIter, [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:
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index 5542ce971bfd5e4316ad26cf2aa1477ffee37250..abe70e5563f8fdf6c82ef81cbd34b72ce08bb953 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -32,6 +32,7 @@ class Interpolator:
self.bndWing = bnd_wing
self.bndWake = bnd_wake
self.invStruct, self.viscStruct = self.__GetWing(bnd_wing, bnd_wake, _cfg)
+ print('Sections created')
self.config = _cfg
self.xxPrev = [[np.zeros(0) for iReg in range(3)] for iSec in range(len(self.viscStruct))]
@@ -123,8 +124,8 @@ class Interpolator:
plt.plot(xLw, UyLw, 'o')
#plt.plot(self.viscStruct[iSec][1][:,0], Uy[iSec][1])
plt.plot(self.invStruct[0][:,0], U_wg[:,1], 'x')
- #plt.xlim([0.4, 1.2])
- #plt.ylim([-0.3, 0.2])
+ plt.xlim([0.4, 1.2])
+ plt.ylim([-0.3, 0.2])
plt.title('Uy = f(x) smooting{:.2e}'.format(self.config['smoothing']))
plt.show()"""
@@ -207,19 +208,24 @@ class Interpolator:
"""plt.plot(self.invStruct[0][:,0], self.invStruct[0][:,1], 'x')
plt.plot(invCgWing[:, 0], invCgWing[:,1], 'o')
- plt.show()"""
+ plt.show()
- invCgWingUp = invCgWing[invCgWing[:,1]>0, :]
- invCgWingLw = invCgWing[invCgWing[:,1]<0, :]
+ plt.plot(self.viscStruct[0][0][:,0], self.viscStruct[0][0][:,1], 'x')
+ plt.plot(viscCgWing[:,0], viscCgWing[:,1], 'o')
+ plt.show()"""
- mappedBlwWing = np.zeros(len(invCgWingUp) + len(invCgWingLw))
self.tms['Interpolation'].start()
+ invCgUp = invCgWing[invCgWing[:,1]>0, :]
+ invCgLw = invCgWing[invCgWing[:,1]<0, :]
+ #mappedBlwWing = np.zeros(len(invCgWing))
+ mappedBlwWing = RBFInterpolator(viscCgWing[:,:2], BlwWing, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(invCgWing[:,:2])
#mappedBlwWing[:self.config['nPoints'][0]] = RBFInterpolator(viscCgWing[:self.config['nPoints'][0],:], BlwWing[:self.config['nPoints'][0]], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=1e-6, degree=self.config['degree'])(invCgWingUp)
#mappedBlwWing[self.config['nPoints'][0]:] = RBFInterpolator(viscCgWing[self.config['nPoints'][0]:, :], BlwWing[self.config['nPoints'][0]:], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=1e-6, degree=self.config['degree'])(invCgWingLw)
- mappedBlwWing[:len(invCgWingUp)] = RBFInterpolator(viscCgWing[:self.config['nPoints'][0],:], BlwWing[:self.config['nPoints'][0]], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWingUp)
- mappedBlwWing[len(invCgWingUp):] = RBFInterpolator(viscCgWing[self.config['nPoints'][0]:, :], BlwWing[self.config['nPoints'][0]:], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWingLw)
+ # Following lines (interp1d) will not work.
+ #mappedBlwWing[:self.config['nPoints'][0]] = interp1d(viscCgWing[:self.config['nPoints'][0],0], BlwWing[:self.config['nPoints'][0]])(invCgUp[:, 0])
+ #mappedBlwWing[self.config['nPoints'][0]:] = interp1d(viscCgWing[self.config['nPoints'][0]:, 0], BlwWing[self.config['nPoints'][0]:])(invCgLw[:, 0])
#mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(invCgWake)
- mappedBlwWake = RBFInterpolator(viscCgWake, BlwWake, neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWake)
+ mappedBlwWake = RBFInterpolator(viscCgWake[:, :2], BlwWake, neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(invCgWake[:, :2])
"""plt.plot(invCgWing[:,0], invCgWing[:,1])
plt.show()"""
@@ -286,43 +292,21 @@ class Interpolator:
delIdx = np.where(completeWing[:,0]==1.0)
delIdx = delIdx[0]
save = completeWing[delIdx, 3]
- print('save', save)
- completeWing = np.delete(completeWing, delIdx[0], axis = 0)
- completeWingUp = completeWing[completeWing[:,1]>0,:]
- completeWingLw = completeWing[completeWing[:,1]<0,:]
-
- completeWingUp = np.vstack((completeWingUp, completeWing[delIdx[0], :]))
- completeWingLw = np.vstack((completeWingLw, completeWing[delIdx[0], :]))
+ completeWing[delIdx[0], 3] = (completeWing[delIdx[0], 3] + completeWing[delIdx[1], 3]) / 2
+ completeWing = np.delete(completeWing, delIdx[1], axis = 0)
mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.viscStruct))]
for iSec in range(len(mappedVar)):
for iReg in range(len(mappedVar[iSec])):
if iReg==0: # Wing
- mappedVar[iSec][iReg] = np.zeros(len(self.viscStruct[iSec][0][:,0]))
- mappedVar[iSec][iReg][:self.config['nPoints'][0]] = RBFInterpolator(completeWingUp[:,:2], completeWingUp[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][0][:self.config['nPoints'][0],:2])
- mappedVar[iSec][iReg][self.config['nPoints'][0]:] = RBFInterpolator(completeWingLw[:,:2], completeWingLw[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][0][self.config['nPoints'][0]:,:2])
-
- #interpWWW = NearestNDInterpolator(list(zip(completeWing[:,0], completeWing[:,1], completeWing[:,2])), completeWing[:,3])
- #mappedVar[iSec][iReg] = interpWWW(self.viscStruct[iSec][0][:,0], self.viscStruct[iSec][0][:,1], self.viscStruct[iSec][0][:,2])
-
- #interp = CloughTocher2DInterpolator(list(zip(completeWing[:,0], completeWing[:,1])), completeWing[:,3])
- #mappedVar[iSec][iReg] = interp(self.viscStruct[0][0][:,0], self.viscStruct[0][0][:,1])
- #mappedVar[iSec][iReg] = griddata(completeWing[:,:2], completeWing[:,3], self.viscStruct[iSec][0][:,:2])
+ mappedVar[iSec][iReg] = RBFInterpolator(completeWing[:,:2], completeWing[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][0][:,:2])
if iReg==1: # Wake
- #mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=self.config['neighbors'], kernel=self.config['rbftype'], smoothing=self.config['smoothing'], degree=self.config['degree'])(self.viscStruct[iSec][1])
mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:3], completeWake[:,3], neighbors=1, kernel='linear', smoothing=1e-6, degree=0)(self.viscStruct[iSec][1])
-
- #interpWWW = NearestNDInterpolator(list(zip(completeWake[:,0], completeWake[:,1], completeWake[:,2])), completeWake[:,3])
- #mappedVar[iSec][iReg] = interpWWW(self.viscStruct[iSec][1][:,0], self.viscStruct[iSec][1][:,1], self.viscStruct[iSec][1][:,2])
-
- #mappedVar[iSec][iReg] = griddata(completeWake[:,:2], completeWake[:,3], self.viscStruct[iSec][1][:,:2], method='linear')
- pass
-
#mappedVar[iSec][0][0] = save[0]
#mappedVar[iSec][0][-1] = save[-1]
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 6762667511945d2018035a55e4c49afc18d51721..1560e34e5751e31007d948bf287ed5ffa5e03e71 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -39,6 +39,8 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
CFL0 = _CFL0;
+ solverExitCode = 0;
+
/* Initialzing boundary layer */
Sections.resize(nSections);
for (size_t iSec=0; iSec < nSections; ++iSec)
@@ -97,7 +99,6 @@ ViscSolver::~ViscSolver()
int ViscSolver::Run(unsigned int couplIter)
{
bool lockTrans;
- int solverExitCode = 0;
int pointExitCode;
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
@@ -126,6 +127,14 @@ int ViscSolver::Run(unsigned int couplIter)
Sections[iSec][1]->SetStagBC(Re);
+ bool reset = false;
+ if (solverExitCode != 0)
+ {
+ reset = true;
+ }
+ solverExitCode = 0;
+
+
/* Loop over the regions */
for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
{
@@ -135,7 +144,7 @@ int ViscSolver::Run(unsigned int couplIter)
/* Check if we need to enter safe mode */
- if (couplIter == 0 || Sections[iSec][iRegion]->mesh->GetDiffnElm() != 0 || stagPtMvmt[iSec] == true)
+ if (couplIter == 0 || Sections[iSec][iRegion]->mesh->GetDiffnElm() != 0 || stagPtMvmt[iSec] == true || reset)
{
tSolver->SetCFL0(1);
tSolver->SetitFrozenJac(1);
@@ -203,7 +212,7 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
/* Time integration */
- /*if (iRegion == 1 && iPoint == 1)
+ if (iRegion == 1 && iPoint == 1)
{
std::scientific;
std::cout << " iSec " << iSec
@@ -215,7 +224,7 @@ int ViscSolver::Run(unsigned int couplIter)
<< " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
<< std::endl;
Sections[iSec][iRegion]->PrintSolution(iPoint);
- }*/
+ }
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index a998acf0ebe9a3bf51fcf31a2daf26139a29ed94..bf0bd953683c1665411e5ee99b59f962166b3e83 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -20,6 +20,8 @@ public:
private:
double Re,
CFL0;
+
+ int solverExitCode;
std::vector<bool> stagPtMvmt;
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index ce4f843ff16942ac4a8d8a0d6c0db878b580dc15..527fe3c424e4c64926cedd9f933fb1fede3dce51 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -48,6 +48,7 @@ import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+from matplotlib import pyplot as plt
import math
@@ -95,18 +96,40 @@ def main():
isolver = floD.newton(pbl)
vsolver = floD.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}
+ iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
+ coupler = viscC.Coupler(iSolverAPI, vsolver)
if rbf:
- config = {'nDim': dim, 'nSections':nSections, 'nPoints':[100, 50], 'eps':0.02, 'rbftype': 'cubic', 'smoothing': 0.0, 'degree': 1, 'neighbors': 10}
- iSolverAPI = Interpol.Interpolator(isolver, blw[0], blw[1], config)
- coupler = viscC.Coupler(iSolverAPI, vsolver)
+ fig, axs = plt.subplots(2)
+ for i in range(3):
+ aeroCoeff = coupler.Run()
+ isolver.reset()
+ axs[0].plot(aeroCoeff[:,0], 'x-')
+ axs[1].plot(aeroCoeff[:,1], 'x-')
+ axs[0].set(xlabel='iters', ylabel='$c_l$')
+ axs[1].set(xlabel='iters', ylabel='$c_d$')
+ plt.show()
else:
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ 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)
+ isolver = floD.newton(pbl)
+ vsolver = floD.initBL(Re, M_inf, CFL0, nSections)
coupler = floVC.Coupler(isolver, vsolver)
coupler.CreateProblem(blw[0], blw[1])
+ fig, axs = plt.subplots(2)
+ for i in range(1):
+ aeroCoeff = coupler.Run()
+ isolver.reset()
+ axs[0].plot(aeroCoeff[:,0], 'x-')
+ axs[1].plot(aeroCoeff[:,1], 'x-')
+ axs[0].set(xlabel='iters', ylabel='$c_l$')
+ axs[1].set(xlabel='iters', ylabel='$c_d$')
+ plt.show()
+ #coupler.Run()
#coupler.RunPolar(alphaSeq)
- coupler.Run()
tms['solver'].stop()
# extract Cp
diff --git a/dart/viscous/Master/DAirfoil.py b/dart/viscous/Master/DAirfoil.py
index 6773c30a7e2a1b5ff1fde5da2da28db1bfec85e0..a45ec4679592651401d2211bd5c421aa0efc9009 100755
--- a/dart/viscous/Master/DAirfoil.py
+++ b/dart/viscous/Master/DAirfoil.py
@@ -47,7 +47,7 @@ class Airfoil(Boundary):
return self.T0, self.H0, self.n0, self.Ct0
- def connectList(self):
+ 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
@@ -74,12 +74,12 @@ class Airfoil(Boundary):
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 it == 0:
+ 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))
+ 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
diff --git a/dart/viscous/Master/DWake.py b/dart/viscous/Master/DWake.py
index 9ec5b93069843918b59bffe6ccd6e6b96ce025d2..85a50a48db9a7455ad0276947c0dab22e07e3dc2 100755
--- a/dart/viscous/Master/DWake.py
+++ b/dart/viscous/Master/DWake.py
@@ -32,7 +32,7 @@ class Wake(Boundary):
self.n0 = 9 # wake is always turbulent
self.Ct0 = 0
- def connectList(self):
+ 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