diff --git a/dart/default.py b/dart/default.py
index e12c4ab8c0118f8dfcf4b2061cc41f840f94da4b..51b3757950f3a815b34009c09fd726d09e76179b 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -24,9 +24,7 @@ import dart
import tbox
def initBL(Re, Minf, CFL0):
- print('Re default.py', Re, Minf, CFL0)
solver = dart.ViscSolver(Re, Minf, CFL0)
- print('ok')
return solver
def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index f71ee9fc7605e1f06ff347585f606a543ef6945f..7de8cbdd2721cdca4f3389dc9774525c1564dfb8 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -10,8 +10,6 @@ using namespace dart;
BLRegion::BLRegion()
{
- std::cout << "Coucou de BLRegion\n" << std::endl;
-
closSolver = new Closures();
@@ -29,6 +27,7 @@ void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
LocMarkers=locM;
GlobMarkers=globM;
nMarkers = locM.size();
+ std::cout << "nMarkers" << nMarkers << std::endl;
/* Resize all variables */
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 7eb71544ff0ca171a7e0b743beef6ba58f37ae09..e3169f49babcdb2935bd6e2bcff0a1c79fa9f1f4 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -9,8 +9,6 @@ using namespace dart;
Closures::Closures()
{
- std::cout << "Coucou de Closures\n" << std::endl;
-
} // end Closures
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 7af50a9e0da44aee5a3368f05989dd107ca499d7..5378e868958f8fc1efe9ed70de5c9e0180c5490a 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -56,6 +56,9 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] == 0){
bl->U[iPoint * nVar + 4] = 0.03;
} // End if
+ else if (bl->Regime[iPoint] == 0){
+ bl->U[iPoint*nVar+4] = 0;
+ } // End else if
} // End InitialCondition
int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
@@ -69,6 +72,11 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
for (size_t i=0; i<nVar; ++i){
temp_U.push_back(bl->U[iPoint*nVar+i]);
} // End for
+ /*if (bl->name=="upper" && iPoint==67){
+ for (size_t k=0; k<nVar; ++k){
+ std::cout << "k" << temp_U[k];
+ }
+ }*/
// Initialize time integration variables
@@ -94,6 +102,11 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
while (innerIters < maxIt){
+
+ /*if (bl->name=="upper" && iPoint==67){
+ bl->PrintSolution(iPoint);
+ std::cout << "CFL " << CFL << std::endl;
+ }*/
if(innerIters % itFrozenJac == 0){
JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
}
@@ -103,6 +116,8 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
for(size_t k=0; k<nVar; ++k){
bl->U[iPoint*nVar+k] += slnIncr(k);
+ //bl->U[iPoint*nVar+0] = std::max(bl->U[iPoint*nVar+0], 1e-7);
+ //bl->U[iPoint*nVar+1] = std::max(bl->U[iPoint*nVar+1], 1.0005);
}
SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
@@ -110,7 +125,6 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
if (normSysRes <= tol * normSysRes0){
/* Successfull exit */
- //std::cout << "Point" << iPoint << " Convergence successfull. Niter: " << innerIters << " NormsysRes" << normSysRes << std::endl;
return 0;
}
@@ -124,7 +138,9 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
innerIters++;
} // End while (innerIters < maxIt)
- //std::cout << "Point" << iPoint << " Max iterations number reached. normSysRes/normSysRes0 " << normSysRes/normSysRes0 << std::endl;
+
+ ResetSolution(iPoint, bl, temp_U);
+ bl->PrintSolution(iPoint);
return innerIters;
} // End Integration
@@ -146,6 +162,54 @@ double TimeSolver::ComputeSoundSpeed(double gradU2)
return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
}
+void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U){
+
+ std::cout << "Reseting point" << iPoint << std::endl;
+ unsigned int nVar = bl->GetnVar();
+
+ /* Reset solution */
+
+ for(size_t k=0; k<nVar; ++k){
+ if (std::isnan(bl->U[iPoint*nVar+k])){
+ for (size_t j = 0; k<bl->GetnVar(); ++j){
+ if(!std::isnan(temp_U[j])){
+ bl->U[iPoint*bl->GetnVar()+j] = temp_U[j];
+ }
+ else{
+ bl->U[iPoint*nVar+j] = bl->U[(iPoint-1)*nVar+j];
+ }
+ break;
+ } // End if
+ } // End for
+
+ /* Reset closures */
+
+ 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->GetMe(iPoint), bl->name);
+ bl->Hstar[iPoint] = lamParam[0];
+ bl->Hstar2[iPoint] = lamParam[1];
+ bl->Hk[iPoint] = lamParam[2];
+ bl->Cf[iPoint] = lamParam[3];
+ bl->Cd[iPoint] = lamParam[4];
+ bl->Delta[iPoint] = lamParam[5];
+ bl->Cteq[iPoint] = lamParam[6];
+ bl->Us[iPoint] = lamParam[7];
+ } // 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->GetMe(iPoint), bl->name);
+ bl->Hstar[iPoint] = turbParam[0];
+ bl->Hstar2[iPoint] = turbParam[1];
+ bl->Hk[iPoint] = turbParam[2];
+ bl->Cf[iPoint] = turbParam[3];
+ bl->Cd[iPoint] = turbParam[4];
+ bl->Delta[iPoint] = turbParam[5];
+ bl->Cteq[iPoint] = turbParam[6];
+ bl->Us[iPoint] = turbParam[7];
+ } // End else if
+
+ } // End for
+}
+
double TimeSolver::GetCFL0(){return CFL0;}
void TimeSolver::SetCFL0(double _CFL0){CFL0 = _CFL0;}
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index e676590293fedd701dafa11dcc67553f6fbe12db..07c160dc85184c9da2f4aa6e56ca2ee4a58ef729 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -40,6 +40,8 @@ public:
private:
double SetTimeStep(double CFL, double dx, double invVel);
double ComputeSoundSpeed(double gradU2);
+ void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U);
+
};
} // namespace dart
#endif //WTimeSolver_H
\ No newline at end of file
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 633378a5b050f8b20c09b8be663fcd3bf596b3f6..dc2ec99107de7d8c4d8808f73947bffefba34ce5 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -20,15 +20,13 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
CFL0 = _CFL0;
std::cout << "--- Viscous solver setup ---\n"
- << "Reynolds number: " << Re
- << "Freestream Mach number: " << Minf
- << "Initial CFL number: " << CFL0
+ << "Reynolds number: " << Re << " \n"
+ << "Freestream Mach number: " << Minf << " \n"
+ << "Initial CFL number: " << CFL0 << " \n"
<< std::endl;
/* Initialzing boundary layer */
- //std::vector<dart::BLRegion> bl;
-
for (size_t i=0; i<3; ++i)
{
BLRegion *region = new BLRegion();
@@ -41,7 +39,10 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
bl[1]->name = "lower";
bl[2]->name = "wake";
- std::cout << " Boundary layer initialized.\n" << std::endl;
+ nbElmUpper = 0;
+ stagPtMvmt = false;
+
+ std::cout << "Initializing boundary layer" << std::endl;
} // end ViscSolver
@@ -78,21 +79,33 @@ void ViscSolver::SendExtVariables(size_t region, std::vector<double> _xxExt, std
int ViscSolver::Run(unsigned int couplIter){
- std::cout << "Entered run"
- << std::endl;
- /* Prepare time solver */
+ /* Prepare time integration */
- if (couplIter == 0){
- tSolver->SetCFL0(1);
+ //std::cout << "couplIter" << couplIter << "nbElmUpper: " << nbElmUpper << "bl[0]->GetnMarkers()" << bl[0]->GetnMarkers() << "stagPtMvmt" << stagPtMvmt << std::endl;
+
+ if ((couplIter == 0) || (nbElmUpper != bl[0]->GetnMarkers()) || (stagPtMvmt)){
+ tSolver->SetCFL0(0.5);
tSolver->SetitFrozenJac(1);
tSolver->SetinitSln(true);
+ if (couplIter != 0 && (nbElmUpper != bl[0]->GetnMarkers())){
+ stagPtMvmt = true;
+ std::cout << "Stagnation point moved" << std::endl;
+ }
+ else{
+ stagPtMvmt = false;
+ }
+ std::cout << "Restart solution" << std::endl;
}
+
else{
- tSolver->SetCFL0(1);
+ tSolver->SetCFL0(0.5);
tSolver->SetitFrozenJac(1);
tSolver->SetinitSln(false);
+ stagPtMvmt = false;
+ std::cout << "Updating solution" << std::endl;
}
+ nbElmUpper = bl[0]->GetnMarkers();
/* Set boundary condition */
@@ -116,10 +129,9 @@ int ViscSolver::Run(unsigned int couplIter){
bl[iRegion]->Regime[k]=1;
}
}
-
lockTrans = false;
- std::cout << "Starting " << bl[iRegion]->name << std::endl;
+ std::cout << "- " << bl[iRegion]->name << " region";
if (iRegion == 2){
/* Impose wake boundary condition */
@@ -128,6 +140,7 @@ int ViscSolver::Run(unsigned int couplIter){
UpperCond = std::vector<double>(bl[0]->U.end() - bl[0]->GetnVar(), bl[0]->U.end()); /* Base std::vector constructor to slice */
LowerCond = std::vector<double>(bl[1]->U.end() - bl[1]->GetnVar(), bl[1]->U.end()); /* Base std::vector constructor to slice */
bl[iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
+ std::cout << "\n" << std::endl;
lockTrans = true;
}
@@ -152,7 +165,7 @@ int ViscSolver::Run(unsigned int couplIter){
// Free transition.
if (bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 2] >= bl[iRegion]->GetnCrit()){
AverageTransition(iPoint, bl[iRegion], 0);
- std::cout << "Free transition x/c = " << bl[iRegion]->xtr << std::endl;
+ std::cout << ", free transition x/c = " << bl[iRegion]->xtr << ". Marker " << bl[iRegion]->transMarker << std::endl;
lockTrans = true;
} // End if
} // End if (!lockTrans).
@@ -274,6 +287,7 @@ void ViscSolver::ComputeDrag()
{
unsigned int nVar = bl[0]->GetnVar();
+ unsigned int lastWkPt = (bl[2]->GetnMarkers() - 1)*nVar;
/* Normalize friction and dissipation coefficients. */
/* std::vector<double> frictioncoeff_upper;
@@ -302,7 +316,7 @@ void ViscSolver::ComputeDrag()
/* Compute total drag coefficient. */
- Cdt = 2 * bl[2]->U[bl[2]->U.back()-5] * pow(bl[2]->U[bl[2]->U.back()-2],(bl[2]->U[bl[2]->U.back()-4]+5)/2);
+ Cdt = 2 * bl[2]->U[lastWkPt+0] * pow(bl[2]->U[lastWkPt+3],(bl[2]->U[lastWkPt+1]+5)/2);
/* Compute pressure drag coefficient. */
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index 0cf54920e1e85d6eff0b5d78dfad3b2cce9630cb..c1e6b07ffa77e2bc916dd3e728a61b6484dcec08 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -18,6 +18,8 @@ private:
double Re;
double Minf;
double CFL0;
+ size_t nbElmUpper;
+ bool stagPtMvmt;
TimeSolver *tSolver;
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 80f4942b5b2b87262a6c07876d3b9580e00b52a6..551f771a908724a19d8cfe5b123c92d8f0630fe7 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -60,7 +60,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 0.*math.pi/180
+ alpha = 5.*math.pi/180
M_inf = 0.
#user_xtr=[0.41,0.41]
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
index f65132bf5140085e838e3bd3c6c985e7cb43c8d5..18ad37fa44b4f2b213b467697b67e4454710d297 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -64,7 +64,7 @@ def main():
alpha = 2*math.pi/180
M_inf = 0.715
- user_xtr=[0.43, None]
+ user_xtr=[None, None]
user_Ti=None
mapMesh = 1
@@ -72,7 +72,7 @@ def main():
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 2
+ CFL0 = 1
# ========================================================================================
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index 6543ced2e408a8b12f6b8e6e38bb183e45ee1d2b..76f57dd1b14c2613015a6f733d0c983a45f618e2 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -55,8 +55,8 @@ def main():
# define flow variables
Re = 1e7
- alpha = 0.*math.pi/180
- M_inf = 0.
+ alpha = 2*math.pi/180
+ M_inf = 0.715
#user_xtr=[0.41,0.41]
#user_xtr=[0,None]
@@ -81,9 +81,9 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.005}
#pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
# set the problem and add medium, initial/boundary conditions
@@ -123,7 +123,7 @@ def main():
# visualize solution and plot results
- floD.initViewer(pbl)
+ #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)
#val=validation.Validation(case)
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 18a07c0a00f8ea27612f39bb488d4e39a501285a..d5e0fa8914cdf897b3b9fc452374f66260912046 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -50,12 +50,11 @@ class Coupler:
# Initilize meshes in c++ objects
couplIter = 0
+ cdPrev = 0
while couplIter < self.maxCouplIter:
- cdPrev = 0
# Run inviscid solver
-
self.iSolver.run()
# Extract inviscid boundary
@@ -74,8 +73,6 @@ class Coupler:
for n in range(0,len(self.group)):
self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList()
fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, 0, 1)
- #print('group[1].deltaStar', self.group[1].deltaStar)
- #print('group[1].xx', self.group[1].xx)
if couplIter == 0:
# Initialize mesh
self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
@@ -85,7 +82,6 @@ class Coupler:
self.vSolver.SendExtVariables(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
self.vSolver.SendExtVariables(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.SendExtVariables(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
- print('fro mpy : Mesh ok')
else:
self.vSolver.SendExtVariables(0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]], fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
self.vSolver.SendExtVariables(1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
@@ -94,13 +90,13 @@ class Coupler:
self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
- print('from py: InvBC ok')
+ print('Inviscid boundary imposed')
# Run viscous solver
exitCode = self.vSolver.Run(couplIter)
- print('From python exitcode: ', exitCode)
+ print('Viscous ops terminated with exit code ', exitCode)
- error = (self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt if self.vSolver.Cdt != 0 else 1
+ error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
if error <= self.couplTol:
print('')
@@ -109,6 +105,8 @@ class Coupler:
print('')
return 0
+ cdPrev = self.vSolver.Cdt
+
print('')
print(ccolors.ANSI_BLUE, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} > {:<2.3f}'
.format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index 780e4cc676527ab8d8b2918a523ae65ec706d5ea..4d83869c771a6eec9dbaf8dca1f3a9adbbe7eefa 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -54,7 +54,7 @@ class Coupler:
print('+---------------------------- VII Solver begin --------------------------------+')
print('')
- while it < 100:
+ while converged == 0:
# Run inviscid solver
print('+----------------------------- Inviscid solver --------------------------------+')
@@ -90,7 +90,7 @@ class Coupler:
# Check convergence and output coupling iteration.
- if 1 <= 0:
+ if error <= self.tol:
print('')
print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
.format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
index 624b9551b08ec6c7c78d5bc325c651d04bd39e82..4da9d5f729ea2db9c2d7d955f2dd8b1b00c917c1 100644
--- a/dart/viscous/Solvers/DIntegration.py
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -139,8 +139,8 @@ class Integration:
print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
quit()
- except Exception:
- print('it ', innerIters, ' Entered Exception')
+ except Exception as e:
+ print('it ', innerIters, ' Entered Exception', e)
nErrors += 1