diff --git a/dart/pyVII/vAPI.py b/dart/pyVII/vAPI.py
new file mode 100644
index 0000000000000000000000000000000000000000..59a9599bfc30cc1cd990a0c31ae80d759304a76c
--- /dev/null
+++ b/dart/pyVII/vAPI.py
@@ -0,0 +1,102 @@
+import numpy as np
+import dart.pyVII.vStruct as vStruct
+from matplotlib import pyplot as plt
+from scipy import interpolate
+
+class dartAPI:
+
+ def __init__(self, _dartSolver, bnd_airfoil, bnd_wake, _nSections):
+ self.nSections = _nSections
+ self.iSolver = _dartSolver
+ self.invStruct = [[vStruct.Airfoil(bnd_airfoil), vStruct.Wake(bnd_wake)] for _ in range(self.nSections)] # Parameters passing protocols
+ self.viscStruct = [[vStruct.viscousIPS("upper"), vStruct.viscousIPS("lower"), vStruct.viscousIPS("wake")] for _ in range(self.nSections)]
+
+ def GetInvBC(self, vSolver, couplIter):
+
+ # Extract parameters
+ for iSec in range(self.nSections):
+ dataIn = [None, None]
+ for n in range(0, len(self.invStruct[iSec])):
+ for i in range(0, len(self.invStruct[iSec][n].boundary.nodes)):
+
+ self.invStruct[iSec][n].v[i,0] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][0]
+ self.invStruct[iSec][n].v[i,1] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][1]
+ self.invStruct[iSec][n].v[i,2] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][2]
+ self.invStruct[iSec][n].Me[i] = self.iSolver.M[self.invStruct[iSec][n].boundary.nodes[i].row]
+ self.invStruct[iSec][n].rhoe[i] = self.iSolver.rho[self.invStruct[iSec][n].boundary.nodes[i].row]
+
+ self.invStruct[iSec][n].connectListNodes, self.invStruct[iSec][n].connectListElems, dataIn[n] = self.invStruct[iSec][n].connectList()
+
+ if len(dataIn[n]) == 2:
+ self.viscStruct[iSec][0].SetValues(dataIn[n][0])
+ self.viscStruct[iSec][1].SetValues(dataIn[n][1])
+ else:
+ self.viscStruct[iSec][2].SetValues(dataIn[n])
+
+ for iRegion in range(len(self.viscStruct[iSec])):
+ if couplIter==0:
+ self.viscStruct[iSec][iRegion].xPrev = self.viscStruct[iSec][iRegion].x
+ vSolver.SetGlobMesh(iSec, iRegion, self.viscStruct[iSec][iRegion].x,
+ self.viscStruct[iSec][iRegion].y,
+ self.viscStruct[iSec][iRegion].z)
+
+ vSolver.SetVelocities(iSec, iRegion, self.viscStruct[iSec][iRegion].vx,
+ self.viscStruct[iSec][iRegion].vy,
+ self.viscStruct[iSec][iRegion].vz)
+
+ vSolver.SetMe(iSec, iRegion, self.viscStruct[iSec][iRegion].Me)
+ vSolver.SetRhoe(iSec, iRegion, self.viscStruct[iSec][iRegion].Rhoe)
+ vSolver.SetxxExt(iSec, iRegion, self.viscStruct[iSec][iRegion].xxExt)
+ if len(self.viscStruct[iSec][iRegion].x) != len(self.viscStruct[iSec][iRegion].xPrev):
+ f = interpolate.interp1d(self.viscStruct[iSec][iRegion].xPrev, self.viscStruct[iSec][iRegion].DeltaStarExt, bounds_error=False, fill_value="extrapolate")
+ newDS = f(self.viscStruct[iSec][iRegion].x)
+ plt.plot(self.viscStruct[iSec][iRegion].xPrev, self.viscStruct[iSec][iRegion].DeltaStarExt, '--')
+ plt.plot(self.viscStruct[iSec][iRegion].x, newDS)
+ plt.show()
+ self.viscStruct[iSec][iRegion].DeltaStarExt = newDS
+ vSolver.SetDeltaStarExt(iSec, iRegion, self.viscStruct[iSec][iRegion].DeltaStarExt)
+ self.viscStruct[iSec][iRegion].xPrev = self.viscStruct[iSec][iRegion].x
+
+ def SetBlowingVel(self, vSolver):
+ """ Collects blowing velocities from inviscid solver, maps them to inviscid mesh and """
+
+ for iSec in range(self.nSections):
+ for iRegion in range(len(self.viscStruct[iSec])):
+ self.viscStruct[iSec][iRegion].BlowingVel = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
+ self.viscStruct[iSec][iRegion].DeltaStarExt = vSolver.ExtractDeltaStar(iSec, iRegion)
+ self.viscStruct[iSec][iRegion].xxExt = vSolver.Extractxx(iSec, iRegion)
+
+ blwVelAF = np.concatenate((self.viscStruct[iSec][0].BlowingVel, self.viscStruct[iSec][1].BlowingVel))
+ blwVelWK = self.viscStruct[iSec][2].BlowingVel
+
+ self.invStruct[iSec][0].u = blwVelAF[self.invStruct[iSec][0].connectListElems.argsort()]
+ self.invStruct[iSec][1].u = blwVelWK[self.invStruct[iSec][1].connectListElems.argsort()]
+
+ for n in range(0, len(self.invStruct[iSec])):
+ for i in range(0, self.invStruct[iSec][n].nE):
+ self.invStruct[iSec][n].boundary.setU(i, self.invStruct[iSec][n].u[i])
+
+ def RunSolver(self):
+ self.iSolver.run()
+
+ def ResetSolver(self):
+ self.iSolver.reset()
+
+ def GetCp(self,bnd):
+ import dart.utils as floU
+
+ Cp = floU.extract(bnd, self.iSolver.Cp)
+ return Cp
+
+ def GetAlpha(self):
+ return self.iSolver.pbl.alpha
+
+ def GetMinf(self):
+ return self.iSolver.pbl.M_inf
+
+ def GetCl(self):
+ return self.iSolver.Cl
+
+ def GetCm(self):
+ return self.iSolver.Cm
+
\ No newline at end of file
diff --git a/dart/vCoupler.py b/dart/pyVII/vCoupler.py
similarity index 93%
rename from dart/vCoupler.py
rename to dart/pyVII/vCoupler.py
index bb82b95736333ad634de34e6eef44a091f1da1c7..2b69eca0b9b50969d57677bf2eba9190f4b2615e 100644
--- a/dart/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -19,15 +19,13 @@
## VII Coupler
# Paul Dechamps
-from dart.viscous.Master.DAirfoil import Airfoil
-from dart.viscous.Master.DWake import Wake
from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
from fwk.coloring import ccolors
import math
import fwk
-import dart.viscUtils as viscU
+import dart.pyVII.vUtils as viscU
import numpy as np
class Coupler:
@@ -50,7 +48,7 @@ class Coupler:
couplIter = 0
cdPrev = 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))
+ .format(self.iSolverAPI.GetAlpha()*180/math.pi, self.iSolverAPI.GetMinf(), self.vSolver.GetRe()/1e6))
print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
print(' --------------------------------------------------------')
while couplIter < self.maxCouplIter:
@@ -64,7 +62,7 @@ class Coupler:
# Extract inviscid boundary
self.tms['processing'].start()
- self.__ImposeInvBC(couplIter)
+ self.iSolverAPI.GetInvBC(self.vSolver, couplIter)
self.tms['processing'].stop()
# Run viscous solver
@@ -74,7 +72,7 @@ class Coupler:
# print('Viscous ops terminated with exit code ', exitCode)
self.tms['processing'].start()
- self.__ImposeBlwVel()
+ self.iSolverAPI.SetBlowingVel(self.vSolver)
self.tms['processing'].stop()
error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
@@ -82,18 +80,19 @@ class Coupler:
if error <= self.couplTol and exitCode==0:
print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
- self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
+ self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
return 0
cdPrev = self.vSolver.Cdt
if couplIter%5==0:
- xtr=self.vSolver.Getxtr()
print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
- self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
+ 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)
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
new file mode 100644
index 0000000000000000000000000000000000000000..01518f4e15fe148f518c932ec342b0ecee15f6b1
--- /dev/null
+++ b/dart/pyVII/vInterpolator.py
@@ -0,0 +1,55 @@
+import numpy as np
+import math as m
+from matplotlib import pyplot as plt
+from scipy import interpolate
+
+class Interpolator:
+
+ def __init__(self, bnd) -> None:
+ arf, up, lw = self.GetArf(bnd)
+
+ def GetArf(self, bnd):
+ "Extract and sort the node coordinates of the airfoil."
+ arf = np.zeros((len(self.bnd.nodes),3))
+
+ # Extract coordinates.
+ for iNode, n in enumerate(bnd.nodes):
+ for iDir in range(len(arf[0,:])):
+ arf[iNode,iDir] = n.pos[iDir]
+
+ upper = np.empty((0,3))
+ lower = np.empty((0,3))
+
+ # Split between upper and lower surfaces.
+ for i in range(len(arf[:,0])):
+ if arf[i,1]>=0:
+ upper=np.vstack([upper, arf[i,:]])
+ else:
+ lower=np.vstack([lower, arf[i,:]])
+
+ # Sort according to x.
+ upper = upper[upper[:, 0].argsort()]
+ lower = lower[lower[:, 0].argsort()]
+
+ # Look for doublons.
+ for i in range(len(upper[:,0])-1):
+ if upper[i+1,0] == 0 and upper[i,0] == 1:
+ lower = np.vstack([lower, upper[i,:]])
+ upper = np.delete(upper, i, 0)
+ break
+ return arf, upper, lower
+
+
+ def meshGeneration(self):
+ print(arf[:,0])
+ x_up = np.linspace(0,1,100)
+ x_lw = np.linspace(0,1,100)
+ func_up = interpolate.interp1d(arf[:self.sep,0], arf[:self.sep,1], bounds_error=False, fill_value="extrapolate")
+ func_lw = interpolate.interp1d(arf[self.sep:,0], arf[self.sep:,1], bounds_error=False, fill_value="extrapolate")
+ y_up = func_up(x_up)
+ y_lw = func_lw(x_lw)
+
+ plt.plot(arf[:self.sep-1,0], arf[:self.sep-1,1])
+ plt.plot(x_up,y_up, 'o')
+ plt.plot(x_lw,y_lw, 'o')
+ plt.show()
\ No newline at end of file
diff --git a/dart/viiIPS.py b/dart/pyVII/vStruct.py
similarity index 94%
rename from dart/viiIPS.py
rename to dart/pyVII/vStruct.py
index 61f5d957fe6f1011ae109f723ea91e238580a2d4..65a895cbbb5bbad4ff171e7124cb0d3c7f449818 100644
--- a/dart/viiIPS.py
+++ b/dart/pyVII/vStruct.py
@@ -23,22 +23,23 @@ import numpy as np
class viscousIPS:
def __init__(self, _name) -> None:
self.name = _name
+ self.DeltaStarExt = np.zeros(1)
+ self.xxExt = np.zeros(1)
def SetValues(self, data):
- self.x = data[:,1]
- self.y = data[:,2]
- self.z = data[:,3]
- self.vx = data[:,4]
- self.vy = data[:,5]
- self.Me = data[:,6]
- self.Rhoe = data[:,7]
- self.DeltaStarExt = data[:,8]
- self.xxExt = data[:,9]
-
- self.BlowingVel = np.zeros(len(data[:,1]))
-
+ self.x = data[:,0]
+ self.y = data[:,1]
+ self.z = data[:,2]
+ self.vx = data[:,3]
+ self.vy = data[:,4]
+ self.vz = np.zeros(len(data[:,4]))
+ self.Me = data[:,5]
+ self.Rhoe = data[:,6]
+ #self.DeltaStarExt = data[:,7]
+ #self.xxExt = data[:,8]
+ self.BlowingVel = np.zeros(len(data[:,0]))
class inviscidIPS:
def __init__(self, _boundary):
@@ -52,6 +53,7 @@ class inviscidIPS:
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
class Airfoil(inviscidIPS):
def __init__(self, _boundary):
diff --git a/dart/viscUtils.py b/dart/pyVII/vUtils.py
similarity index 97%
rename from dart/viscUtils.py
rename to dart/pyVII/vUtils.py
index 74b523eff13fd159fbcb8f0559c63af6861d7552..af3336d01337783f35a432c095b454cc338475b4 100644
--- a/dart/viscUtils.py
+++ b/dart/pyVII/vUtils.py
@@ -1,7 +1,7 @@
import dart
-def initBL(Re, Minf, CFL0, meshFactor):
- solver = dart.ViscSolver(Re, Minf, CFL0, meshFactor)
+def initBL(Re, Minf, CFL0, nSections):
+ solver = dart.ViscSolver(Re, Minf, CFL0, nSections)
return solver
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index f370150c56cec068027c1e820b540ada956036d8..60b71287e7ffc320e5e44832197094d9cb826dd9 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -12,7 +12,7 @@ BLRegion::BLRegion()
{
closSolver = new Closures();
mesh = new ViscMesh();
- invBC = new InvBnd();
+ invBnd = new InvBnd();
} // end BLRegion
@@ -23,53 +23,37 @@ BLRegion::~BLRegion()
std::cout << "~BLRegion()\n";
} // end ~BLRegion
-void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
+void BLRegion::SetMesh(std::vector<double> x,
+ std::vector<double> y,
+ std::vector<double> z)
{
- unsigned int nMarkers = locM.size();
-
- mesh->SetLoc(locM);
- mesh->SetGlob(globM);
-
- /* Resize all variables */
-
- U.resize(nVar * nMarkers, 0);
- Ret.resize(nMarkers, 0);
- Cd.resize(nMarkers);
- Cf.resize(nMarkers,0);
- Hstar.resize(nMarkers,0);
- Hstar2.resize(nMarkers,0);
- Hk.resize(nMarkers,0);
- Cteq.resize(nMarkers,0);
- Us.resize(nMarkers,0);
- DeltaStar.resize(nMarkers,0);
- Delta.resize(nMarkers,0);
- Regime.resize(nMarkers,0);
-} // End SetMesh
+ unsigned int nMarkers = x.size();
-void BLRegion::SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
-{
- if (_xxExt.size() != mesh->GetnMarkers() || _deltaStarExt.size() != mesh->GetnMarkers())
- {
- throw std::runtime_error("dart::BLRegion: Miss match between mesh and external variables.");
- }
+ mesh->SetGlob(x, y, z);
- mesh->SetExt(_xxExt);
- invBC->SetDeltaStar(_deltaStarExt);
-}
+ /* Resize and reset all variables if needed */
-void BLRegion::SetInvBC(std::vector<double> _Ue,
- std::vector<double> _Me,
- std::vector<double> _Rhoe)
-{
- invBC->SetMe(_Me);
- invBC->SetVt(_Ue);
- invBC->SetRhoe(_Rhoe);
-} // End SetInvBC
+ if (mesh->GetDiffnElm()!=0)
+ {
+ U.resize(nVar * nMarkers, 0);
+ Ret.resize(nMarkers, 0);
+ Cd.resize(nMarkers);
+ Cf.resize(nMarkers,0);
+ Hstar.resize(nMarkers,0);
+ Hstar2.resize(nMarkers,0);
+ Hk.resize(nMarkers,0);
+ Cteq.resize(nMarkers,0);
+ Us.resize(nMarkers,0);
+ DeltaStar.resize(nMarkers,0);
+ Delta.resize(nMarkers,0);
+ Regime.resize(nMarkers,0);
+ }
+} // End SetMesh
void BLRegion::SetStagBC(double Re)
{
- double dv0 = (invBC->GetVt(1) - invBC->GetVt(0)) / (mesh->GetLoc(1) - mesh->GetLoc(0));
+ double dv0 = (invBnd->GetVt(1) - invBnd->GetVt(0)) / (mesh->GetLoc(1) - mesh->GetLoc(0));
if (dv0 > 0)
{
@@ -81,7 +65,7 @@ void BLRegion::SetStagBC(double Re)
} // end else
U[1] = 2.23; // H
U[2] = 0; // N
- U[3] = invBC->GetVt(0); // ue
+ U[3] = invBnd->GetVt(0); // ue
U[4] = 0; // Ctau
Ret[0] = U[3] * U[0] * Re;
@@ -96,7 +80,7 @@ void BLRegion::SetStagBC(double Re)
DeltaStar[0] = U[0] * U[1];
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBC->GetMe(0), name);
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBnd->GetMe(0), name);
Hstar[0] = ClosParam[0];
Hstar2[0] = ClosParam[1];
@@ -115,7 +99,7 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
U[0] = UpperCond[0] + LowerCond[0]; /* (Theta_up+Theta_low). */
U[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / U[0]; /* ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low). */
U[2] = 9;
- U[3] = invBC->GetVt(0); /* Inviscid velocity. */
+ U[3] = invBnd->GetVt(0); /* Inviscid velocity. */
U[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / U[0]; /* ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). */
Ret[0] = U[3] * U[0] * Re; /* Reynolds number based on the momentum thickness. */
@@ -129,7 +113,7 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
/* Laminar closures. */
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBC->GetMe(0), name);
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBnd->GetMe(0), name);
Hstar[0] = ClosParam[0];
Hstar2[0] = ClosParam[1];
@@ -146,6 +130,19 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
}
}
+void BLRegion::ComputeBlwVel()
+{
+ DeltaStar[0] = U[0] * U[1];
+ Blowing.resize(mesh->GetnMarkers()-1,0);
+ for (size_t iPoint = 1; iPoint<mesh->GetnMarkers(); ++iPoint)
+ {
+ DeltaStar[iPoint] = U[iPoint*nVar + 0] * U[iPoint*nVar + 1];
+ Blowing[iPoint-1] = (invBnd->GetRhoe(iPoint) * invBnd->GetVt(iPoint) * DeltaStar[iPoint]
+ - invBnd->GetRhoe(iPoint-1) * invBnd->GetVt(iPoint-1) * DeltaStar[iPoint-1])/(invBnd->GetRhoe(iPoint)*(mesh->GetLoc(iPoint) - mesh->GetLoc(iPoint-1)));
+ }
+}
+
+
void BLRegion::PrintSolution(size_t iPoint)
{
std::cout << "\nPoint " << iPoint << " \n"
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index ec24d0e2de44694750f13993610899118630fe2d..b81d15e6edc87537a5ef19a56987b59e1998758c 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -4,6 +4,7 @@
#include "dart.h"
#include "wClosures.h"
#include "wViscMesh.h"
+#include "wInvBnd.h"
#include <vector>
#include <string>
@@ -24,7 +25,7 @@ private:
public:
ViscMesh *mesh;
- InvBnd *invBC;
+ InvBnd *invBnd;
Closures *closSolver;
std::string name; /* Name of the region */
@@ -57,11 +58,9 @@ public:
~BLRegion();
/* Distribute BC */
- void SetMesh(std::vector<double> locM, std::vector<double> globM);
- void SetInvBC(std::vector<double> _Ue,
- std::vector<double> _Me,
- std::vector<double> _Rhoe);
- void SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt);
+ void SetMesh(std::vector<double> x,
+ std::vector<double> y,
+ std::vector<double> z);
/* Boundary conditions */
void SetStagBC(double Re);
diff --git a/dart/src/wInterpolator.h b/dart/src/wInterpolator.h
index 45a74bc2602be4aa89d680721e5eb25d62bd2e9c..95453bb05a790eb5097d294c2e1664d95137cc83 100644
--- a/dart/src/wInterpolator.h
+++ b/dart/src/wInterpolator.h
@@ -15,8 +15,8 @@ class DART_API Interpolator
ViscMesh *vmesh;
public:
-
- void GenerateMesh();
+ Interpolator();
+ ~Interpolator();
};
} // namespace dart
#endif //WINTERPOLATOR_H
\ No newline at end of file
diff --git a/dart/src/wInvBnd.cpp b/dart/src/wInvBnd.cpp
index 4fd1aceb80353075e861c3f13eefbfd46564a9e0..9b19d860a6f3b0dceb4a14a54a36dc6f656cff98 100644
--- a/dart/src/wInvBnd.cpp
+++ b/dart/src/wInvBnd.cpp
@@ -13,6 +13,9 @@
#include "wInvBnd.h"
#include <iostream>
#include <vector>
+#include <unordered_set>
+#include <iomanip>
+#include <fstream>
/* --- Namespaces -------------------------------------------------------------------- */
@@ -23,12 +26,3 @@ using namespace dart;
InvBnd::InvBnd(){}
InvBnd::~InvBnd(){}
-
-void InvBnd::ComputeBlowingVel(std::vector<double> DeltaStar, std::vector<double> LocMarkers)
-{
- BlowingVel.resize(DeltaStar.size()-1, 0);
- for (size_t iPoint = 1; iPoint < DeltaStar.size(); ++iPoint)
- {
- BlowingVel[iPoint - 1] = (Rhoe[iPoint] * Vt[iPoint] * DeltaStar[iPoint] - Rhoe[iPoint - 1] * Vt[iPoint - 1] * DeltaStar[iPoint - 1]) / (Rhoe[iPoint] * (LocMarkers[iPoint] - LocMarkers[iPoint - 1]));
- }
-}
diff --git a/dart/src/wInvBnd.h b/dart/src/wInvBnd.h
index 9a546fa9e61d89afa2d5eeab7ce0fdbb6b7e1856..77a88cbbe0ca5bebfe51c3abf7b309be82492399 100644
--- a/dart/src/wInvBnd.h
+++ b/dart/src/wInvBnd.h
@@ -4,6 +4,7 @@
#include "dart.h"
#include "wBLRegion.h"
#include <vector>
+#include <iostream>
namespace dart
{
@@ -30,11 +31,9 @@ class DART_API InvBnd
void SetMe(std::vector<double> _Me) {Me = _Me;};
void SetVt(std::vector<double> _Vt) {Vt = _Vt;};
- void SetRhoe(std::vector<double> _Rhoe) {DeltaStar = _Rhoe;};
+ void SetRhoe(std::vector<double> _Rhoe) {Rhoe = _Rhoe;};
void SetDeltaStar(std::vector<double> _DeltaStar) {DeltaStar = _DeltaStar;};
- void ComputeBlowingVel(std::vector<double> deltaStar, std::vector<double> LocMarkers);
-
};
} // namespace dart
#endif //WINVBND_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 43d837ced7956634426f6aff0a01b2c9b4210b04..a3c0c7d65271a53172ef550d5be3b1001f1d390e 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -39,7 +39,7 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
bl->U[iPoint * nVar + 0] = bl->U[(iPoint - 1) * nVar + 0];
bl->U[iPoint * nVar + 1] = bl->U[(iPoint - 1) * nVar + 1];
bl->U[iPoint * nVar + 2] = bl->U[(iPoint - 1) * nVar + 2];
- bl->U[iPoint * nVar + 3] = bl->invBC->GetVt(iPoint);
+ bl->U[iPoint * nVar + 3] = bl->invBnd->GetVt(iPoint);
bl->U[iPoint * nVar + 4] = bl->U[(iPoint - 1) * nVar + 4];
bl->Ret[iPoint] = bl->Ret[iPoint - 1];
@@ -80,7 +80,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
double CFL = CFL0;
unsigned int itFrozenJac = itFrozenJac0;
- double timeStep = SetTimeStep(CFL, dx, bl->invBC->GetVt(iPoint));
+ double timeStep = SetTimeStep(CFL, dx, bl->invBnd->GetVt(iPoint));
Matrix<double, 5, 5> IMatrix;
IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
@@ -122,7 +122,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
/* CFL adaptation. */
CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
- timeStep = SetTimeStep(CFL, dx, bl->invBC->GetVt(iPoint));
+ timeStep = SetTimeStep(CFL, dx, bl->invBnd->GetVt(iPoint));
IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
innerIters++;
@@ -176,7 +176,7 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
if (bl->Regime[iPoint] == 0)
{
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->Ret[iPoint], bl->invBC->GetMe(iPoint), bl->name);
+ std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -188,7 +188,7 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
} // End if
else if (bl->Regime[iPoint] == 1)
{
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->U[iPoint * nVar + 4], bl->Ret[iPoint], bl->invBC->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[iPoint * nVar + 0], bl->U[iPoint * nVar + 1], bl->U[iPoint * nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 8c02dd5d2d92ae97170029fae39e9ad210f4415a..709040486c8ce4ef75996f19bb684b3e3f9af11c 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -83,7 +83,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
if (bl->Regime[iPoint] == 0)
{
/* Laminar closure relations */
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[1], bl->Ret[iPoint], bl->invBC->GetMe(iPoint), bl->name);
+ std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[1], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
@@ -97,7 +97,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
else if (bl->Regime[iPoint] == 1)
{
/* Turbulent closure relations */
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(u[0], u[1], u[4], bl->Ret[iPoint], bl->invBC->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam = bl->closSolver->TurbulentClosures(u[0], u[1], u[4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
@@ -122,8 +122,8 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
double dCt_dx = (u[4] - bl->U[(iPoint - 1) * nVar + 4]) / dx;
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
- double dueExt_dx = (bl->invBC->GetVt(iPoint) - bl->invBC->GetVt(iPoint - 1)) / dxExt;
- double ddeltaStarExt_dx = (bl->invBC->GetDeltaStar(iPoint) - bl->invBC->GetDeltaStar(iPoint - 1)) / dxExt;
+ double dueExt_dx = (bl->invBnd->GetVt(iPoint) - bl->invBnd->GetVt(iPoint - 1)) / dxExt;
+ double ddeltaStarExt_dx = (bl->invBnd->GetDeltaStar(iPoint) - bl->invBnd->GetDeltaStar(iPoint - 1)) / dxExt;
double c = 4 / (M_PI * dx);
double cExt = 4 / (M_PI * dxExt);
@@ -144,7 +144,7 @@ Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<do
dN_dx = 0;
} // End else
- double Mea = bl->invBC->GetMe(iPoint);
+ double Mea = bl->invBnd->GetMe(iPoint);
double Hstara = bl->Hstar[iPoint];
double Hstar2a = bl->Hstar2[iPoint];
diff --git a/dart/src/wViscMesh.cpp b/dart/src/wViscMesh.cpp
index b1476011d42d70a1fc27fbf15ae61cbf871978ed..69546187ae84ef7da3753f76de5761acb72cbc4e 100644
--- a/dart/src/wViscMesh.cpp
+++ b/dart/src/wViscMesh.cpp
@@ -21,31 +21,24 @@ using namespace dart;
/* ----------------------------------------------------------------------------------- */
-ViscMesh::ViscMesh(){}
+ViscMesh::ViscMesh(){prevnMarkers = 0;}
ViscMesh::~ViscMesh(){}
-void ViscMesh::SetGlob(std::vector<double> x, std::vector<double> y, std::vector<double> z)
+void ViscMesh::SetGlob(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
{
if (x.size() != y.size() || y.size() != z.size() || x.size() != z.size())
{
throw std::runtime_error("dart::ViscMesh Wrong mesh size.");
}
- nMarkers = x.size();
+ nMarkers = _x.size();
nElm = nMarkers - 1;
- Glob.resize(nMarkers);
- for (size_t iPoint=0; iPoint<Glob.size(); ++iPoint)
- {
- Glob[iPoint].resize(3, 0);
- }
+ x = _x;
+ y = _y;
+ z = _z;
- for (size_t iPoint=0; iPoint<x.size(); ++iPoint)
- {
- Glob[iPoint][0] = x[iPoint];
- Glob[iPoint][1] = y[iPoint];
- Glob[iPoint][2] = z[iPoint];
- }
+ ComputeBLcoord();
}
void ViscMesh::ComputeBLcoord()
@@ -56,9 +49,9 @@ void ViscMesh::ComputeBLcoord()
for (size_t iPoint=0; iPoint<nElm; ++iPoint)
{
- alpha[iPoint] = std::atan2(Glob[iPoint+1][1], Glob[iPoint+1][0]);
+ alpha[iPoint] = std::atan2(y[iPoint+1]-y[iPoint], x[iPoint+1]-x[iPoint]-1);
/* dx = sqrt(Delta x ^2 + Delta y ^2) */
- dx[iPoint] = std::sqrt((Glob[iPoint+1][0]-Glob[iPoint][0])*(Glob[iPoint+1][0]-Glob[iPoint][0]) + (Glob[iPoint+1][1]-Glob[iPoint][1])*(Glob[iPoint+1][1]-Glob[iPoint][1]));
+ dx[iPoint] = std::sqrt((x[iPoint+1]-x[iPoint])*(x[iPoint+1]-x[iPoint]) + (y[iPoint+1]-y[iPoint])*(y[iPoint+1]-y[iPoint]));
Loc[iPoint + 1] = Loc[iPoint] + dx[iPoint];
}
}
diff --git a/dart/src/wViscMesh.h b/dart/src/wViscMesh.h
index 1497faa25ef18041e7799566685f30f7ec366738..decfbe845e2cc54877b9bfd4b5a982c1d5f14be3 100644
--- a/dart/src/wViscMesh.h
+++ b/dart/src/wViscMesh.h
@@ -4,6 +4,7 @@
#include "dart.h"
#include "wBLRegion.h"
#include <vector>
+#include <iostream>
namespace dart
{
@@ -11,8 +12,10 @@ namespace dart
class DART_API ViscMesh
{
private:
- std::vector<std::vector<double>> Glob;
- std::vector<double> Loc, /* BL (local) coordinates */
+ std::vector<double> x,
+ y,
+ z,
+ Loc, /* BL (local) coordinates */
Ext, /* Coordinates at the edge of the BL */
dx,
alpha; /* Angle of the element for the rotation matrix */
@@ -20,18 +23,23 @@ class DART_API ViscMesh
unsigned int nMarkers;
unsigned int nElm;
+
public:
+ unsigned int prevnMarkers;
+
ViscMesh();
~ViscMesh();
unsigned int GetnMarkers() const {return nMarkers;};
double GetLoc(size_t iPoint) const {return Loc[iPoint];};
- double Getx(size_t iPoint) const {return Glob[iPoint][0];};
- double Gety(size_t iPoint) const {return Glob[iPoint][1];};
- double Getz(size_t iPoint) const {return Glob[iPoint][2];};
+ double Getx(size_t iPoint) const {return x[iPoint];};
+ double Gety(size_t iPoint) const {return y[iPoint];};
+ double Getz(size_t iPoint) const {return z[iPoint];};
double GetExt(size_t iPoint) const {return Ext[iPoint];};
- void SetGlob(std::vector<double> x, std::vector<double> y, std::vector<double> z);
+ unsigned int GetDiffnElm() const {return nMarkers-prevnMarkers;};
+
+ void SetGlob(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
void SetExt(std::vector<double> ExtM) {Ext = ExtM;};
void ComputeBLcoord();
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 3028c974e4b9c56ca4d89d5373c1912bfef8b2d3..e2d6e91a381a0dc6e59bbe00b44193f9dbfd04b6 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -13,9 +13,10 @@
#include "wViscSolver.h"
#include "wBLRegion.h"
#include "wTimeSolver.h"
-#include "wSolutionInterp.h"
+#include "wInterpolator.h"
#include <iostream>
#include <vector>
+#include <cmath>
/* --- Namespaces -------------------------------------------------------------------- */
@@ -51,15 +52,12 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSec
Sections[iSec][2]->name = "wake";
}
+ stagPtMvmt.resize(Sections.size(), false);
+
/* Temporal solver initialization */
tSolver = new TimeSolver(_CFL0, _Minf, Re);
- /* Stagnation point movements monitoring and control */
-
- nbElmUpper.resize(Sections.size(), 0);
- stagPtMvmt.resize(Sections.size(), false);
-
/* User input numerical parameters */
std::cout << "--- Viscous solver setup ---\n"
@@ -86,76 +84,11 @@ ViscSolver::~ViscSolver()
}
delete tSolver;
- delete meshConverter;
-
PrintBanner();
std::cout << "~ViscSolver()\n";
} // end ~ViscSolver
-/**
- * @brief Initialize 1D mesh of specified region.
- */
-void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM)
-{
- bl[region]->coarseSize = locM.size();
- if (mapMesh)
- {
- std::vector<double> fLocMesh = meshConverter->CreateFineMesh(locM);
- std::vector<double> fGlobMesh = meshConverter->CreateFineMesh(globM);
- bl[region]->fineSize = fLocMesh.size();
- bl[region]->SetMesh(fLocMesh, fGlobMesh);
- //std::cout << "- " << bl[region]->name << ": Mesh mapped from " << locM.size() << " to " << bl[region]->GetnMarkers() << " elements." << std::endl;
- }
- else
- {
- bl[region]->SetMesh(locM, globM);
- bl[region]->fineSize = locM.size();
- //std::cout << "- " << bl[region]->name << ": Mesh elements " << bl[region]->GetnMarkers() << std::endl;
- }
-}
-
-/**
- * @brief Set inviscid boundary conditions of the coupled problem.
- */
-void ViscSolver::SendInvBC(size_t region, std::vector<double> _Ue,
- std::vector<double> _Me,
- std::vector<double> _Rhoe)
-{
-
- if (mapMesh)
- {
- std::vector<double> UeFine = meshConverter->InterpolateToFineMesh(_Ue, bl[region]->fineSize);
- std::vector<double> MeFine = meshConverter->InterpolateToFineMesh(_Me, bl[region]->fineSize);
- std::vector<double> RhoeFine = meshConverter->InterpolateToFineMesh(_Rhoe, bl[region]->fineSize);
-
- bl[region]->SetInvBC(UeFine, MeFine, RhoeFine);
- //std::cout << "- " << bl[region]->name << ": Inviscid boundary interpolated." << std::endl;
- }
- else
- {
- bl[region]->SetInvBC(_Ue, _Me, _Rhoe);
- }
-}
-
-/**
- * @brief Set mesh and displacement thickness at the edge of the boundary layer.
- * Used for the interaction law.
- */
-void ViscSolver::SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
-{
- if (mapMesh)
- {
- std::vector<double> xxExtFine = meshConverter->InterpolateToFineMesh(_xxExt, bl[region]->fineSize);
- std::vector<double> deltaStarExtFine = meshConverter->InterpolateToFineMesh(_deltaStarExt, bl[region]->fineSize);
- bl[region]->SetExtVariables(xxExtFine, deltaStarExtFine);
- }
- else
- {
- bl[region]->SetExtVariables(_xxExt, _deltaStarExt);
- }
-}
-
/**
* @brief Runs viscous operations. Temporal solver parametrisation is first adapted,
* Solutions are initialized than time marched towards steady state successively for each points.
@@ -171,11 +104,42 @@ int ViscSolver::Run(unsigned int couplIter)
/* Prepare time integration */
- if ((couplIter == 0) || Sections[iSec][0]->mesh->GetDiffnElm() || (stagPtMvmt[iSec]))
+ bool resetSolution = false;
+ if ((couplIter == 0) ||
+ Sections[iSec][0]->mesh->GetDiffnElm()!=0 ||
+ Sections[iSec][1]->mesh->GetDiffnElm()!=0 ||
+ Sections[iSec][2]->mesh->GetDiffnElm()!=0 ||
+ (stagPtMvmt[iSec]))
+ {
+ resetSolution = true;
+ }
+ else
+ {
+ resetSolution = false;
+ }
+ /*std::cout << "resetSolution " << resetSolution << std::endl;
+
+ std::cout << Sections[iSec][0]->mesh->GetDiffnElm() <<
+ Sections[iSec][1]->mesh->GetDiffnElm() <<
+ Sections[iSec][2]->mesh->GetDiffnElm() << std::endl;*/
+
+ if (resetSolution)
{
tSolver->SetCFL0(1);
tSolver->SetitFrozenJac(1);
tSolver->SetinitSln(true);
+
+ /* Special case where the external flow mesh cannot be used for the interaction law */
+ for (size_t iRegion=0; iRegion<Sections[iSec].size(); ++iRegion)
+ {
+ std::vector<double> xxExtCurr;
+ for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
+ {
+ xxExtCurr.push_back(Sections[iSec][iRegion]->mesh->GetLoc(iPoint));
+ }
+ Sections[iSec][iRegion]->mesh->SetExt(xxExtCurr);
+ }
+
if (couplIter != 0 && (Sections[iSec][0]->mesh->GetDiffnElm()))
{
stagPtMvmt[iSec] = true;
@@ -205,11 +169,11 @@ int ViscSolver::Run(unsigned int couplIter)
std::cout << "Updating solution" << std::endl;
}
}
-
- for (size_t iSec=0; iSec<Sections.size(); ++iSec)
- {
- nbElmUpper[iSec] = Sections[iSec][0]->GetnMarkers();
+ for (size_t iRegion=0; iRegion<Sections[iSec].size(); ++iRegion)
+ {
+ Sections[iSec][iRegion]->mesh->prevnMarkers = Sections[iSec][iRegion]->mesh->GetnMarkers();
+ }
/* Set boundary condition */
@@ -220,11 +184,16 @@ int ViscSolver::Run(unsigned int couplIter)
Sections[iSec][0]->SetStagBC(Re);
Sections[iSec][1]->SetStagBC(Re);
- /* Loop over the regions */
+ /* Loop over the regions */
for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
{
-
+ if (printOn)
+ {
+ std::cout << "- " << Sections[iSec][iRegion]->name << " region";
+ }
+ lockTrans = false;
+
/* Reset regime */
if (iRegion < 2)
{
@@ -239,33 +208,17 @@ int ViscSolver::Run(unsigned int couplIter)
{
Sections[iSec][iRegion]->Regime[k] = 1;
}
- }
- lockTrans = false;
-
- if (printOn)
- {
- std::cout << "- " << Sections[iSec][iRegion]->name << " region";
- }
- if (iRegion == 2)
- {
/* Impose wake boundary condition */
std::vector<double> UpperCond;
std::vector<double> LowerCond;
UpperCond = std::vector<double>(Sections[iSec][0]->U.end() - Sections[iSec][0]->GetnVar(), Sections[iSec][0]->U.end()); /* Base std::vector constructor to slice */
LowerCond = std::vector<double>(Sections[iSec][1]->U.end() - Sections[iSec][1]->GetnVar(), Sections[iSec][1]->U.end()); /* Base std::vector constructor to slice */
Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
- if (printOn)
- {
- std::cout << "\n"
- << std::endl;
- }
-
lockTrans = true;
}
/* Loop over points */
-
for (size_t iPoint = 1; iPoint < Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
{
@@ -274,6 +227,15 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
/* Time integration */
+ /*std::cout << " iSec " << iSec
+ << " iPoint " << iPoint
+ << " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
+ << " Vt " << Sections[iSec][iRegion]->invBnd->GetVt(iPoint)
+ << " Me " << Sections[iSec][iRegion]->invBnd->GetMe(iPoint)
+ << " xxExt " << Sections[iSec][iRegion]->mesh->GetExt(iPoint)
+ << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
+ << std::endl;
+ Sections[iSec][iRegion]->PrintSolution(iPoint);*/
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
@@ -308,14 +270,16 @@ int ViscSolver::Run(unsigned int couplIter)
} // End if (!lockTrans).
} // End for iPoints
+ Sections[iSec][iRegion]->ComputeBlwVel();
} // End for iRegion
- }
- ComputeDrag();
- for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
+ } // End for iSections
+
+ ComputeDrag(Sections[0]);
+ if (printOn)
{
- Sections[iSec][iRegion]->ComputeBlwVel();
- } // End for iRegion
+ std::cout << "\n" << std::endl;
+ }
return solverExitCode; // exit code
}
@@ -358,20 +322,20 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->transMarker = iPoint; /* Save transition marker */
/* Loop over remaining points in the region */
- for (size_t k = 0; k < bl->mesh->GernMarkers() - iPoint; ++k)
+ for (size_t k = 0; k < bl->mesh->GetnMarkers() - iPoint; ++k)
{
bl->Regime[iPoint + k] = 1; // Set Turbulent regime.
}
/* Compute transition location */
- bl->xtr = (bl->GetnCrit() - bl->U[(iPoint - 1) * nVar + 2]) * ((bl->mesh->GetGlob(iPoint) - bl->mesh->GetGlob(iPoint - 1)) / (bl->U[iPoint * nVar + 2] - bl->U[(iPoint - 1) * nVar + 2])) + bl->mesh->GetGlob(iPoint - 1);
+ bl->xtr = (bl->GetnCrit() - bl->U[(iPoint - 1) * nVar + 2]) * ((bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1)) / (bl->U[iPoint * nVar + 2] - bl->U[(iPoint - 1) * nVar + 2])) + bl->mesh->Getx(iPoint - 1);
/* Percentage of the subsection corresponding to a turbulent flow */
- double avTurb = (bl->mesh->GetGlob(iPoint) - bl->xtr) / (bl->mesh->GetGlob(iPoint) - bl->mesh->GetGlob(iPoint - 1));
+ double avTurb = (bl->mesh->Getx(iPoint) - bl->xtr) / (bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1));
double avLam = 1 - avTurb;
/* Impose boundary condition */
- std::vector<double> turbParam1 = bl->closSolver->TurbulentClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->GetMe(iPoint - 1), bl->name);
+ std::vector<double> turbParam1 = bl->closSolver->TurbulentClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], 0.03, bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
bl->Cteq[iPoint - 1] = turbParam1[6];
bl->U[(iPoint - 1) * nVar + 4] = 0.7 * bl->Cteq[iPoint - 1];
@@ -396,7 +360,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->U[iPoint * nVar + 4] = avTurb * bl->U[(iPoint)*nVar + 4]; /* Ct */
/* Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations). */
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->GetMe(iPoint - 1), bl->name);
+ std::vector<double> lamParam = bl->closSolver->LaminarClosures(bl->U[(iPoint - 1) * nVar + 0], bl->U[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->invBnd->GetMe(iPoint - 1), bl->name);
bl->Hstar[iPoint - 1] = lamParam[0];
bl->Hstar2[iPoint - 1] = lamParam[1];
bl->Hk[iPoint - 1] = lamParam[2];
@@ -406,7 +370,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->Cteq[iPoint - 1] = lamParam[6];
bl->Us[iPoint - 1] = lamParam[7];
- std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[(iPoint)*nVar + 0], bl->U[(iPoint)*nVar + 1], bl->U[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->GetMe(iPoint), bl->name);
+ std::vector<double> turbParam = bl->closSolver->TurbulentClosures(bl->U[(iPoint)*nVar + 0], bl->U[(iPoint)*nVar + 1], bl->U[(iPoint)*nVar + 4], bl->Ret[iPoint], bl->invBnd->GetMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
@@ -417,66 +381,6 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->Us[iPoint] = turbParam[7];
}
-/**
- * @brief Extracts blowing velocity distribution of the corresponding region.
- */
-std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
-{
- std::vector<double> cBlwVel;
- if (mapMesh)
- {
- cBlwVel = meshConverter->BlowingFineToCoarse(bl[iRegion]->Blowing, bl[iRegion]->coarseSize - 1);
- }
- else
- {
- cBlwVel = bl[iRegion]->Blowing;
- }
- return cBlwVel;
-}
-
-/**
- * @brief Extracts displacement thickness distribution of the corresponding region.
- */
-std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
-{
- std::vector<double> cDStar;
- if (mapMesh)
- {
- // std::cout << "DeltaStar size" << bl[iRegion]->DeltaStar.size() << std::endl;
- cDStar = meshConverter->VarsFineToCoarse(bl[iRegion]->DeltaStar, bl[iRegion]->coarseSize);
- }
- else
- {
- cDStar = bl[iRegion]->DeltaStar;
- }
- return cDStar;
-}
-
-/**
- * @brief Extracts mesh distribution of the corresponding region.
- *
- * @param iRegion Region ID.
- * @return LocMarkersOut Local markers of the BL mesh.
- */
-std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
-{
- std::vector<double> LocMarkersOut;
- for (size_t k = 0; k < bl[iRegion]->GetnMarkers(); ++k)
- {
- LocMarkersOut.push_back(bl[iRegion]->GetLocMarkers(k));
- }
- std::vector<double> cXxExt;
- if (mapMesh)
- {
- cXxExt = meshConverter->VarsFineToCoarse(LocMarkersOut, bl[iRegion]->coarseSize);
- }
- else
- {
- cXxExt = LocMarkersOut;
- }
- return cXxExt;
-}
-
/**
* @brief Friction, pressure and total drag computation of the configuration.
*
@@ -484,21 +388,21 @@ std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
* Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
* Cdp = Cdtot - Cdf.
*/
-void ViscSolver::ComputeDrag()
+void ViscSolver::ComputeDrag(std::vector<BLRegion *> bl)
{
unsigned int nVar = bl[0]->GetnVar();
- unsigned int lastWkPt = (bl[2]->GetnMarkers() - 1) * nVar;
+ unsigned int lastWkPt = (bl[2]->mesh->GetnMarkers() - 1) * nVar;
/* Normalize friction and dissipation coefficients. */
std::vector<double> frictioncoeff_upper;
- for (size_t k = 0; k < bl[0]->GetnMarkers(); ++k)
+ for (size_t k = 0; k < bl[0]->mesh->GetnMarkers(); ++k)
{
frictioncoeff_upper.push_back(bl[0]->U[k * nVar + 3] * bl[0]->U[k * nVar + 3] * bl[0]->Cf[k]);
}
std::vector<double> frictioncoeff_lower;
- for (size_t k = 0; k < bl[1]->GetnMarkers(); ++k)
+ for (size_t k = 0; k < bl[1]->mesh->GetnMarkers(); ++k)
{
frictioncoeff_lower.push_back(bl[1]->U[k * nVar + 3] * bl[1]->U[k * nVar + 3] * bl[1]->Cf[k]);
}
@@ -506,14 +410,14 @@ void ViscSolver::ComputeDrag()
/* Compute friction drag coefficient. */
double Cdf_upper = 0.0;
- for (size_t k = 1; k < bl[0]->GetnMarkers(); ++k)
+ for (size_t k = 1; k < bl[0]->mesh->GetnMarkers(); ++k)
{
- Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k - 1)) / 2;
+ Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->mesh->Getx(k) - bl[0]->mesh->Getx(k - 1)) / 2;
}
double Cdf_lower = 0.0;
- for (size_t k = 1; k < bl[1]->GetnMarkers(); ++k)
+ for (size_t k = 1; k < bl[1]->mesh->GetnMarkers(); ++k)
{
- Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k - 1)) / 2;
+ Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->mesh->Getx(k) - bl[1]->mesh->Getx(k - 1)) / 2;
}
Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
@@ -527,45 +431,66 @@ void ViscSolver::ComputeDrag()
Cdp = Cdt - Cdf;
}
-std::vector<std::vector<double>> ViscSolver::ExtractSolution()
+void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z)
{
- // [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
- std::vector<std::vector<double>> blVectors(18);
- unsigned int nVar = bl[0]->GetnVar();
+ std::cout << "len x " << _x.size() << std::endl;
+ Sections[iSec][iRegion]->SetMesh(_x, _y, _z);
+}
+
+void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
+{
+ std::vector<double> Vt(vx.size());
+ double alpha=0;
- for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion)
+ for (size_t iPoint=0; iPoint<vx.size()-1; ++iPoint)
{
- for (size_t iPoint = 0; iPoint < bl[iRegion]->GetnMarkers(); ++iPoint)
- {
- blVectors[0].push_back(bl[iRegion]->GetLocMarkers(iPoint));
- blVectors[1].push_back(bl[iRegion]->GetGlobMarkers(iPoint));
- blVectors[2].push_back(bl[iRegion]->DeltaStar[iPoint]); // DeltaStar
- blVectors[3].push_back(bl[iRegion]->U[iPoint * nVar + 0]); // Theta
- blVectors[4].push_back(bl[iRegion]->U[iPoint * nVar + 1]); // H
- blVectors[5].push_back(bl[iRegion]->U[iPoint * nVar + 2]); // N
- blVectors[6].push_back(bl[iRegion]->U[iPoint * nVar + 3]); // Ue
- blVectors[7].push_back(bl[iRegion]->U[iPoint * nVar + 4]); // Ct
- blVectors[8].push_back(bl[iRegion]->Hk[iPoint]); // Hk
- blVectors[9].push_back(bl[iRegion]->Hstar[iPoint]); // H*
- blVectors[10].push_back(bl[iRegion]->Hstar2[iPoint]); // H**
- blVectors[11].push_back(bl[iRegion]->Cf[iPoint]); // Cf
- blVectors[12].push_back(bl[iRegion]->Cd[iPoint]); // Cd
- blVectors[13].push_back(bl[iRegion]->Delta[iPoint]); // delta
- blVectors[14].push_back(bl[iRegion]->GetVtExt(iPoint)); // VtExt
- blVectors[15].push_back(bl[iRegion]->GetMe(iPoint)); // Me
- blVectors[16].push_back(bl[iRegion]->GetRhoe(iPoint)); // Rhoe
- blVectors[17].push_back(bl[iRegion]->Blowing[iPoint]); // Blowing velocity
- }
+ alpha = std::atan2(Sections[iSec][iRegion]->mesh->Gety(iPoint+1) - Sections[iSec][iRegion]->mesh->Gety(iPoint),
+ Sections[iSec][iRegion]->mesh->Getx(iPoint+1) - Sections[iSec][iRegion]->mesh->Getx(iPoint));
+ Vt[iPoint] = vx[iPoint]*std::cos(alpha) + vy[iPoint]*std::sin(alpha);
}
- return blVectors;
+ Vt.back() = vx.back()*cos(alpha) + vy.back()*std::sin(alpha);
+
+ Sections[iSec][iRegion]->invBnd->SetVt(Vt);
+}
+
+void ViscSolver::SetMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
+{
+ Sections[iSec][iRegion]->invBnd->SetMe(_Me);
+}
+
+void ViscSolver::SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
+{
+ Sections[iSec][iRegion]->invBnd->SetRhoe(_Rhoe);
+}
+
+void ViscSolver::SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
+{
+ std::cout << "len xxExt " << _xxExt.size() << std::endl;
+ Sections[iSec][iRegion]->mesh->SetExt(_xxExt);
+}
+
+void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
+{
+ std::cout << "len DeltaStarExt " << _DeltaStarExt.size() << std::endl;
+ Sections[iSec][iRegion]->invBnd->SetDeltaStar(_DeltaStarExt);
}
-std::vector<double> ViscSolver::Getxtr()
+std::vector<double> ViscSolver::ExtractBlowingVel(size_t iSec, size_t iRegion)
+{
+ return Sections[iSec][iRegion]->Blowing;
+}
+std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion)
+{
+ std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0);
+ for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
+ {
+ xx[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
+ }
+ return xx;
+}
+std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion)
{
- std::vector<double> xtr(2);
- xtr[0] = bl[0]->xtr;
- xtr[1] = bl[1]->xtr;
- return xtr;
+ return Sections[iSec][iRegion]->DeltaStar;
}
void ViscSolver::PrintBanner()
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index 42a326ab048ecb15484ec22766dc710ad865b25c..d9486092be26593449bbfb1dfddbdb1d5b9fcbdd 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -3,7 +3,6 @@
#include "dart.h"
#include "wBLRegion.h"
-
namespace dart
{
@@ -21,35 +20,34 @@ public:
private:
double Re,
CFL0;
- bool mapMesh = false;
std::vector<bool> stagPtMvmt;
- std::vector<unsigned int> nbElmUpper; /* Nb of elements on the upper side. Used for stagnation point movement control. */
TimeSolver *tSolver;
- SolutionInterp *meshConverter;
public:
- ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor);
+ ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections);
~ViscSolver();
-
- void InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM);
- void SendInvBC(size_t region, std::vector<double> _Ue,
- std::vector<double> _Me,
- std::vector<double> _Rhoe);
- void SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt);
-
int Run(unsigned int couplIter);
- std::vector<double> ExtractBlowingVel(size_t iRegion);
- std::vector<double> ExtractDeltaStar(size_t iRegion);
- std::vector<double> ExtractXx(size_t iRegion);
- std::vector<std::vector<double>> ExtractSolution();
- std::vector<double> Getxtr();
+ 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);
+ void SetVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz);
+ void SetMe(size_t iSec, size_t iRegion, std::vector<double> _Me);
+ void SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe);
+ void SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt);
+ void SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt);
+
+ std::vector<double> ExtractBlowingVel(size_t iSec, size_t iRegion);
+ std::vector<double> Extractxx(size_t iSec, size_t iRegion);
+ std::vector<double> ExtractDeltaStar(size_t iSec, size_t iRegion);
+
+ double Getxtr(size_t iSec, size_t iRegion) const {return Sections[iSec][iRegion]->xtr;};
private:
void CheckWaves(size_t iPoint, BLRegion *bl);
void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
- void ComputeDrag();
+ void ComputeDrag(std::vector<BLRegion *> bl);
void ComputeBlwVel();
void PrintBanner();
};
diff --git a/dart/tests/bli_Polar.py b/dart/tests/bli_Polar.py
index f5907b4567ed1d0728cecbe988382511c68ce379..3e08f1380f598bfee417e25bd366612cb4bfb3a0 100644
--- a/dart/tests/bli_Polar.py
+++ b/dart/tests/bli_Polar.py
@@ -36,9 +36,9 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscUtils as viscU
+import dart.VII.viscUtils as viscU
-import dart.vCoupler as floVC
+import dart.VII.vCoupler as floVC
import numpy as np
import tbox
diff --git a/dart/tests/bli_Sep.py b/dart/tests/bli_Sep.py
index 85b9092ae60916604d5c7af7442d60b5399021ce..363809b55869ba500b58863197d547854a3def0f 100644
--- a/dart/tests/bli_Sep.py
+++ b/dart/tests/bli_Sep.py
@@ -36,9 +36,9 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscUtils as viscU
+import dart.VII.viscUtils as viscU
-import dart.vCoupler as floVC
+import dart.VII.vCoupler as floVC
import numpy as np
import tbox
diff --git a/dart/tests/bli_Trans.py b/dart/tests/bli_Trans.py
index 21ceb6fd9340b770b7c80f2c02d8b6248b502c8a..9e27a0a12ea13501c48c16e1432b01297f9c490a 100644
--- a/dart/tests/bli_Trans.py
+++ b/dart/tests/bli_Trans.py
@@ -36,9 +36,9 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscUtils as viscU
+import dart.VII.viscUtils as viscU
-import dart.vCoupler as floVC
+import dart.VII.vCoupler as floVC
import numpy as np
import tbox
diff --git a/dart/tests/bli_dualMesh.py b/dart/tests/bli_dualMesh.py
index 20f135802802efe5e4e0ea3983ffcf80214ad65b..ac96d71565c5808a2b2b82bee2a32e4419ed77b1 100644
--- a/dart/tests/bli_dualMesh.py
+++ b/dart/tests/bli_dualMesh.py
@@ -36,9 +36,9 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscUtils as viscU
+import dart.VII.viscUtils as viscU
-import dart.vCoupler as floVC
+import dart.VII.vCoupler as floVC
import numpy as np
import tbox
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index b250e08a9a9cfb29130ce6fa5e36eb8c3cc0bdc7..1494838939216721161c747968ab7c32bcfdc4aa 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -36,10 +36,12 @@
import dart.utils as floU
import dart.default as floD
-import dart.viscUtils as viscU
-import dart.viscAPI as viscAPI
-import dart.vCoupler as floVC
+import dart.pyVII.vUtils as viscU
+import dart.pyVII.vAPI as viscAPI
+import dart.pyVII.vCoupler as viscC
+import dart.pyVII.vInterpolator as vInterpol
+
import numpy as np
import tbox
@@ -58,7 +60,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5*math.pi/180
+ alpha = 10.*math.pi/180
#alphaSeq = np.arange(-5, 5.5, 0.5)
M_inf = 0.
@@ -73,14 +75,14 @@ def main():
dim = 2
if dim == 2:
- nSection = 1
+ nSections = 1
else:
- nSection = 10
+ nSections = 10
# 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}
+ 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()
@@ -94,23 +96,23 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
- isolverAPI = viscAPI.dartAPI(floD.newton(pbl), blw[0], blw[1])
- coupler = floVC.Coupler(isolverAPI, vsolver)
-
- coupler.CreateProblem(blw[0], blw[1])
+ vInterp = vInterpol.Interpolator(bnd)
+ quit()
+ vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
+ iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), blw[0], blw[1], nSections)
+ coupler = viscC.Coupler(iSolverAPI, vSolver)
#coupler.RunPolar(alphaSeq)
coupler.Run()
tms['solver'].stop()
# extract Cp
- Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ Cp = iSolverAPI.GetCp(bnd.groups[0].tag.elems)
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))
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, iSolverAPI.GetCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.GetCm()))
# display timers
tms['total'].stop()
@@ -120,13 +122,11 @@ def main():
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)
+ 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)
+ blScalars, blVectors = viscU.ExtractVars(vSolver)
wData = viscU.Dictionary(blScalars, blVectors, xtr)
viscU.PlotVars(wData, plotVar)"""
@@ -135,19 +135,19 @@ def main():
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))
+ tests.add(CTest('Cl', iSolverAPI.GetCl(), 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', vSolver.Getxtr(0,0), 0.20, 0.03, forceabs=True))
+ tests.add(CTest('xtr_bot', vSolver.Getxtr(0,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.add(CTest('Cl', iSolverAPI.GetCl(), 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', vSolver.Getxtr(0,0), 0.06, 0.03, forceabs=True))
+ tests.add(CTest('xtr_bot', vSolver.Getxtr(0,1), 0.74, 0.03, forceabs=True))
+ """else:
+ raise Exception('Test not defined for this flow')"""
tests.run()
diff --git a/dart/viscAPI.py b/dart/viscAPI.py
deleted file mode 100644
index 12a99ea795d4a2c186dbcb95f21ee8c1de84a41b..0000000000000000000000000000000000000000
--- a/dart/viscAPI.py
+++ /dev/null
@@ -1,66 +0,0 @@
-import numpy as np
-import dart.viiIPS as viiIPS
-
-class dartAPI:
-
- def __init__(self, _dartSolver, bnd_airfoil, bnd_wake, _nSections):
- self.nSections = _nSections
- self.dartSolver = _dartSolver
- self.parampp = [[viiIPS.Airfoil(bnd_airfoil), viiIPS.Wake(bnd_wake)] for _ in range(self.nSections)] # Parameters passing protocols
- self.regions = [[viiIPS.Region("upper"), viiIPS.Region("lower"), viiIPS.Region("wake")] for _ in range(self.nSection)]
-
- def ExtractInvBC(self):
-
- # Extract parameters
- for iSec in range(self.nSections):
- dataIn = [None, None]
- for n in range(0, len(self.parampp)):
- for i in range(0, len(self.parampp[n].boundary.nodes)):
-
- self.parampp[iSec][n].v[i,0] = self.iSolver.U[self.parampp[iSec][n].boundary.nodes[i].row][0]
- self.parampp[iSec][n].v[i,1] = self.iSolver.U[self.parampp[iSec][n].boundary.nodes[i].row][1]
- self.parampp[iSec][n].v[i,2] = self.iSolver.U[self.parampp[iSec][n].boundary.nodes[i].row][2]
- self.parampp[iSec][n].Me[i] = self.iSolver.M[self.parampp[iSec][n].boundary.nodes[i].row]
- self.parampp[iSec][n].rhoe[i] = self.iSolver.rho[self.parampp[iSec][n].boundary.nodes[i].row]
-
- self.parampp[iSec][n].connectListNodes, self.parampp[iSec][n].connectListElems, dataIn[n] = self.parampp[iSec][n].connectList()
- if len(self.parampp[iSec][n] == 2):
- self.regions[iSec][0].SetValues(dataIn[n][0])
- self.regions[iSec][1].SetValues(dataIn[n][1])
- else:
- self.region[iSec][n].SetValues(dataIn[n])
-
- def ImposeBlwVel(self, vSolver):
- """ Collects blowing velocities from inviscid solver, maps them to inviscid mesh and """
-
- for iSec in range(self.nSections):
- for iRegion in range(len(self.regions)):
- self.regions[iSec][iRegion].BlowingVel = vSolver.ExtractBlwVel(iSec, iRegion)
- self.regions[iSec][iRegion].DeltaStarExt = vSolver.ExtractDeltaStarExt(iSec, iRegion)
- self.regions[iSec][iRegion].xxExt = vSolver.ExtractxxExt(iSec, iRegion)
-
- blwVelAF = np.concatenate((self.regions[iSec][0].BlowingVel, self.regions[iSec][1].BlowingVel))
- deltaStarAF = np.concatenate((self.regions[iSec][0].BlowingVel, self.regions[iSec][1].DeltaStarExt[1:])) # Remove stagnation point doublon.
- xxAF = np.concatenate((self.regions[iSec][0].xxExt, self.regions[iSec][1].xxExt[1:])) # Remove stagnation point doublon.
- blwVelWK = self.regions[iSec][2].BlowingVel
- deltaStarWK = self.regions[2].DeltaStarExt
- xxWK = self.regions[2].xxExt
-
- self.parampp[iSec][0].u = blwVelAF[self.parampp[iSec][0].connectListElems.argsort()]
- self.parampp[iSec][1].u = blwVelWK[self.parampp[iSec][1].connectListElems.argsort()]
-
- self.parampp[iSec][0].deltaStar = deltaStarAF[self.parampp[iSec][0].connectListNodes.argsort()]
- self.parampp[iSec][1].deltaStar = deltaStarWK[self.parampp[iSec][1].connectListNodes.argsort()]
-
- self.parampp[iSec][0].xx = xxAF[self.parampp[iSec][0].connectListNodes.argsort()]
- self.parampp[iSec][1].xx = xxWK[self.parampp[iSec][1].connectListNodes.argsort()]
-
- for n in range(0, len(self.parampp[iSec])):
- for i in range(0, self.parampp[iSec][n].nE):
- self.parampp[iSec][n].boundary.setU(i, self.parampp[iSec][n].u[i])
-
- def RunInvSolver(self):
- self.dartSolver.run()
-
- def ResetSolver(self):
- self.dartSolver.reset()
\ No newline at end of file