diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index ffe14f8058309d6f5b2462101488514518308589..5f7c560002d6f9da2d0b186223c1edf1c035170e 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -31,4 +31,6 @@
#include "wTimeSolver.h"
#include "wViscFluxes.h"
#include "wClosures.h"
-#include "wSolutionInterp.h"
+#include "wInterpolator.h"
+#include "wViscMesh.h"
+#include "wInvBnd.h"
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index 29154d99b35476438d3686513978d8055ee9d180..531a1e5fc5820068259751f0e632049242bf2380 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -63,7 +63,9 @@ threads="1"
%include "wTimeSolver.h"
%include "wViscFluxes.h"
%include "wClosures.h"
-%include "wSolutionInterp.h"
+%include "wInterpolator.h"
+%include "wViscMesh.h"
+%include "wInvBnd.h"
%include "wFluid.h"
%include "wAssign.h"
diff --git a/dart/default.py b/dart/default.py
index 81e5153ec78fac11692a4e555ca09c384e32442f..3a03df5626170c58ef3be646579d518954072a69 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,10 +23,6 @@ from fwk.wutils import parseargs
import dart
import tbox
-def initBL(Re, Minf, CFL0, meshFactor):
- solver = dart.ViscSolver(Re, Minf, CFL0, meshFactor)
- return solver
-
def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
"""Initialize mesh and mesh writer
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 21a46657c2baf959cd560e78dded1dda357e26bb..91c6427dec6eb5197e9ad39e54d651666c3c9fc0 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -64,12 +64,16 @@ class Picard;
class Newton;
class FpeLSFunction;
class Adjoint;
+
+// Viscous
class ViscSolver;
class BLRegion;
class TimeSolver;
class ViscSolver;
class Closures;
-class SolutionInterp;
+class Interpolator;
+class ViscMesh;
+class InvBnd;
// General
class F0Ps;
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 3b394a79bea12ee2a20daa9e97b5e270543f670c..f370150c56cec068027c1e820b540ada956036d8 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -1,177 +1,165 @@
#include "wBLRegion.h"
#include "wClosures.h"
+#include "wViscMesh.h"
+#include "wInvBnd.h"
#include <iostream>
#include <vector>
#include <cmath>
-using namespace tbox;
using namespace dart;
-
BLRegion::BLRegion()
{
- closSolver = new Closures();
+ closSolver = new Closures();
+ mesh = new ViscMesh();
+ invBC = new InvBnd();
-
} // end BLRegion
-
BLRegion::~BLRegion()
-{
- delete closSolver;
- std::cout << "~BLRegion()\n";
+{
+ delete closSolver;
+ delete mesh;
+ std::cout << "~BLRegion()\n";
} // end ~BLRegion
void BLRegion::SetMesh(std::vector<double> locM, std::vector<double> globM)
{
- LocMarkers=locM;
- GlobMarkers=globM;
- nMarkers = locM.size();
-
- /* Resize all variables */
-
- U.resize(nVar * nMarkers);
- Ret.resize(nMarkers);
- Cd.resize(nMarkers);
- Cf.resize(nMarkers);
- Hstar.resize(nMarkers);
- Hstar2.resize(nMarkers);
- Hk.resize(nMarkers);
- Cteq.resize(nMarkers);
- Us.resize(nMarkers);
- DeltaStar.resize(nMarkers);
- Delta.resize(nMarkers);
- Me.resize(nMarkers);
- Rhoe.resize(nMarkers);
- VtExt.resize(nMarkers);
- DeltaStarExt.resize(nMarkers);
- XxExt.resize(nMarkers);
- ReExt.resize(nMarkers);
- Regime.resize(nMarkers);
- Blowing.resize(nMarkers-1);
+
+ 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
void BLRegion::SetExtVariables(std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
{
- XxExt = _xxExt;
- DeltaStarExt = _deltaStarExt;
+ if (_xxExt.size() != mesh->GetnMarkers() || _deltaStarExt.size() != mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::BLRegion: Miss match between mesh and external variables.");
+ }
+
+ mesh->SetExt(_xxExt);
+ invBC->SetDeltaStar(_deltaStarExt);
}
void BLRegion::SetInvBC(std::vector<double> _Ue,
std::vector<double> _Me,
std::vector<double> _Rhoe)
{
- VtExt = _Ue;
- Me = _Me;
- Rhoe = _Rhoe;
-}// End SetInvBC
+ invBC->SetMe(_Me);
+ invBC->SetVt(_Ue);
+ invBC->SetRhoe(_Rhoe);
+} // End SetInvBC
void BLRegion::SetStagBC(double Re)
{
- double dv0 = (VtExt[1] - VtExt[0])/(LocMarkers[1] - LocMarkers[0]);
-
- if (dv0 > 0){
- U[0] = sqrt(0.075/(Re*dv0)); // Theta
- } // End if
- else{
- U[0] = 1e-6; // Theta
- } // end else
- U[1] = 2.23; // H
- U[2] = 0; // N
- U[3] = VtExt[0]; // ue
- U[4] = 0; // Ctau
-
- Ret[0] = U[3]*U[0]*Re;
-
- if (Ret[0] < 1)
- {
- Ret[0] = 1;
- U[3] = Ret[0]/(U[0]*Re);
- } // End if
-
- //std::cout << "name " << name << " dv0 " << dv0 << " U3 " << U[3] << " U[0] " << U[0] << " \n" << std::scientific;
-
-
- DeltaStar[0] = U[0] * U[1];
-
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
-
- Hstar[0] = ClosParam[0];
- Hstar2[0] = ClosParam[1];
- Hk[0] = ClosParam[2];
- Cf[0] = ClosParam[3];
- Cd[0] = ClosParam[4];
- Delta[0] = ClosParam[5];
- Cteq[0] = ClosParam[6];
- Us[0] = ClosParam[7];
+ double dv0 = (invBC->GetVt(1) - invBC->GetVt(0)) / (mesh->GetLoc(1) - mesh->GetLoc(0));
+
+ if (dv0 > 0)
+ {
+ U[0] = sqrt(0.075 / (Re * dv0)); // Theta
+ } // End if
+ else
+ {
+ U[0] = 1e-6; // Theta
+ } // end else
+ U[1] = 2.23; // H
+ U[2] = 0; // N
+ U[3] = invBC->GetVt(0); // ue
+ U[4] = 0; // Ctau
+
+ Ret[0] = U[3] * U[0] * Re;
+
+ if (Ret[0] < 1)
+ {
+ Ret[0] = 1;
+ U[3] = Ret[0] / (U[0] * Re);
+ } // End if
+
+ //std::cout << "name " << name << " dv0 " << dv0 << " U3 " << U[3] << " U[0] " << U[0] << " \n" << std::scientific;
+
+ DeltaStar[0] = U[0] * U[1];
+
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBC->GetMe(0), name);
+
+ Hstar[0] = ClosParam[0];
+ Hstar2[0] = ClosParam[1];
+ Hk[0] = ClosParam[2];
+ Cf[0] = ClosParam[3];
+ Cd[0] = ClosParam[4];
+ Delta[0] = ClosParam[5];
+ Cteq[0] = ClosParam[6];
+ Us[0] = ClosParam[7];
} // End SetStagBC
void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
{
- 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] = VtExt[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. */
-
- if (Ret[0] < 1){
- Ret[0] = 1;
- U[3] = Ret[0] / (U[0]*Re);
- }
- DeltaStar[0] = U[0]*U[1]; /* deltaStar */
-
- /* Laminar closures. */
-
- std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
-
- Hstar[0] = ClosParam[0];
- Hstar2[0] = ClosParam[1];
- Hk[0] = ClosParam[2];
- Cf[0] = ClosParam[3];
- Cd[0] = ClosParam[4];
- Delta[0] = ClosParam[5];
- Cteq[0] = ClosParam[6];
- Us[0] = ClosParam[7];
-
- for (size_t k=0; k<nMarkers; ++k)
- {
- Regime[k] = 1;
- }
+ 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[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. */
+
+ if (Ret[0] < 1)
+ {
+ Ret[0] = 1;
+ U[3] = Ret[0] / (U[0] * Re);
+ }
+ DeltaStar[0] = U[0] * U[1]; /* deltaStar */
+
+ /* Laminar closures. */
+
+ std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], invBC->GetMe(0), name);
+
+ Hstar[0] = ClosParam[0];
+ Hstar2[0] = ClosParam[1];
+ Hk[0] = ClosParam[2];
+ Cf[0] = ClosParam[3];
+ Cd[0] = ClosParam[4];
+ Delta[0] = ClosParam[5];
+ Cteq[0] = ClosParam[6];
+ Us[0] = ClosParam[7];
+
+ for (size_t k = 0; k < mesh->GetnMarkers(); ++k)
+ {
+ Regime[k] = 1;
+ }
}
-void BLRegion::ComputeBlwVel()
-{
- for (size_t iPoint=1; iPoint<nMarkers; ++iPoint){
-
- DeltaStar[iPoint] = U[iPoint*nVar+0] * U[iPoint*nVar+1];
- Blowing[iPoint-1] = (Rhoe[iPoint]*VtExt[iPoint]*DeltaStar[iPoint]
- -Rhoe[iPoint-1]*VtExt[iPoint-1]*DeltaStar[iPoint-1])/(Rhoe[iPoint]*(LocMarkers[iPoint]-LocMarkers[iPoint-1]));
-
- }
-}
-
-unsigned int BLRegion::GetnVar(){return nVar;}
-size_t BLRegion::GetnMarkers(){return nMarkers;}
-double BLRegion::GetnCrit(){return nCrit;}
-double BLRegion::GetLocMarkers(size_t iPoint){return LocMarkers[iPoint];}
-double BLRegion::GetGlobMarkers(size_t iPoint){return GlobMarkers[iPoint];}
-double BLRegion::GetVtExt(size_t iPoint){return VtExt[iPoint];}
-double BLRegion::GetMe(size_t iPoint){return Me[iPoint];}
-double BLRegion::GetRhoe(size_t iPoint){return Rhoe[iPoint];}
-double BLRegion::GetXxExt(size_t iPoint){return XxExt[iPoint];}
-double BLRegion::GetDeltaStarExt(size_t iPoint){return DeltaStarExt[iPoint];}
-
void BLRegion::PrintSolution(size_t iPoint)
{
- std::cout << "\nPoint " << iPoint << " \n" << std::scientific;
- std::cout << "Theta " << U[iPoint * nVar + 0] << " \n" << std::scientific;
- std::cout << "H " << U[iPoint * nVar + 1] << " \n" << std::scientific;
- std::cout << "N " << U[iPoint * nVar + 2] << " \n" << std::scientific;
- std::cout << "ue " << U[iPoint * nVar + 3] << " \n" << std::scientific;
- std::cout << "Ct " << U[iPoint * nVar + 4] << " \n" << std::scientific;
- std::cout << "Regime " << Regime[iPoint] << " \n" << std::scientific;
- }
\ No newline at end of file
+ std::cout << "\nPoint " << iPoint << " \n"
+ << std::scientific;
+ std::cout << "Theta " << U[iPoint * nVar + 0] << " \n"
+ << std::scientific;
+ std::cout << "H " << U[iPoint * nVar + 1] << " \n"
+ << std::scientific;
+ std::cout << "N " << U[iPoint * nVar + 2] << " \n"
+ << std::scientific;
+ std::cout << "ue " << U[iPoint * nVar + 3] << " \n"
+ << std::scientific;
+ std::cout << "Ct " << U[iPoint * nVar + 4] << " \n"
+ << std::scientific;
+ std::cout << "Regime " << Regime[iPoint] << " \n"
+ << std::scientific;
+}
\ No newline at end of file
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index 244fde997b9d0d02d4cd7f459b4a9430c7c22307..ec24d0e2de44694750f13993610899118630fe2d 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -3,6 +3,7 @@
#include "dart.h"
#include "wClosures.h"
+#include "wViscMesh.h"
#include <vector>
#include <string>
@@ -13,88 +14,67 @@ class DART_API BLRegion
{
private:
+ unsigned int nVar = 5;
- unsigned int nVar = 5;
+ /* Transition related variables */
- /* Mesh parameters */
+ double nCrit = 9.0;
- size_t nMarkers;
- std::vector<double> GlobMarkers; /* Position of the mesh elements in the global frame of reference */
- std::vector<double> LocMarkers; /* Position of the mesh elements in the local frame of reference */
-
- /* External flow variables */
-
- std::vector<double> Me;
- std::vector<double> Rhoe;
- std::vector<double> VtExt;
- std::vector<double> DeltaStarExt;
- std::vector<double> XxExt;
- std::vector<double> ReExt;
-
- /* Transition related variables */
-
- double nCrit=9.0;
public:
- Closures *closSolver;
-
- std::string name;
-
- size_t coarseSize; // Size of the coarse mesh on the corresponding region (only used if mesh mapping is enabled).
- size_t fineSize; // Size of the coarse mesh on the corresponding region (only used if mesh mapping is enabled).
-
- /* Boundary layer variables */
-
- std::vector<double> U; /* Solution vector (theta, H, N, ue, Ct) */
- std::vector<double> Ret; /* Reynolds number based on the momentum thickness (theta) */
- std::vector<double> Cd; /* Local dissipation coefficient */
- std::vector<double> Cf; /* Local friction coefficient */
- std::vector<double> Hstar;
- std::vector<double> Hstar2;
- std::vector<double> Hk;
- std::vector<double> Cteq;
- std::vector<double> Us;
- std::vector<double> Delta;
-
- std::vector<double> DeltaStar;
-
- /* Transition related variables */
+ ViscMesh *mesh;
+ InvBnd *invBC;
+ Closures *closSolver;
+
+ std::string name; /* Name of the region */
+
+ /* Boundary layer variables */
+
+ std::vector<double> U; /* Solution vector (theta, H, N, ue, Ct) */
+ std::vector<double> Ret; /* Reynolds number based on the momentum thickness (theta) */
+ std::vector<double> Cd; /* Local dissipation coefficient */
+ std::vector<double> Cf; /* Local friction coefficient */
+ std::vector<double> Hstar;
+ std::vector<double> Hstar2;
+ std::vector<double> Hk;
+ std::vector<double> Cteq;
+ std::vector<double> Us;
+ std::vector<double> Delta;
+
+ std::vector<double> DeltaStar;
+
+ /* Transition related variables */
+
+ std::vector<int> Regime;
+ std::vector<double> Blowing;
+
+ unsigned int transMarker;
+ double xtrF;
+ double xtr;
+
+ BLRegion();
+ ~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);
+
+ /* Boundary conditions */
+ void SetStagBC(double Re);
+ void SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond);
+
+ /* General access */
+ unsigned int GetnVar() {return nVar;};
+ double GetnCrit() {return nCrit;};
+
+ /* Others */
+ void PrintSolution(size_t iPoint);
+ void ComputeBlwVel();
- std::vector<int> Regime;
- std::vector<double> Blowing;
-
- unsigned int transMarker;
- double xtrF;
- double xtr;
-
- BLRegion();
- ~BLRegion();
- 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 SetStagBC(double Re);
- void SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond);
-
- unsigned int GetnVar();
- size_t GetnMarkers();
- double GetnCrit();
- double GetLocMarkers(size_t iPoint);
- double GetGlobMarkers(size_t iPoint);
- double GetVtExt(size_t iPoint);
- double GetMe(size_t iPoint);
- double GetRhoe(size_t iPoint);
- double GetXxExt(size_t iPoint);
- double GetDeltaStarExt(size_t iPoint);
-
-
- void PrintSolution(size_t iPoint);
- void ComputeBlwVel();
-
-private:
};
} // namespace dart
#endif //WBLREGION_H
\ No newline at end of file
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index bc83e72f5196fd2dc441306b9ed1f17c52f17b68..edaf5109f34b53add692ab6cdd99bd8c563ce5ef 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -7,11 +7,7 @@ using namespace tbox;
using namespace dart;
// This is an empty constructor
-Closures::Closures()
-{
-
-} // end Closures
-
+Closures::Closures() {}
Closures::~Closures()
{
@@ -21,139 +17,151 @@ Closures::~Closures()
std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret, double Me, std::string name)
{
- double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me); // Kinematic shape factor
- double Hstar2 = (0.064/(Hk-0.8)+0.251)*Me*Me; // Density shape parameter
- double delta = (3.15+ H + (1.72/(Hk-1)))*theta;
- double Hks = Hk - 4.35;
- double Hstar;
- if (Hk <= 4.35){
- double dens = Hk + 1;
- Hstar = 0.0111*Hks*Hks/dens - 0.0278*Hks*Hks*Hks/dens + 1.528 -0.0002*(Hks*Hk)*(Hks*Hk);
- Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
- }
- else if(Hk > 4.35)
- {
- Hstar = 1.528 + 0.015*Hks*Hks/Hk;
- Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); // Whitfield's minor additional compressibility correction
- }
-
- double tmp;
- double Cf;
- double Cd;
- if (Hk < 5.5){
- tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1);
- Cf = (-0.07 + 0.0727*tmp)/Ret;}
- else if (Hk >= 5.5){
- tmp = 1 - 1/(Hk-4.5);
- Cf = (-0.07 + 0.015*tmp*tmp) /Ret;}
- if (Hk < 4){
- Cd = ( 0.00205*pow(4.0-Hk, 5.5) + 0.207 ) * (Hstar/(2*Ret));}
- else if (Hk >= 4){
- double HkCd = (Hk-4)*(Hk-4);
- double denCd = 1+0.02*HkCd;
- Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret));}
- if (name == "wake"){
- Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
- Cf = 0;}
-
- double Cteq = 0;
- double Us = 0;
-
- std::vector<double> ClosureVars;
-
- ClosureVars.push_back(Hstar);
- ClosureVars.push_back(Hstar2);
- ClosureVars.push_back(Hk);
- ClosureVars.push_back(Cf);
- ClosureVars.push_back(Cd);
- ClosureVars.push_back(delta);
- ClosureVars.push_back(Cteq);
- ClosureVars.push_back(Us);
- return ClosureVars;
+ double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me); // Kinematic shape factor
+ double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me; // Density shape parameter
+ double delta = (3.15 + H + (1.72 / (Hk - 1))) * theta;
+ double Hks = Hk - 4.35;
+ double Hstar;
+ if (Hk <= 4.35)
+ {
+ double dens = Hk + 1;
+ Hstar = 0.0111 * Hks * Hks / dens - 0.0278 * Hks * Hks * Hks / dens + 1.528 - 0.0002 * (Hks * Hk) * (Hks * Hk);
+ Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); // Whitfield's minor additional compressibility correction
+ }
+ else if (Hk > 4.35)
+ {
+ Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
+ Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); // Whitfield's minor additional compressibility correction
+ }
+
+ double tmp;
+ double Cf;
+ double Cd;
+ if (Hk < 5.5)
+ {
+ tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1);
+ Cf = (-0.07 + 0.0727 * tmp) / Ret;
+ }
+ else if (Hk >= 5.5)
+ {
+ tmp = 1 - 1 / (Hk - 4.5);
+ Cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
+ }
+ if (Hk < 4)
+ {
+ Cd = (0.00205 * pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
+ }
+ else if (Hk >= 4)
+ {
+ double HkCd = (Hk - 4) * (Hk - 4);
+ double denCd = 1 + 0.02 * HkCd;
+ Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
+ }
+ if (name == "wake")
+ {
+ Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2 * Ret));
+ Cf = 0;
+ }
+
+ double Cteq = 0;
+ double Us = 0;
+
+ std::vector<double> ClosureVars;
+
+ ClosureVars.push_back(Hstar);
+ ClosureVars.push_back(Hstar2);
+ ClosureVars.push_back(Hk);
+ ClosureVars.push_back(Cf);
+ ClosureVars.push_back(Cd);
+ ClosureVars.push_back(delta);
+ ClosureVars.push_back(Cteq);
+ ClosureVars.push_back(Us);
+ return ClosureVars;
}
std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name)
{
- double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
- Hk = std::max(Hk, 1.0005);
- double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
-
- double Fc = sqrt(1+0.2*Me*Me);
-
- double H0;
- if(Ret <= 400)
- {
- H0 = 4;
- }
- else
- {
- H0 = 3 + (400/Ret);
- }
-
- double Hstar;
- if (Hk <= H0)
- {
- Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))*((H0-Hk)/(H0-1)))*(1.5/(Hk+0.5))+1.5+4/Ret;
- }
- else
- {
- Hstar = (Hk-H0)*(Hk-H0)*(0.007*log(Ret)/((Hk-H0+4/log(Ret))*(Hk-H0+4/log(Ret))) + 0.015/Hk) + 1.5 + 4/Ret;
- }
- Hstar = (Hstar + 0.028*Me*Me)/(1+0.014*Me*Me); /* Whitfield's minor additional compressibility correction */
-
- double Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)); /* Equivalent normalized wall slip velocity */
- double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta);
- double Ctcon = 0.5/(6.7*6.7*0.75);
-
- double Cf;
- double Cd;
- double Hkc;
- double Cdw;
- double Cdd;
- double Cdl;
- double Cteq;
- int n;
- if (name == "wake")
- {
- Us = std::min(Us, 0.99995);
- n = 2; // two wake halves
- Cf = 0; // no friction inside the wake
- Hkc = Hk - 1;
- Cdw = 0; // wall contribution
- Cdd = (0.995-Us)*Ct*Ct; // turbulent outer layer contribution
- Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret; // laminar stress contribution to outer layer
- Cteq = sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
- }
- else
- {
- if (Us > 0.95){
- Us = 0.98;
+ double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
+ Hk = std::max(Hk, 1.0005);
+ double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
+
+ double Fc = sqrt(1 + 0.2 * Me * Me);
+
+ double H0;
+ if (Ret <= 400)
+ {
+ H0 = 4;
+ }
+ else
+ {
+ H0 = 3 + (400 / Ret);
+ }
+
+ double Hstar;
+ if (Hk <= H0)
+ {
+ Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
+ }
+ else
+ {
+ Hstar = (Hk - H0) * (Hk - H0) * (0.007 * log(Ret) / ((Hk - H0 + 4 / log(Ret)) * (Hk - H0 + 4 / log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
+ }
+ Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me); /* Whitfield's minor additional compressibility correction */
+
+ double Us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H)); /* Equivalent normalized wall slip velocity */
+ double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
+ double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
+
+ double Cf;
+ double Cd;
+ double Hkc;
+ double Cdw;
+ double Cdd;
+ double Cdl;
+ double Cteq;
+ int n;
+ if (name == "wake")
+ {
+ Us = std::min(Us, 0.99995);
+ n = 2; // two wake halves
+ Cf = 0; // no friction inside the wake
+ Hkc = Hk - 1;
+ Cdw = 0; // wall contribution
+ Cdd = (0.995 - Us) * Ct * Ct; // turbulent outer layer contribution
+ Cdl = 0.15 * (0.995 - Us) * (0.995 - Us) / Ret; // laminar stress contribution to outer layer
+ Cteq = sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - Us) * (Hk * Hk) * H)); // Here it is Cteq^0.5
+ }
+ else
+ {
+ if (Us > 0.95)
+ {
+ Us = 0.98;
+ }
+ n = 1;
+ Cf = 1 / Fc * (0.3 * exp(-1.33 * Hk) * pow((log10(Ret / Fc)), -1.74 - 0.31 * Hk) + 0.00011 * (tanh(4 - Hk / 0.875) - 1));
+ Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
+
+ // Dissipation coefficient
+ double Hmin = 1 + 2.1 / log(Ret);
+ double Fl = (Hk - 1) / (Hmin - 1);
+ double Dfac = 0.5 + 0.5 * tanh(Fl);
+ Cdw = 0.5 * (Cf * Us) * Dfac;
+ Cdd = (0.995 - Us) * Ct * Ct;
+ Cdl = 0.15 * (0.995 - Us) * (0.995 - Us) / Ret;
+ Cteq = sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - Us) * (Hk * Hk) * H)); // Here it is Cteq^0.5
}
- n = 1;
- Cf = 1/Fc*(0.3*exp(-1.33*Hk)*pow((log10(Ret/Fc)),-1.74-0.31*Hk) + 0.00011*(tanh(4-Hk/0.875) - 1));
- Hkc = std::max(Hk-1-18/Ret, 0.01);
-
- // Dissipation coefficient
- double Hmin = 1 + 2.1/log(Ret);
- double Fl = (Hk-1)/(Hmin-1);
- double Dfac = 0.5+0.5*tanh(Fl);
- Cdw = 0.5*(Cf*Us) * Dfac;
- Cdd = (0.995-Us)*Ct*Ct;
- Cdl = 0.15*(0.995-Us)*(0.995-Us)/Ret;
- Cteq = sqrt(Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)); // Here it is Cteq^0.5
- }
- Cd = n*(Cdw+ Cdd + Cdl);
-
- std::vector<double> ClosureVars;
-
- ClosureVars.push_back(Hstar);
- ClosureVars.push_back(Hstar2);
- ClosureVars.push_back(Hk);
- ClosureVars.push_back(Cf);
- ClosureVars.push_back(Cd);
- ClosureVars.push_back(delta);
- ClosureVars.push_back(Cteq);
- ClosureVars.push_back(Us);
- return ClosureVars;
+ Cd = n * (Cdw + Cdd + Cdl);
+
+ std::vector<double> ClosureVars;
+
+ ClosureVars.push_back(Hstar);
+ ClosureVars.push_back(Hstar2);
+ ClosureVars.push_back(Hk);
+ ClosureVars.push_back(Cf);
+ ClosureVars.push_back(Cd);
+ ClosureVars.push_back(delta);
+ ClosureVars.push_back(Cteq);
+ ClosureVars.push_back(Us);
+ return ClosureVars;
}
diff --git a/dart/src/wClosures.h b/dart/src/wClosures.h
index a51871730f564bc689473f3729f776c0730d1250..2946cf76748751ff5b246f89a425c36eca6c2627 100644
--- a/dart/src/wClosures.h
+++ b/dart/src/wClosures.h
@@ -11,10 +11,10 @@ namespace dart
class DART_API Closures
{
public:
- Closures();
- ~Closures();
- std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name);
- std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name);
+ Closures();
+ ~Closures();
+ std::vector<double> LaminarClosures(double theta, double H, double Ret, double Me, std::string name);
+ std::vector<double> TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name);
};
} // namespace dart
#endif //WClosures_H
\ No newline at end of file
diff --git a/dart/src/wInterpolator.cpp b/dart/src/wInterpolator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..afef795e64dbf40310acbb880b08c92f1748a29c
--- /dev/null
+++ b/dart/src/wInterpolator.cpp
@@ -0,0 +1,18 @@
+#include "wInterpolator.h"
+#include "wClosures.h"
+#include "wViscMesh.h"
+#include "wInvBnd.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+using namespace dart;
+
+Interpolator::Interpolator()
+{
+} // end Interpolator
+
+Interpolator::~Interpolator()
+{
+} // end ~Interpolator
+
diff --git a/dart/src/wInterpolator.h b/dart/src/wInterpolator.h
new file mode 100644
index 0000000000000000000000000000000000000000..45a74bc2602be4aa89d680721e5eb25d62bd2e9c
--- /dev/null
+++ b/dart/src/wInterpolator.h
@@ -0,0 +1,22 @@
+#ifndef WINTERPOLATOR_H
+#define WINTERPOLATOR_H
+
+#include "dart.h"
+#include "wClosures.h"
+#include "wViscMesh.h"
+#include <vector>
+#include <string>
+
+namespace dart
+{
+
+class DART_API Interpolator
+{
+ ViscMesh *vmesh;
+
+public:
+
+ void GenerateMesh();
+};
+} // namespace dart
+#endif //WINTERPOLATOR_H
\ No newline at end of file
diff --git a/dart/src/wInvBnd.cpp b/dart/src/wInvBnd.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4fd1aceb80353075e861c3f13eefbfd46564a9e0
--- /dev/null
+++ b/dart/src/wInvBnd.cpp
@@ -0,0 +1,34 @@
+/* ----------------------------------------------------------------------------------- */
+/* */
+/* Coupled integral boundary layer solver */
+/* */
+/* Université de Liège */
+/* February 2022 */
+/* Author: Paul Dechamps */
+/* v0.1 */
+/* ----------------------------------------------------------------------------------- */
+
+/* --- Project includes -------------------------------------------------------------- */
+
+#include "wInvBnd.h"
+#include <iostream>
+#include <vector>
+
+/* --- Namespaces -------------------------------------------------------------------- */
+
+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
new file mode 100644
index 0000000000000000000000000000000000000000..9a546fa9e61d89afa2d5eeab7ce0fdbb6b7e1856
--- /dev/null
+++ b/dart/src/wInvBnd.h
@@ -0,0 +1,40 @@
+#ifndef WINVBND_H
+#define WINVBND_H
+
+#include "dart.h"
+#include "wBLRegion.h"
+#include <vector>
+
+namespace dart
+{
+
+class DART_API InvBnd
+{
+ private:
+ std::vector<double> Me,
+ Vt,
+ Rhoe,
+ DeltaStar;
+
+ public:
+
+ std::vector<double> BlowingVel;
+
+ InvBnd();
+ ~InvBnd();
+
+ double GetMe(size_t iPoint) const {return Me[iPoint];};
+ double GetVt(size_t iPoint) const {return Vt[iPoint];};
+ double GetRhoe(size_t iPoint) const {return Rhoe[iPoint];};
+ double GetDeltaStar(size_t iPoint) const {return DeltaStar[iPoint];};
+
+ 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 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/wNewton.cpp b/dart/src/wNewton.cpp
index 28efba8418e0a040d7100405ff4a5f9f630188e8..993095bb097df8f7d136967d86bb3a910f431902 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -89,11 +89,11 @@ bool Newton::run()
// Display current freestream conditions
if (printOn)
{
- std::cout << std::setprecision(2)
- << "- Mach " << pbl->M_inf << ", "
- << pbl->alpha * 180 / 3.14159 << " deg AoA, "
- << pbl->beta * 180 / 3.14159 << " deg AoS"
- << std::endl;
+ std::cout << std::setprecision(2)
+ << "- Mach " << pbl->M_inf << ", "
+ << pbl->alpha * 180 / 3.14159 << " deg AoA, "
+ << pbl->beta * 180 / 3.14159 << " deg AoS"
+ << std::endl;
}
// Initialize solver loop
@@ -122,21 +122,21 @@ bool Newton::run()
// Display residual
if (printOn)
{
- std::cout << std::setw(8) << "N_Iter"
- << std::setw(8) << "L_Iter"
- << std::setw(8) << "f_eval"
- << std::setw(12) << "Cl"
- << std::setw(12) << "Cd"
- << std::setw(15) << "Rel_Res[phi]"
- << std::setw(15) << "Abs_Res[phi]" << std::endl;
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(8) << nIt
- << std::setw(8) << 0
- << std::setw(8) << 1
- << std::setw(12) << "-"
- << std::setw(12) << "-"
- << std::setw(15) << log10(relRes)
- << std::setw(15) << log10(absRes) << "\n";
+ std::cout << std::setw(8) << "N_Iter"
+ << std::setw(8) << "L_Iter"
+ << std::setw(8) << "f_eval"
+ << std::setw(12) << "Cl"
+ << std::setw(12) << "Cd"
+ << std::setw(15) << "Rel_Res[phi]"
+ << std::setw(15) << "Abs_Res[phi]" << std::endl;
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(8) << nIt
+ << std::setw(8) << 0
+ << std::setw(8) << 1
+ << std::setw(12) << "-"
+ << std::setw(12) << "-"
+ << std::setw(15) << log10(relRes)
+ << std::setw(15) << log10(absRes) << "\n";
}
do
@@ -243,13 +243,14 @@ bool Newton::run()
// Check the solution
if (relRes < relTol || absRes < absTol)
{
- if (printOn){
- std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
- if (verbose > 1)
- std::cout << "Newton solver CPU" << std::endl
- << tms;
- std::cout << std::endl;
- }
+ if (printOn)
+ {
+ std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
+ if (verbose > 1)
+ std::cout << "Newton solver CPU" << std::endl
+ << tms;
+ std::cout << std::endl;
+ }
return true;
}
else if (std::isnan(relRes))
diff --git a/dart/src/wNewton.h b/dart/src/wNewton.h
index f7b632f6090c1aef6216571229f76c5c7e4286b6..3b6f403b067810009b57ede374a1eb37aa01253b 100644
--- a/dart/src/wNewton.h
+++ b/dart/src/wNewton.h
@@ -42,7 +42,7 @@ public:
int maxLsIt; ///< max number of line search
double avThrsh; ///< residual threshold to deacrease artificial viscosity
- bool printOn=true;
+ bool printOn = true;
private:
double mCOv; ///< variable cut-off Mach number
diff --git a/dart/src/wSolutionInterp.cpp b/dart/src/wSolutionInterp.cpp
deleted file mode 100644
index ba32bfb24687c21f45bb3283ad0b88cbe234cbf4..0000000000000000000000000000000000000000
--- a/dart/src/wSolutionInterp.cpp
+++ /dev/null
@@ -1,121 +0,0 @@
-#define _USE_MATH_DEFINES
-
-#include "wSolutionInterp.h"
-#include <iostream>
-#include <unordered_set>
-#include <fstream>
-#include <iomanip>
-
-using namespace dart;
-
-SolutionInterp::SolutionInterp(unsigned int _meshFactor)
-{
- meshFactor = _meshFactor;
- CoarseMeshInit = false;
-
-} // end SolutionInterp
-
-
-SolutionInterp::~SolutionInterp()
-{
- std::cout << "~SolutionInterp()\n";
-} // end ~SolutionInterp
-
-void SolutionInterp::InputCoarseMesh(std::vector<double> _CoarseLocMesh, std::vector<double> _CoarseGlobMesh)
-{
- CoarseLocMesh = _CoarseLocMesh;
- CoarseGlobMesh = _CoarseGlobMesh;
- CoarseMeshInit = true;
-}
-void SolutionInterp::InputFineMesh(std::vector<double> _fineLocMesh, std::vector<double> _fineGlobMesh)
-{
- FineLocMesh=_fineLocMesh;
- FineGlobMesh=_fineGlobMesh;
-}
-
-std::vector<double> SolutionInterp::InterpolateToFineMesh(std::vector<double> cVector, size_t fineSize)
-{
- if (meshFactor<=1)
- {
- std::cout << "Can not create finer mesh while meshFactor is " << meshFactor << std::endl;
- return cVector;
- }
- else
- {
- std::vector<double> fVector(fineSize-1);
- double dv;
- for (size_t cPoint=0; cPoint < cVector.size()-1; ++cPoint)
- {
- dv = cVector[cPoint + 1] - cVector[cPoint];
- for (size_t fPoint=0; fPoint<meshFactor; ++fPoint)
- {
- fVector[meshFactor * cPoint + fPoint] = cVector[cPoint] + fPoint * dv / meshFactor;
- }
- }
- fVector.push_back(cVector.back());
- return fVector;
- }
-}
-
-std::vector<double> SolutionInterp::BlowingFineToCoarse(std::vector<double> fVector, size_t coarseSize)
-{
- std::vector<double> cVector(coarseSize);
-
- for (size_t cCell=0; cCell<cVector.size(); ++cCell)
- {
- for (size_t k=0; k<meshFactor; ++k)
- {
- cVector[cCell] += fVector[cCell*meshFactor+k];
- }
- cVector[cCell]/=meshFactor;
- }
- cVector.back() = fVector.back();
- return cVector;
-}
-
-std::vector<double> SolutionInterp::VarsFineToCoarse(std::vector<double> fVector, size_t coarseSize)
-{
- std::vector<double> cVector(coarseSize);
-
- for (size_t cPoint=0; cPoint<cVector.size(); ++cPoint)
- {
- cVector[cPoint] = fVector[cPoint * meshFactor];
- }
- cVector.back() = fVector.back();
- return cVector;
-}
-
-
-std::vector<double> SolutionInterp::CreateFineMesh(std::vector<double> cMesh)
-{
- if (meshFactor<1)
- {
- throw std::runtime_error("SolutionInterp::CreateFineMesh: Mesh factor should be greater than 1.");
- }
-
- else
- {
-
- double dx;
-
- /* Local markers */
-
- std::vector<double> fMesh;
- for (size_t iPoint=0; iPoint<cMesh.size()-1; ++iPoint)
- {
- dx = cMesh[iPoint + 1] - cMesh[iPoint];
-
- for(size_t iRefinedPoint = 0; iRefinedPoint < meshFactor; ++iRefinedPoint)
- {
- fMesh.push_back(cMesh[iPoint] + iRefinedPoint * dx/meshFactor);
- }
- }
- fMesh.push_back(cMesh.back());
- return fMesh;
- }
-}
-
-void SolutionInterp::ResetmeshFactor(unsigned int newmeshFactor)
-{
- meshFactor = newmeshFactor;
-}
\ No newline at end of file
diff --git a/dart/src/wSolutionInterp.h b/dart/src/wSolutionInterp.h
deleted file mode 100644
index a16b2444ce057bdc54e09fa054e91c3c30865819..0000000000000000000000000000000000000000
--- a/dart/src/wSolutionInterp.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef WSolutionInterp_H
-#define WSolutionInterp_H
-
-#include "dart.h"
-#include "wBLRegion.h"
-#include <vector>
-#include <string>
-
-namespace dart
-{
-
-class DART_API SolutionInterp
-{
-private:
- unsigned int meshFactor;
-public:
- std::vector<double> FineLocMesh;
- std::vector<double> FineGlobMesh;
- std::vector<double> CoarseLocMesh;
- std::vector<double> CoarseGlobMesh;
- bool CoarseMeshInit;
-public:
- SolutionInterp(unsigned int _meshFactor=1);
- ~SolutionInterp();
-
- void InputFineMesh(std::vector<double> _fineLocMesh, std::vector<double> _fineGlobMesh);
- std::vector<double> BlowingFineToCoarse(std::vector<double> fVector, size_t coarseSize);
- std::vector<double> VarsFineToCoarse(std::vector<double> fVector, size_t coarseSize);
- std::vector<double> CreateFineMesh(std::vector<double> cMesh);
- void ResetmeshFactor(unsigned int newmeshFactor);
- std::vector<double> InterpolateToFineMesh(std::vector<double> cVector, size_t fineSize);
- void InputCoarseMesh(std::vector<double> _CoarseLocMesh, std::vector<double> _CoarseGlobMesh);
-};
-} // namespace dart
-#endif //WSolutionInterp_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index 916764c256db144782179057c8a895b2026c57ca..43d837ced7956634426f6aff0a01b2c9b4210b04 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -1,6 +1,7 @@
#include "wTimeSolver.h"
#include "wBLRegion.h"
#include "wViscFluxes.h"
+#include "wInvBnd.h"
#include <Eigen/Dense>
#include <iostream>
#include <vector>
@@ -10,262 +11,201 @@ using namespace Eigen;
using namespace tbox;
using namespace dart;
-
TimeSolver::TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac)
{
- CFL0 = _CFL0;
- maxIt = _maxIt;
- tol = _tol;
- itFrozenJac0 = _itFrozenJac;
- Minf = std::max(_Minf, 0.1);
+ CFL0 = _CFL0;
+ maxIt = _maxIt;
+ tol = _tol;
+ itFrozenJac0 = _itFrozenJac;
+ Minf = std::max(_Minf, 0.1);
- SysMatrix = new ViscFluxes(Re);
-
-} // end TimeSolver
+ SysMatrix = new ViscFluxes(Re);
+} // end TimeSolver
TimeSolver::~TimeSolver()
-{
- delete SysMatrix;
- std::cout << "~TimeSolver()\n";
+{
+ delete SysMatrix;
+ std::cout << "~TimeSolver()\n";
} // end ~TimeSolver
void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
{
- unsigned int nVar = bl->GetnVar();
- if (initSln){
-
- 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->GetVtExt(iPoint);
- bl->U[iPoint * nVar + 4] = bl->U[(iPoint-1) * nVar + 4];
-
- bl->Ret[iPoint] = bl->Ret[iPoint - 1];
- bl->Cd[iPoint] = bl->Cd[iPoint - 1];
- bl->Cf[iPoint] = bl->Cf[iPoint - 1];
- bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
- bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
- bl->Hk[iPoint] = bl->Hk[iPoint - 1];
- bl->Cteq[iPoint] = bl->Cteq[iPoint - 1];
- bl->Us[iPoint] = bl->Us[iPoint - 1];
- bl->Delta[iPoint] = bl->Delta[iPoint - 1];
- bl->DeltaStar[iPoint] = bl->DeltaStar[iPoint - 1];
- }
- 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
+ unsigned int nVar = bl->GetnVar();
+ if (initSln)
+ {
+
+ 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 + 4] = bl->U[(iPoint - 1) * nVar + 4];
+
+ bl->Ret[iPoint] = bl->Ret[iPoint - 1];
+ bl->Cd[iPoint] = bl->Cd[iPoint - 1];
+ bl->Cf[iPoint] = bl->Cf[iPoint - 1];
+ bl->Hstar[iPoint] = bl->Hstar[iPoint - 1];
+ bl->Hstar2[iPoint] = bl->Hstar2[iPoint - 1];
+ bl->Hk[iPoint] = bl->Hk[iPoint - 1];
+ bl->Cteq[iPoint] = bl->Cteq[iPoint - 1];
+ bl->Us[iPoint] = bl->Us[iPoint - 1];
+ bl->Delta[iPoint] = bl->Delta[iPoint - 1];
+ bl->DeltaStar[iPoint] = bl->DeltaStar[iPoint - 1];
+ }
+ if (bl->Regime[iPoint] == 1 && bl->U[iPoint * nVar + 4] == 0)
+ {
+ bl->U[iPoint * nVar + 4] = 0.03;
+ } // End if
} // End InitialCondition
int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
{
- unsigned int nVar = bl->GetnVar();
-
- /* Save initial condition */
+ unsigned int nVar = bl->GetnVar();
+
+ /* Save initial condition */
- std::vector<double> Sln0;
- for (size_t i=0; i<nVar; ++i){
- Sln0.push_back(bl->U[iPoint*nVar+i]);
- } // End for
+ std::vector<double> Sln0;
+ for (size_t i = 0; i < nVar; ++i)
+ {
+ Sln0.push_back(bl->U[iPoint * nVar + i]);
+ } // End for
- /* Initialize time integration variables */
+ /* Initialize time integration variables */
- double dx = bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1);
+ double dx = bl->mesh->GetLoc(iPoint) - bl->mesh->GetLoc(iPoint - 1);
- /* Initial time step */
+ /* Initial time step */
- double CFL = CFL0;
- unsigned int itFrozenJac = itFrozenJac0;
- double timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
- Matrix<double, 5, 5> IMatrix;
- IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
+ double CFL = CFL0;
+ unsigned int itFrozenJac = itFrozenJac0;
+ double timeStep = SetTimeStep(CFL, dx, bl->invBC->GetVt(iPoint));
+ Matrix<double, 5, 5> IMatrix;
+ IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
- /* Initial system residuals */
+ /* Initial system residuals */
- Vector<double, 5> SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
- double normSysRes0 = SysRes.norm();
- double normSysRes = normSysRes0;
+ Vector<double, 5> SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
+ double normSysRes0 = SysRes.norm();
+ double normSysRes = normSysRes0;
- Matrix<double, 5, 5> JacMatrix;
- Vector<double, 5> slnIncr;
-
- nErrors = 0; // Number of errors encountered
- unsigned int innerIters = 0; // Inner (non-linear) iterations
+ Matrix<double, 5, 5> JacMatrix;
+ Vector<double, 5> slnIncr;
- while (innerIters < maxIt){
+ nErrors = 0; // Number of errors encountered
+ unsigned int innerIters = 0; // Inner (non-linear) iterations
- if(innerIters % itFrozenJac == 0)
+ while (innerIters < maxIt)
{
- if (nErrors >= 5)
- {
- ResetSolution(iPoint, bl, Sln0);
- return -1; // Convergence failed
- }
- if (std::isnan(CFL))
- {
- CFL = 0.5;
- }
- CFL = ControlIntegration(iPoint, bl, Sln0, CFL);
- itFrozenJac = itFrozenJac0;
- JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
- }
- slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-SysRes);
+ if (innerIters % itFrozenJac == 0)
+ {
+ JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
+ }
- for(size_t k=0; k<nVar; ++k)
- {
- bl->U[iPoint*nVar+k] += slnIncr(k);
- }
+ slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(-SysRes);
- SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
- normSysRes = (SysRes + slnIncr/timeStep).norm();
+ for (size_t k = 0; k < nVar; ++k)
+ {
+ bl->U[iPoint * nVar + k] += slnIncr(k);
+ }
- if (normSysRes <= tol * normSysRes0)
- {
- return 0; /* Successfull exit */
- }
+ SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
+ normSysRes = (SysRes + slnIncr / timeStep).norm();
- // CFL adaptation.
+ if (normSysRes <= tol * normSysRes0)
+ {
+ return 0; /* Successfull exit */
+ }
- CFL = CFL0* pow(normSysRes0/normSysRes, 0.7);
- CFL = std::max(CFL, 0.3);
- timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
- IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
+ /* CFL adaptation. */
- innerIters++;
- } // End while (innerIters < maxIt)
+ CFL = std::max(CFL0 * pow(normSysRes0 / normSysRes, 0.7), 0.1);
+ timeStep = SetTimeStep(CFL, dx, bl->invBC->GetVt(iPoint));
+ IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
- for (size_t k=0; k<nVar; ++k)
- {
- if (std::isnan(bl->U[iPoint*nVar + k]))
- {
- ResetSolution(iPoint, bl, Sln0);
- return innerIters; // New CFL.
- }
- }
- if (normSysRes <= 1e-4 * normSysRes0)
- {
- return 0;
- }
- else
- {
- ResetSolution(iPoint, bl, Sln0, true);
- return innerIters;
- }
-
+ innerIters++;
+ } // End while (innerIters < maxIt)
- ResetSolution(iPoint, bl, Sln0);
- return innerIters;
+ ResetSolution(iPoint, bl, Sln0);
+ return innerIters;
} // End Integration
double TimeSolver::SetTimeStep(double CFL, double dx, double invVel)
{
- /* Speed of sound. */
- double gradU2 = (invVel*invVel);
- double soundSpeed = ComputeSoundSpeed(gradU2);
+ /* Speed of sound. */
+ double gradU2 = (invVel * invVel);
+ double soundSpeed = ComputeSoundSpeed(gradU2);
- /* Velocity of the fastest travelling wave */
- double advVel = soundSpeed + invVel;
+ /* Velocity of the fastest travelling wave */
+ double advVel = soundSpeed + invVel;
- /* Time step computation */
- return CFL*dx/advVel;
+ /* Time step computation */
+ return CFL * dx / advVel;
}
double TimeSolver::ComputeSoundSpeed(double gradU2)
{
- return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
+ return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
}
-double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, double CFLIn){
-
- /* Check if there are NaN in the solution */
- unsigned int nVar = bl->GetnVar();
-
- /* Check for eventual errors */
+void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint)
+{
- bl->U[iPoint*nVar + 0] = std::max(bl->U[iPoint*nVar + 0], 1e-6);
- bl->U[iPoint*nVar + 1] = std::max(bl->U[iPoint*nVar + 1], 1.0005);
+ unsigned int nVar = bl->GetnVar();
- if (bl->U[iPoint*nVar + 3] <= 0)
- {
- ResetSolution(iPoint, bl, Sln0, true);
- nErrors+=1;
- return 0.2;
- }
+ /* Reset solution */
- for (size_t k=0; k<nVar; ++k)
- {
- if (std::isnan(bl->U[iPoint*nVar + k]))
+ if (usePrevPoint)
{
- ResetSolution(iPoint, bl, Sln0, true);
- //itFrozenJac0 = 1; // Impose continuous Jacobian evaluation.
- nErrors+=1;
- return 0.1; // New CFL.
- }
- }
- return CFLIn; // Keep CFL.
-}
-
-void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint)
-{
-
- unsigned int nVar = bl->GetnVar();
+ for (size_t k = 0; k < nVar; ++k)
+ {
+ bl->U[iPoint * nVar + k] = bl->U[(iPoint - 1) * nVar + k];
+ } // End for
+ } // End if (usePrevPoint)
+ else
+ {
+ for (size_t k = 0; k < nVar; ++k)
+ {
+ bl->U[iPoint * nVar + k] = Sln0[k];
+ } // End for
+ } // End else (usePrevPoint)
- /* Reset solution */
+ /* Reset closures */
- if (usePrevPoint)
- {
- for(size_t k=0; k<nVar; ++k)
- {
- bl->U[iPoint*nVar + k] = bl->U[(iPoint-1)*nVar + k];
- } // End for
- } // End if (usePrevPoint)
- else
- {
- for(size_t k=0; k<nVar; ++k)
+ if (bl->Regime[iPoint] == 0)
{
- bl->U[iPoint*nVar + k] = Sln0[k];
- } // End for
- } // End else (usePrevPoint)
-
- /* 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
+ 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);
+ 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->invBC->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 ResetSolution
+double TimeSolver::GetCFL0() { return CFL0; }
+void TimeSolver::SetCFL0(double _CFL0) { CFL0 = _CFL0; }
+unsigned long TimeSolver::GetmaxIt() { return maxIt; }
+void TimeSolver::SetmaxIt(unsigned long _maxIt) { maxIt = _maxIt; }
+unsigned int TimeSolver::GetitFrozenJac() { return itFrozenJac0; }
+void TimeSolver::SetitFrozenJac(unsigned int _itFrozenJac) { itFrozenJac0 = _itFrozenJac; }
-double TimeSolver::GetCFL0(){return CFL0;}
-void TimeSolver::SetCFL0(double _CFL0){CFL0 = _CFL0;}
-unsigned long TimeSolver::GetmaxIt(){return maxIt;}
-void TimeSolver::SetmaxIt(unsigned long _maxIt){maxIt = _maxIt;}
-unsigned int TimeSolver::GetitFrozenJac(){return itFrozenJac0;}
-void TimeSolver::SetitFrozenJac(unsigned int _itFrozenJac){itFrozenJac0 = _itFrozenJac;}
-
-void TimeSolver::SetinitSln(bool _initSln){initSln = _initSln;}
\ No newline at end of file
+void TimeSolver::SetinitSln(bool _initSln) { initSln = _initSln; }
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index dcf8baebf38a1481161c49548622792e7b9e34ec..ff0d6900e45609594c4d5704688a2cc04eba6c31 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -13,37 +13,34 @@ class DART_API TimeSolver
{
private:
+ double CFL0;
+ unsigned long maxIt;
+ double tol;
+ unsigned int itFrozenJac0;
+ bool initSln;
+ double gamma = 1.4;
+ double Minf;
+ unsigned int nErrors;
- double CFL0;
- unsigned long maxIt;
- double tol;
- unsigned int itFrozenJac0;
- bool initSln;
- double gamma = 1.4;
- double Minf;
- unsigned int nErrors;
-
- ViscFluxes *SysMatrix;
+ ViscFluxes *SysMatrix;
public:
- TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 50, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
- ~TimeSolver();
- void InitialCondition(size_t iPoint, BLRegion *bl);
- int Integration(size_t iPoint, BLRegion *bl);
- double GetCFL0();
- void SetCFL0(double _CFL0);
- unsigned long GetmaxIt();
- void SetmaxIt(unsigned long _maxIt);
- unsigned int GetitFrozenJac();
- void SetitFrozenJac(unsigned int _itFrozenJac);
- void SetinitSln(bool _initSln);
+ TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _maxIt = 50, double _tol = 1e-8, unsigned int _itFrozenJac = 5);
+ ~TimeSolver();
+ void InitialCondition(size_t iPoint, BLRegion *bl);
+ int Integration(size_t iPoint, BLRegion *bl);
+ double GetCFL0();
+ void SetCFL0(double _CFL0);
+ unsigned long GetmaxIt();
+ void SetmaxIt(unsigned long _maxIt);
+ unsigned int GetitFrozenJac();
+ void SetitFrozenJac(unsigned int _itFrozenJac);
+ void SetinitSln(bool _initSln);
private:
- double SetTimeStep(double CFL, double dx, double invVel);
- double ComputeSoundSpeed(double gradU2);
- double ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, double CFLIn);
- void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint=false);
-
+ double SetTimeStep(double CFL, double dx, double invVel);
+ double ComputeSoundSpeed(double gradU2);
+ void ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> Sln0, bool usePrevPoint = false);
};
} // namespace dart
#endif //WTimeSolver_H
\ No newline at end of file
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 1c68fee80b27ae0c509398a59d98aa6be947b225..8c02dd5d2d92ae97170029fae39e9ad210f4415a 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -1,7 +1,8 @@
#define _USE_MATH_DEFINES
-
+
#include "wViscFluxes.h"
#include "wBLRegion.h"
+#include "wInvBnd.h"
#include <iostream>
#include <Eigen/Dense>
#include <cmath>
@@ -11,249 +12,286 @@
using namespace dart;
using namespace Eigen;
-
ViscFluxes::ViscFluxes(double _Re)
{
- Re = _Re;
-
-} // end ViscFluxes
+ Re = _Re;
+} // end ViscFluxes
ViscFluxes::~ViscFluxes()
-{
- std::cout << "~ViscFluxes()\n";
+{
+ std::cout << "~ViscFluxes()\n";
} // end ~ViscFluxes
Vector<double, 5> ViscFluxes::ComputeFluxes(size_t iPoint, BLRegion *bl)
{
- unsigned int nVar = bl->GetnVar();
- std::vector<double> u;
- for (size_t k=0; k<nVar; ++k){
- u.push_back(bl->U[iPoint*nVar + k]);
- } // End for
- return BLlaws(iPoint, bl, u);
+ unsigned int nVar = bl->GetnVar();
+ std::vector<double> u;
+ for (size_t k = 0; k < nVar; ++k)
+ {
+ u.push_back(bl->U[iPoint * nVar + k]);
+ } // End for
+ return BLlaws(iPoint, bl, u);
}
Matrix<double, 5, 5> ViscFluxes::ComputeJacobian(size_t iPoint, BLRegion *bl, Vector<double, 5> &sysRes, double eps)
{
- unsigned int nVar = bl->GetnVar();
- Matrix<double, 5, 5> JacMatrix;
- std::vector<double> uUp;
- for (size_t k=0; k<nVar; ++k){
- uUp.push_back(bl->U[iPoint*nVar + k]);
- } // End for
-
- double varSave;
- for (size_t iVar=0; iVar<nVar; ++iVar){
- varSave = uUp[iVar];
- uUp[iVar] += eps;
-
- JacMatrix.col(iVar) = (BLlaws(iPoint, bl, uUp) - sysRes) / eps;
-
- uUp[iVar] = varSave;
- }
- return JacMatrix;
+ unsigned int nVar = bl->GetnVar();
+ Matrix<double, 5, 5> JacMatrix;
+ std::vector<double> uUp;
+ for (size_t k = 0; k < nVar; ++k)
+ {
+ uUp.push_back(bl->U[iPoint * nVar + k]);
+ } // End for
+
+ double varSave;
+ for (size_t iVar = 0; iVar < nVar; ++iVar)
+ {
+ varSave = uUp[iVar];
+ uUp[iVar] += eps;
+
+ JacMatrix.col(iVar) = (BLlaws(iPoint, bl, uUp) - sysRes) / eps;
+
+ uUp[iVar] = varSave;
+ }
+ return JacMatrix;
}
Vector<double, 5> ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
{
- unsigned int nVar = bl->GetnVar();
-
- Matrix<double, 5, 5> timeMatrix;
- Vector<double, 5> spaceVector;
-
- double dissipRatio;
- if (bl->name == "wake"){
- dissipRatio = 0.9;
- } // End if
- else{
- dissipRatio = 1;
- } // End else
-
- 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 */
-
- /* Boundary layer closure relations */
-
- if (bl->Regime[iPoint] == 0){
- /* Laminar closure relations */
- std::vector<double> lamParam = bl->closSolver->LaminarClosures(u[0], u[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];
+ unsigned int nVar = bl->GetnVar();
+
+ Matrix<double, 5, 5> timeMatrix;
+ Vector<double, 5> spaceVector;
+
+ double dissipRatio;
+ if (bl->name == "wake")
+ {
+ dissipRatio = 0.9;
+ } // End if
+ else
+ {
+ dissipRatio = 1;
+ } // End else
+
+ bl->Ret[iPoint] = u[3] * u[0] * Re; /* Reynolds number based on theta */
+ bl->DeltaStar[iPoint] = u[0] * u[1]; /* Displacement thickness */
+
+ /* Boundary layer closure relations */
+
+ 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);
+ 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)
+ {
+ /* 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);
+ 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
+ else
+ {
+ std::cout << "Wrong regime\n"
+ << std::endl;
+ } // End else
+
+ /* Gradients computation */
+ double dx = (bl->mesh->GetLoc(iPoint) - bl->mesh->GetLoc(iPoint - 1));
+ double dxExt = (bl->mesh->GetExt(iPoint) - bl->mesh->GetExt(iPoint - 1));
+ double dT_dx = (u[0] - bl->U[(iPoint - 1) * nVar + 0]) / dx;
+ double dH_dx = (u[1] - bl->U[(iPoint - 1) * nVar + 1]) / dx;
+ double due_dx = (u[3] - bl->U[(iPoint - 1) * nVar + 3]) / dx;
+ 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 c = 4 / (M_PI * dx);
+ double cExt = 4 / (M_PI * dxExt);
+
+ /* Amplification routine */
+ double dN_dx;
+ double Ax;
+
+ if (bl->xtrF != -1 && bl->Regime[iPoint] == 0)
+ {
+ Ax = AmplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
+ dN_dx = (u[2] - bl->U[(iPoint - 1) * nVar + 2]) / dx;
+
+ } // End if
+ else
+ {
+ Ax = 0;
+ dN_dx = 0;
+ } // End else
+
+ double Mea = bl->invBC->GetMe(iPoint);
+
+ double Hstara = bl->Hstar[iPoint];
+ double Hstar2a = bl->Hstar2[iPoint];
+ double Hka = bl->Hk[iPoint];
+ double Cfa = bl->Cf[iPoint];
+ double Cda = bl->Cd[iPoint];
+ double deltaa = bl->Delta[iPoint];
+ double Cteqa = bl->Cteq[iPoint];
+ double Usa = bl->Us[iPoint];
+
+ /* if (bl->name == "upper" && iPoint == 168)
+ {
+ std::cout << "Regime " << bl->Regime[iPoint] << " \n"
+ << "Mea " << Mea << " \n"
+ << "Ret " << bl->Ret[iPoint] << " \n"
+ << "Hstar " << Hstara << " \n"
+ << "Hstar2 " << Hstar2a << " \n"
+ << "Hk " << Hka << " \n"
+ << "Cf " << Cfa << " \n"
+ << "Cd " << Cda << " \n"
+ << "delta " << deltaa << " \n"
+ << "Cteq " << Cteqa << " \n"
+ << "Us " << Usa << std::endl;
+ } */
+
+ /* Space part */
+
+ spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - Cfa / 2;
+ spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2 * Cda + Hstara * Cfa / 2;
+ spaceVector(2) = dN_dx - Ax;
+ spaceVector(3) = due_dx - c * (u[1] * dT_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
+
+ if (bl->Regime[iPoint] == 1)
+ {
+ spaceVector(4) = (2 * deltaa / u[4]) * dCt_dx - 5.6 * (Cteqa - u[4] * dissipRatio) - 2 * deltaa * (4 / (3 * u[1] * u[0]) * (Cfa / 2 - (((Hka - 1) / (6.7 * Hka * dissipRatio)) * ((Hka - 1) / (6.7 * Hka * dissipRatio)))) - 1 / u[3] * due_dx);
} // End if
-
- 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->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
- else{
- std::cout << "Wrong regime\n" << std::endl;
- } // End else
-
- /* Gradients computation */
- double dx = (bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1));
- double dxExt = (bl->GetXxExt(iPoint) - bl->GetXxExt(iPoint-1));
- double dT_dx = (u[0] - bl->U[(iPoint-1)*nVar + 0])/dx;
- double dH_dx = (u[1] - bl->U[(iPoint-1)*nVar + 1])/dx;
- double due_dx = (u[3] - bl->U[(iPoint-1)*nVar + 3])/dx;
- 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->GetVtExt(iPoint) - bl->GetVtExt(iPoint-1))/dxExt;
- double ddeltaStarExt_dx = (bl->GetDeltaStarExt(iPoint) - bl->GetDeltaStarExt(iPoint-1))/dxExt;
-
- double c=4/(M_PI*dx);
- double cExt=4/(M_PI*dxExt);
-
- /* Amplification routine */
- double dN_dx;
- double Ax;
-
- if (bl->xtrF!=-1 && bl->Regime[iPoint] == 0){
- Ax = AmplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
- dN_dx = (u[2] - bl->U[(iPoint-1)*nVar + 2])/dx;
-
- } // End if
- else{
- Ax = 0;
- dN_dx = 0;
- } // End else
-
- double Mea = bl->GetMe(iPoint);
-
- double Hstara = bl->Hstar[iPoint];
- double Hstar2a = bl->Hstar2[iPoint];
- double Hka = bl->Hk[iPoint];
- double Cfa = bl->Cf[iPoint];
- double Cda = bl->Cd[iPoint];
- double deltaa = bl->Delta[iPoint];
- double Cteqa = bl->Cteq[iPoint];
- double Usa = bl->Us[iPoint];
-
- /* Space part */
-
- spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0]/u[3]) * due_dx - Cfa/2;
- spaceVector(1) = u[0] * dHstar_dx + (2*Hstar2a + Hstara * (1-u[1])) * u[0]/u[3] * due_dx - 2*Cda + Hstara * Cfa/2;
- spaceVector(2) = dN_dx - Ax;
- spaceVector(3) = due_dx - c*(u[1] * dT_dx + u[0]*dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
-
- if (bl->Regime[iPoint] == 1){
- spaceVector(4) = (2*deltaa/u[4]) * dCt_dx - 5.6*(Cteqa - u[4]*dissipRatio)
- - 2*deltaa*(4/(3*u[1]*u[0]) * (Cfa/2 - (((Hka-1)/(6.7*Hka * dissipRatio)) * ((Hka-1)/(6.7*Hka * dissipRatio)))) - 1/u[3] * due_dx);
- } // End if
- else{
- spaceVector(4) = 0.0;
- }
-
- /* Time part */
-
- timeMatrix(0,0) = u[1]/u[3];
- timeMatrix(0,1) = u[0]/u[3];
- timeMatrix(0,2) = 0.0;
- timeMatrix(0,3) = u[0]*u[1]/(u[3]*u[3]);
- timeMatrix(0,4) = 0.0;
-
- timeMatrix(1,0) = (1+u[1]*(1-Hstara))/u[3];
- timeMatrix(1,1) = (1-Hstara)*u[0]/u[3];
- timeMatrix(1,2) = 0.0;
- timeMatrix(1,3) = (2-Hstara*u[1])*u[0]/(u[3]*u[3]);
- timeMatrix(1,4) = 0.0;
-
- timeMatrix(2,0) = 0.0;
- timeMatrix(2,1) = 0.0;
- timeMatrix(2,2) = 1.0;
- timeMatrix(2,3) = 0.0;
- timeMatrix(2,4) = 0.0;
-
- timeMatrix(3,0) = -c*u[1];
- timeMatrix(3,1) = -c*u[0];
- timeMatrix(3,2) = 0.0;
- timeMatrix(3,3) = 1.0;
- timeMatrix(3,4) = 0.0;
-
- timeMatrix(4,0) = 0.0;
- timeMatrix(4,1) = 0.0;
- timeMatrix(4,2) = 0.0;
-
- if (bl->Regime[iPoint]==1){
- timeMatrix(4,3) = 2*deltaa/(Usa*u[3]*u[3]);
- timeMatrix(4,4) = 2*deltaa/(u[3]*u[4]*Usa);
- }
- else{
- timeMatrix(4,3) = 0.0;
- timeMatrix(4,4) = 1.0;
- }
-
- return timeMatrix.partialPivLu().solve(spaceVector);
+ else
+ {
+ spaceVector(4) = 0.0;
+ }
+
+ /* Time part */
+
+ timeMatrix(0, 0) = u[1] / u[3];
+ timeMatrix(0, 1) = u[0] / u[3];
+ timeMatrix(0, 2) = 0.0;
+ timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
+ timeMatrix(0, 4) = 0.0;
+
+ timeMatrix(1, 0) = (1 + u[1] * (1 - Hstara)) / u[3];
+ timeMatrix(1, 1) = (1 - Hstara) * u[0] / u[3];
+ timeMatrix(1, 2) = 0.0;
+ timeMatrix(1, 3) = (2 - Hstara * u[1]) * u[0] / (u[3] * u[3]);
+ timeMatrix(1, 4) = 0.0;
+
+ timeMatrix(2, 0) = 0.0;
+ timeMatrix(2, 1) = 0.0;
+ timeMatrix(2, 2) = 1.0;
+ timeMatrix(2, 3) = 0.0;
+ timeMatrix(2, 4) = 0.0;
+
+ timeMatrix(3, 0) = -c * u[1];
+ timeMatrix(3, 1) = -c * u[0];
+ timeMatrix(3, 2) = 0.0;
+ timeMatrix(3, 3) = 1.0;
+ timeMatrix(3, 4) = 0.0;
+
+ timeMatrix(4, 0) = 0.0;
+ timeMatrix(4, 1) = 0.0;
+ timeMatrix(4, 2) = 0.0;
+
+ if (bl->Regime[iPoint] == 1)
+ {
+ timeMatrix(4, 3) = 2 * deltaa / (Usa * u[3] * u[3]);
+ timeMatrix(4, 4) = 2 * deltaa / (u[3] * u[4] * Usa);
+ }
+ else
+ {
+ timeMatrix(4, 3) = 0.0;
+ timeMatrix(4, 4) = 1.0;
+ }
+
+ return timeMatrix.partialPivLu().solve(spaceVector);
}
double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta)
{
- double Ax;
-
- double Dgr = 0.08;
- double Hk1 = 3.5;
- double Hk2 = 4;
- double Hmi = (1/(Hk-1));
- double logRet_crit = 2.492*pow(1/(Hk-1),0.43) + 0.7*(tanh(14*Hmi-9.24)+1.0); // Value at which the amplification starts to grow
- if (Ret <=0){
- Ret = 1;
- }
- if (log10(Ret) < logRet_crit - Dgr){
- Ax = 0;
- }
- else{
- double Rnorm = (log10(Ret) - (logRet_crit-Dgr))/(2*Dgr);
- double Rfac;
- if(Rnorm <=0){
- Rfac = 0;
- }
- else if(Rnorm >= 1){
- Rfac = 1;
- }
- else{
- Rfac = 3*Rnorm*Rnorm - 2*Rnorm*Rnorm*Rnorm;
- }
- // Normal envelope amplification rate
- double m = -0.05+2.7*Hmi-5.5*Hmi*Hmi+3*Hmi*Hmi*Hmi+0.1*exp(-20*Hmi);
- double arg = 3.87*Hmi-2.52;
- double dndRet = 0.028*(Hk-1)-0.0345*exp(-(arg*arg));
- Ax = (m*dndRet/theta)*Rfac;
- // Set correction for separated profiles
- if(Hk > Hk1){
- double Hnorm = (Hk - Hk1)/(Hk2-Hk1);
- double Hfac;
- if(Hnorm >=1){
- Hfac = 1;
- }
- else{
- Hfac = 3*Hnorm*Hnorm - 2*Hnorm*Hnorm*Hnorm;
- }
- double Ax1 = Ax;
- double Gro = 0.3+0.35*exp(-0.15*(Hk-5));
- double Tnr = tanh(1.2*(log10(Ret)-Gro));
- double Ax2 = (0.086*Tnr - 0.25/pow(Hk-1,1.5))/theta;
- if(Ax2 < 0){
- Ax2 = 0;
- }
- Ax = Hfac*Ax2 + (1-Hfac)*Ax1;
- Ax = std::max(Ax, 0.0);
- }
- }
- return Ax;
+ double Ax;
+
+ double Dgr = 0.08;
+ double Hk1 = 3.5;
+ double Hk2 = 4;
+ double Hmi = (1 / (Hk - 1));
+ double logRet_crit = 2.492 * pow(1 / (Hk - 1), 0.43) + 0.7 * (tanh(14 * Hmi - 9.24) + 1.0); // Value at which the amplification starts to grow
+ if (Ret <= 0)
+ {
+ Ret = 1;
+ }
+ if (log10(Ret) < logRet_crit - Dgr)
+ {
+ Ax = 0;
+ }
+ else
+ {
+ double Rnorm = (log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
+ double Rfac;
+ if (Rnorm <= 0)
+ {
+ Rfac = 0;
+ }
+ else if (Rnorm >= 1)
+ {
+ Rfac = 1;
+ }
+ else
+ {
+ Rfac = 3 * Rnorm * Rnorm - 2 * Rnorm * Rnorm * Rnorm;
+ }
+ // Normal envelope amplification rate
+ double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi + 0.1 * exp(-20 * Hmi);
+ double arg = 3.87 * Hmi - 2.52;
+ double dndRet = 0.028 * (Hk - 1) - 0.0345 * exp(-(arg * arg));
+ Ax = (m * dndRet / theta) * Rfac;
+ // Set correction for separated profiles
+ if (Hk > Hk1)
+ {
+ double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
+ double Hfac;
+ if (Hnorm >= 1)
+ {
+ Hfac = 1;
+ }
+ else
+ {
+ Hfac = 3 * Hnorm * Hnorm - 2 * Hnorm * Hnorm * Hnorm;
+ }
+ double Ax1 = Ax;
+ double Gro = 0.3 + 0.35 * exp(-0.15 * (Hk - 5));
+ double Tnr = tanh(1.2 * (log10(Ret) - Gro));
+ double Ax2 = (0.086 * Tnr - 0.25 / pow(Hk - 1, 1.5)) / theta;
+ if (Ax2 < 0)
+ {
+ Ax2 = 0;
+ }
+ Ax = Hfac * Ax2 + (1 - Hfac) * Ax1;
+ Ax = std::max(Ax, 0.0);
+ }
+ }
+ return Ax;
}
\ No newline at end of file
diff --git a/dart/src/wViscFluxes.h b/dart/src/wViscFluxes.h
index f79f87ea3fa05c2d605744f44484a5f3c36c1eaf..403c388e8ce4fd9d6fc53070dcfe40233a558ef8 100644
--- a/dart/src/wViscFluxes.h
+++ b/dart/src/wViscFluxes.h
@@ -13,17 +13,17 @@ namespace dart
class DART_API ViscFluxes
{
private:
- double Re;
-public:
- ViscFluxes(double _Re);
- ~ViscFluxes();
- Eigen::Vector<double, 5> ComputeFluxes(size_t iPoint, BLRegion *bl);
- Eigen::Matrix<double, 5, 5> ComputeJacobian(size_t iPoint, BLRegion *bl, Eigen::Vector<double, 5> &sysRes, double eps);
-private:
- Eigen::Vector<double, 5> BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u);
- double AmplificationRoutine(double hk, double ret, double theta);
+ double Re;
+public:
+ ViscFluxes(double _Re);
+ ~ViscFluxes();
+ Eigen::Vector<double, 5> ComputeFluxes(size_t iPoint, BLRegion *bl);
+ Eigen::Matrix<double, 5, 5> ComputeJacobian(size_t iPoint, BLRegion *bl, Eigen::Vector<double, 5> &sysRes, double eps);
+private:
+ Eigen::Vector<double, 5> BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u);
+ double AmplificationRoutine(double hk, double ret, double theta);
};
} // namespace dart
#endif //WViscFluxes_H
\ No newline at end of file
diff --git a/dart/src/wViscMesh.cpp b/dart/src/wViscMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b1476011d42d70a1fc27fbf15ae61cbf871978ed
--- /dev/null
+++ b/dart/src/wViscMesh.cpp
@@ -0,0 +1,64 @@
+/* ----------------------------------------------------------------------------------- */
+/* */
+/* Coupled integral boundary layer solver */
+/* */
+/* Université de Liège */
+/* February 2022 */
+/* Author: Paul Dechamps */
+/* v0.1 */
+/* ----------------------------------------------------------------------------------- */
+
+/* --- Project includes -------------------------------------------------------------- */
+
+#include "wViscMesh.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+/* --- Namespaces -------------------------------------------------------------------- */
+
+using namespace dart;
+
+/* ----------------------------------------------------------------------------------- */
+
+ViscMesh::ViscMesh(){}
+
+ViscMesh::~ViscMesh(){}
+
+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();
+ nElm = nMarkers - 1;
+ Glob.resize(nMarkers);
+ for (size_t iPoint=0; iPoint<Glob.size(); ++iPoint)
+ {
+ Glob[iPoint].resize(3, 0);
+ }
+
+ for (size_t iPoint=0; iPoint<x.size(); ++iPoint)
+ {
+ Glob[iPoint][0] = x[iPoint];
+ Glob[iPoint][1] = y[iPoint];
+ Glob[iPoint][2] = z[iPoint];
+ }
+}
+
+void ViscMesh::ComputeBLcoord()
+{
+ Loc.resize(nMarkers, 0);
+ std::vector<double> dx(nElm,0);
+ alpha.resize(nElm, 0);
+
+ for (size_t iPoint=0; iPoint<nElm; ++iPoint)
+ {
+ alpha[iPoint] = std::atan2(Glob[iPoint+1][1], Glob[iPoint+1][0]);
+ /* 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]));
+ Loc[iPoint + 1] = Loc[iPoint] + dx[iPoint];
+ }
+}
diff --git a/dart/src/wViscMesh.h b/dart/src/wViscMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..1497faa25ef18041e7799566685f30f7ec366738
--- /dev/null
+++ b/dart/src/wViscMesh.h
@@ -0,0 +1,41 @@
+#ifndef WVISCMESH_H
+#define WVISCMESH_H
+
+#include "dart.h"
+#include "wBLRegion.h"
+#include <vector>
+
+namespace dart
+{
+
+class DART_API ViscMesh
+{
+ private:
+ std::vector<std::vector<double>> Glob;
+ std::vector<double> Loc, /* BL (local) coordinates */
+ Ext, /* Coordinates at the edge of the BL */
+ dx,
+ alpha; /* Angle of the element for the rotation matrix */
+
+ unsigned int nMarkers;
+ unsigned int nElm;
+
+ public:
+ 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 GetExt(size_t iPoint) const {return Ext[iPoint];};
+
+ void SetGlob(std::vector<double> x, std::vector<double> y, std::vector<double> z);
+ void SetExt(std::vector<double> ExtM) {Ext = ExtM;};
+
+ void ComputeBLcoord();
+
+};
+} // namespace dart
+#endif //WVISCMESH_H
\ No newline at end of file
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index ebdc62d2bb103aa90dcb57f5bf1fe8c2bc59d272..3028c974e4b9c56ca4d89d5373c1912bfef8b2d3 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -19,7 +19,6 @@
/* --- Namespaces -------------------------------------------------------------------- */
-using namespace tbox;
using namespace dart;
/* ----------------------------------------------------------------------------------- */
@@ -27,56 +26,49 @@ using namespace dart;
/**
* @brief Viscous solver and boundary layer(s) initialization
*/
-ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor)
+ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int nSections)
{
- PrintBanner();
+ PrintBanner();
- /* Freestream parameters */
+ /* Freestream parameters */
- Re = _Re; /* Freestream Reynolds number */
- Minf = _Minf; /* Freestream Mach number */
+ Re = _Re; /* Freestream Reynolds number */
- /* User input numerical parameters */
-
- CFL0 = _CFL0; /* Starting CFL number */
+ CFL0 = _CFL0;
- std::cout << "--- Viscous solver setup ---\n"
- << "Reynolds number: " << Re << " \n"
- << "Freestream Mach number: " << Minf << " \n"
- << "Initial CFL number: " << CFL0 << " \n"
- << std::endl;
-
- /* Initialzing boundary layer */
-
- for (size_t i=0; i<3; ++i)
- {
- BLRegion *region = new BLRegion();
- bl.push_back(region);
- } // end for
-
- bl[0]->name = "upper";
- bl[1]->name = "lower";
- bl[2]->name = "wake";
+ /* Initialzing boundary layer */
+ Sections.resize(nSections);
+ for (size_t iSec=0; iSec < nSections; ++iSec)
+ {
+ for (size_t i = 0; i < 3; ++i)
+ {
+ BLRegion *region = new BLRegion();
+ Sections[iSec].push_back(region);
+ }
+ Sections[iSec][0]->name = "upper";
+ Sections[iSec][1]->name = "lower";
+ Sections[iSec][2]->name = "wake";
+ }
- /* Temporal solver initialization */
+ /* Temporal solver initialization */
- tSolver = new TimeSolver(CFL0, Minf, Re);
+ tSolver = new TimeSolver(_CFL0, _Minf, Re);
- /* Mesh conversion setup */
+ /* Stagnation point movements monitoring and control */
- if (meshFactor > 1)
- {
- mapMesh = true;
- }
- meshConverter = new SolutionInterp(meshFactor);
+ nbElmUpper.resize(Sections.size(), 0);
+ stagPtMvmt.resize(Sections.size(), false);
- /* Stagnation point movements monitoring and control */
+ /* User input numerical parameters */
- nbElmUpper = 0;
- stagPtMvmt = false;
+ std::cout << "--- Viscous solver setup ---\n"
+ << "Reynolds number: " << Re << " \n"
+ << "Freestream Mach number: " << _Minf << " \n"
+ << "Initial CFL number: " << CFL0 << " \n"
+ << std::endl;
- std::cout << "Boundary layer initialized" << std::endl;
+ std::cout << "Boundary layer initialized" << std::endl;
} // end ViscSolver
@@ -84,17 +76,21 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int mesh
* @brief Viscous solver and subsequent dependancies destruction
*/
ViscSolver::~ViscSolver()
-{
- for (size_t i=0; i<bl.size(); ++i){
- delete bl[i];
- } // end for
+{
+ for (size_t iSec=0; iSec<Sections.size(); ++iSec)
+ {
+ for (size_t i = 0; i < Sections[iSec].size(); ++i)
+ {
+ delete Sections[iSec][i];
+ }
+ }
- delete tSolver;
- delete meshConverter;
+ delete tSolver;
+ delete meshConverter;
- PrintBanner();
+ PrintBanner();
- std::cout << "~ViscSolver()\n";
+ std::cout << "~ViscSolver()\n";
} // end ~ViscSolver
/**
@@ -102,44 +98,44 @@ ViscSolver::~ViscSolver()
*/
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;
- }
+ 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)
+ 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);
- }
+ 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);
+ }
}
/**
@@ -148,164 +144,198 @@ void ViscSolver::SendInvBC(size_t region, std::vector<double> _Ue,
*/
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);
- }
+ 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.
*/
-int ViscSolver::Run(unsigned int couplIter){
-
- /* Prepare time integration */
-
- if ((couplIter == 0) || (nbElmUpper != bl[0]->GetnMarkers()) || (stagPtMvmt)){
- tSolver->SetCFL0(1);
- tSolver->SetitFrozenJac(1);
- tSolver->SetinitSln(true);
- if (couplIter != 0 && (nbElmUpper != bl[0]->GetnMarkers())){
- stagPtMvmt = true;
- if (printOn){
- std::cout << "Stagnation point moved" << std::endl;
- }
- }
- else{
- stagPtMvmt = false;
- }
- if (printOn){
- std::cout << "Restart solution" << std::endl;
- }
- }
-
- else{
- tSolver->SetCFL0(CFL0);
- tSolver->SetitFrozenJac(5);
- tSolver->SetinitSln(false);
- stagPtMvmt = false;
- if (printOn){
- std::cout << "Updating solution" << std::endl;
- }
- }
-
- nbElmUpper = bl[0]->GetnMarkers();
-
- /* Set boundary condition */
-
- bl[0]->SetStagBC(Re);
- bl[1]->SetStagBC(Re);
-
- /* Loop over the regions */
-
- bool lockTrans;
- int solverExitCode = 0;
- int pointExitCode;
-
- for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion){
+int ViscSolver::Run(unsigned int couplIter)
+{
+ bool lockTrans;
+ int solverExitCode = 0;
+ int pointExitCode;
- /* Reset regime */
- if (iRegion<2){
- for (size_t k=0; k<bl[iRegion]->GetnMarkers(); k++){
- bl[iRegion]->Regime[k]=0;
- }
- }
- else{
- for (size_t k=0; k<bl[iRegion]->GetnMarkers(); k++){
- bl[iRegion]->Regime[k]=1;
- }
- }
- lockTrans = false;
-
- if (printOn){
- std::cout << "- " << bl[iRegion]->name << " region";
- }
-
- if (iRegion == 2){
- /* Impose wake boundary condition */
- std::vector<double> UpperCond;
- std::vector<double> LowerCond;
- 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);
- if (printOn){
- std::cout << "\n" << std::endl;
- }
-
- lockTrans = true;
- }
-
- /* Loop over points */
-
- for (size_t iPoint = 1; iPoint < bl[iRegion]->GetnMarkers(); ++iPoint){
-
- /* Initialize solution at point */
-
- tSolver->InitialCondition(iPoint, bl[iRegion]);
-
- /* Time integration */
+ for (size_t iSec=0; iSec<Sections.size(); ++iSec)
+ {
- pointExitCode = tSolver->Integration(iPoint, bl[iRegion]);
- if (pointExitCode!=0)
- {
- if (printOn){
- std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
+ /* Prepare time integration */
+
+ if ((couplIter == 0) || Sections[iSec][0]->mesh->GetDiffnElm() || (stagPtMvmt[iSec]))
+ {
+ tSolver->SetCFL0(1);
+ tSolver->SetitFrozenJac(1);
+ tSolver->SetinitSln(true);
+ if (couplIter != 0 && (Sections[iSec][0]->mesh->GetDiffnElm()))
+ {
+ stagPtMvmt[iSec] = true;
+ if (printOn)
+ {
+ std::cout << "Stagnation point moved" << std::endl;
+ }
+ }
+ else
+ {
+ stagPtMvmt[iSec] = false;
+ }
+ if (printOn)
+ {
+ std::cout << "Restart solution" << std::endl;
+ }
}
- solverExitCode = -1; // Convergence failed
- }
-
- /* Transition */
- if (!lockTrans)
- {
- CheckWaves(iPoint, bl[iRegion]);
-
- // Free transition.
- if (bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 2] >= bl[iRegion]->GetnCrit()){
- AverageTransition(iPoint, bl[iRegion], 0);
- if (printOn){
- std::cout << ", free transition x/c = " << bl[iRegion]->xtr << ". Marker " << bl[iRegion]->transMarker << std::endl;
- }
- lockTrans = true;
- } // End if
- } // End if (!lockTrans).
-
- } // End for iPoints
- } // End for iRegion
-
- ComputeDrag();
+ else
+ {
+ tSolver->SetCFL0(CFL0);
+ tSolver->SetitFrozenJac(5);
+ tSolver->SetinitSln(false);
+ stagPtMvmt[iSec] = false;
+ if (printOn)
+ {
+ std::cout << "Updating solution" << std::endl;
+ }
+ }
+
+ for (size_t iSec=0; iSec<Sections.size(); ++iSec)
+ {
+
+ nbElmUpper[iSec] = Sections[iSec][0]->GetnMarkers();
+
+ /* Set boundary condition */
+
+ Sections[iSec][0]->xtr = 1.0; /* Upper side initial xtr */
+ Sections[iSec][1]->xtr = 1.0; /* Lower side initial xtr */
+ Sections[iSec][2]->xtr = 0.0; /* Wake initial xtr (full turbulent) */
+
+ Sections[iSec][0]->SetStagBC(Re);
+ Sections[iSec][1]->SetStagBC(Re);
+
+ /* Loop over the regions */
+
+ for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
+ {
+
+ /* Reset regime */
+ if (iRegion < 2)
+ {
+ for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
+ {
+ Sections[iSec][iRegion]->Regime[k] = 0;
+ }
+ }
+ else
+ {
+ for (size_t k = 0; k < Sections[iSec][iRegion]->mesh->GetnMarkers(); k++)
+ {
+ 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)
+ {
+
+ /* Initialize solution at point */
+
+ tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
+
+ /* Time integration */
+
+ pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
+
+ if (pointExitCode != 0)
+ {
+ if (printOn)
+ {
+ std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
+ }
+ solverExitCode = -1; // Convergence failed
+ }
+
+ /* Transition */
+
+ if (!lockTrans)
+ {
+ CheckWaves(iPoint, Sections[iSec][iRegion]);
+
+ // Free transition.
+ if (Sections[iSec][iRegion]->U[iPoint * Sections[iSec][iRegion]->GetnVar() + 2] >= Sections[iSec][iRegion]->GetnCrit())
+ {
+ AverageTransition(iPoint, Sections[iSec][iRegion], 0);
+ if (printOn)
+ {
+ std::cout << ", free transition x/c = "
+ << Sections[iSec][iRegion]->xtr
+ << ". Marker " << Sections[iSec][iRegion]->transMarker
+ << std::endl;
+ }
+ lockTrans = true;
+ } // End if
+ } // End if (!lockTrans).
+
+ } // End for iPoints
+ } // End for iRegion
+ }
+ ComputeDrag();
- for (size_t iRegion=0; iRegion<bl.size(); ++iRegion)
- {
- bl[iRegion]->ComputeBlwVel();
- } // End for iRegion
+ for (size_t iRegion = 0; iRegion < Sections[iSec].size(); ++iRegion)
+ {
+ Sections[iSec][iRegion]->ComputeBlwVel();
+ } // End for iRegion
- return solverExitCode; // exit code
- }
+ return solverExitCode; // exit code
+}
/**
* @brief Check whether instabilities are growing in the laminar boundary layer.
*/
void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
{
- /* Determine if the amplification waves start growing or not. */
+ /* Determine if the amplification waves start growing or not. */
- double logRet_crit = 2.492*pow(1/(bl->Hk[iPoint]-1), 0.43)
- + 0.7*(tanh(14*(1/(bl->Hk[iPoint]-1))-9.24)+1.0);
+ double logRet_crit = 2.492 * pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (tanh(14 * (1 / (bl->Hk[iPoint] - 1)) - 9.24) + 1.0);
- if (bl->Ret[iPoint] > 0)
- { // Avoid log of negative number.
- if (log10(bl->Ret[iPoint]) <= logRet_crit - 0.08){
- bl->U[iPoint*bl->GetnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ if (bl->Ret[iPoint] > 0)
+ { // Avoid log of negative number.
+ if (log10(bl->Ret[iPoint]) <= logRet_crit - 0.08)
+ {
+ bl->U[iPoint * bl->GetnVar() + 2] = 0; // Set N to 0 at that point if waves do not grow.
+ }
}
- }
}
/**
@@ -313,73 +343,78 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
*/
void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
{
- /* Averages solution on transition marker */
-
- unsigned int nVar = bl->GetnVar();
-
- // Save laminar solution.
- std::vector<double> lamSol;
- for (size_t k=0; k<nVar; ++k)
- {
- lamSol.push_back(bl->U[iPoint*nVar + k]);
- }
-
- // Set up turbulent portion boundary condition.
- bl->transMarker = iPoint; /* Save transition marker */
-
- /* Loop over remaining points in the region */
- for (size_t k=0; k<bl->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->GetGlobMarkers(iPoint)-bl->GetGlobMarkers(iPoint-1))/(bl->U[iPoint*nVar + 2]-bl->U[(iPoint-1)*nVar+2]))+bl->GetGlobMarkers(iPoint-1);
-
- /* Percentage of the subsection corresponding to a turbulent flow */
- double avTurb = (bl->GetGlobMarkers(iPoint)-bl->xtr)/(bl->GetGlobMarkers(iPoint)-bl->GetGlobMarkers(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);
- bl->Cteq[iPoint-1] = turbParam1[6];
- bl->U[(iPoint-1)*nVar + 4] = 0.7 * bl->Cteq[iPoint-1];
-
- /* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched line 183) */
- /* These initial guess values do not influence the converged solution */
- bl->U[iPoint*nVar + 4] = bl->U[(iPoint-1)*nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
- bl->U[iPoint*nVar + 1] = 1.515; // Because we expect a significant drop of the shape factor.
-
- /* Solve point in turbulent configuration */
- tSolver->Integration(iPoint, bl);
-
- /* Average both solutions (Now solution @ iPoint is the turbulent solution) */
- bl->U[iPoint*nVar + 0] = avLam * lamSol[0] + avTurb * bl->U[(iPoint)*nVar + 0]; /* Theta */
- bl->U[iPoint*nVar + 1] = avLam * lamSol[1] + avTurb * bl->U[(iPoint)*nVar + 1]; /* H */
- bl->U[iPoint*nVar + 2] = 9; /* N */
- bl->U[iPoint*nVar + 3] = avLam * lamSol[3] + avTurb * bl->U[(iPoint)*nVar + 3]; /* ue */
- 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);
- bl->Hstar[iPoint-1] = lamParam[0];
- bl->Hstar2[iPoint-1] = lamParam[1];
- bl->Hk[iPoint-1] = lamParam[2];
- bl->Cf[iPoint-1] = lamParam[3];
- bl->Cd[iPoint-1] = lamParam[4];
- bl->Delta[iPoint-1] = lamParam[5];
- 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);
- 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];
+ /* Averages solution on transition marker */
+
+ unsigned int nVar = bl->GetnVar();
+
+ // Save laminar solution.
+ std::vector<double> lamSol(nVar);
+ for (size_t k = 0; k < lamSol.size(); ++k)
+ {
+ lamSol[k] = bl->U[iPoint * nVar + k];
+ }
+
+ // Set up turbulent portion boundary condition.
+ bl->transMarker = iPoint; /* Save transition marker */
+
+ /* Loop over remaining points in the region */
+ for (size_t k = 0; k < bl->mesh->GernMarkers() - 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);
+
+ /* 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 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);
+ bl->Cteq[iPoint - 1] = turbParam1[6];
+ bl->U[(iPoint - 1) * nVar + 4] = 0.7 * bl->Cteq[iPoint - 1];
+
+ /* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched line 183) */
+ /* These initial guess values do not influence the converged solution */
+ bl->U[iPoint * nVar + 4] = bl->U[(iPoint - 1) * nVar + 4]; // IC of transition point = turbulent BC (imposed @ previous point).
+ bl->U[iPoint * nVar + 1] = 1.515; // Because we expect a significant drop of the shape factor.
+
+ /* Solve point in turbulent configuration */
+ int exitCode = tSolver->Integration(iPoint, bl);
+
+ if (exitCode != 0)
+ {
+ std::cout << "Warning level 3: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
+ }
+
+ /* Average both solutions (Now solution @ iPoint is the turbulent solution) */
+ bl->U[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->U[(iPoint)*nVar + 0]; /* Theta */
+ bl->U[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->U[(iPoint)*nVar + 1]; /* H */
+ bl->U[iPoint * nVar + 2] = 9; /* N */
+ bl->U[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->U[(iPoint)*nVar + 3]; /* ue */
+ 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);
+ bl->Hstar[iPoint - 1] = lamParam[0];
+ bl->Hstar2[iPoint - 1] = lamParam[1];
+ bl->Hk[iPoint - 1] = lamParam[2];
+ bl->Cf[iPoint - 1] = lamParam[3];
+ bl->Cd[iPoint - 1] = lamParam[4];
+ bl->Delta[iPoint - 1] = lamParam[5];
+ 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);
+ 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];
}
/**
@@ -387,16 +422,16 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
*/
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;
+ std::vector<double> cBlwVel;
+ if (mapMesh)
+ {
+ cBlwVel = meshConverter->BlowingFineToCoarse(bl[iRegion]->Blowing, bl[iRegion]->coarseSize - 1);
+ }
+ else
+ {
+ cBlwVel = bl[iRegion]->Blowing;
+ }
+ return cBlwVel;
}
/**
@@ -404,17 +439,17 @@ std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
*/
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;
+ 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;
}
/**
@@ -425,20 +460,21 @@ std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
*/
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;
+ 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;
}
/**
@@ -451,106 +487,110 @@ std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
void ViscSolver::ComputeDrag()
{
- unsigned int nVar = bl[0]->GetnVar();
- unsigned int lastWkPt = (bl[2]->GetnMarkers() - 1)*nVar;
+ unsigned int nVar = bl[0]->GetnVar();
+ unsigned int lastWkPt = (bl[2]->GetnMarkers() - 1) * nVar;
- /* Normalize friction and dissipation coefficients. */
+ /* Normalize friction and dissipation coefficients. */
- std::vector<double> frictioncoeff_upper;
- for (size_t k=0; k<bl[0]->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){
- frictioncoeff_lower.push_back(bl[1]->U[k*nVar + 3]*bl[1]->U[k*nVar + 3] * bl[1]->Cf[k]);
- }
+ std::vector<double> frictioncoeff_upper;
+ for (size_t k = 0; k < bl[0]->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)
+ {
+ frictioncoeff_lower.push_back(bl[1]->U[k * nVar + 3] * bl[1]->U[k * nVar + 3] * bl[1]->Cf[k]);
+ }
- /* Compute friction drag coefficient. */
+ /* Compute friction drag coefficient. */
- double Cdf_upper = 0.0;
- for (size_t k=1; k<bl[0]->GetnMarkers(); ++k){
- Cdf_upper += (frictioncoeff_upper[k]+frictioncoeff_upper[k-1])*(bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k-1))/2;
- }
- double Cdf_lower = 0.0;
- for (size_t k=1; k<bl[1]->GetnMarkers(); ++k){
- Cdf_lower += (frictioncoeff_lower[k]+frictioncoeff_lower[k-1])*(bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k-1))/2;
+ double Cdf_upper = 0.0;
+ for (size_t k = 1; k < bl[0]->GetnMarkers(); ++k)
+ {
+ Cdf_upper += (frictioncoeff_upper[k] + frictioncoeff_upper[k - 1]) * (bl[0]->GetGlobMarkers(k) - bl[0]->GetGlobMarkers(k - 1)) / 2;
+ }
+ double Cdf_lower = 0.0;
+ for (size_t k = 1; k < bl[1]->GetnMarkers(); ++k)
+ {
+ Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (bl[1]->GetGlobMarkers(k) - bl[1]->GetGlobMarkers(k - 1)) / 2;
+ }
+
+ Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
- }
+ /* Compute total drag coefficient. */
- Cdf = Cdf_upper + Cdf_lower; // No friction in the wake
-
- /* Compute total drag coefficient. */
+ Cdt = 2 * bl[2]->U[lastWkPt + 0] * pow(bl[2]->U[lastWkPt + 3], (bl[2]->U[lastWkPt + 1] + 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. */
+ /* Compute pressure drag coefficient. */
- Cdp = Cdt - Cdf;
+ Cdp = Cdt - Cdf;
}
std::vector<std::vector<double>> ViscSolver::ExtractSolution()
{
- // [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();
+ // [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();
- for(size_t iRegion=0; iRegion<bl.size(); ++iRegion)
- {
- for (size_t iPoint=0; iPoint<bl[iRegion]->GetnMarkers(); ++iPoint)
+ for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion)
{
- 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
+ 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
+ }
}
- }
- return blVectors;
+ return blVectors;
}
std::vector<double> ViscSolver::Getxtr()
{
- std::vector<double> xtr(2);
- xtr[0] = bl[0]->xtr;
- xtr[1] = bl[1]->xtr;
- return xtr;
+ std::vector<double> xtr(2);
+ xtr[0] = bl[0]->xtr;
+ xtr[1] = bl[1]->xtr;
+ return xtr;
}
void ViscSolver::PrintBanner()
{
- std::cout << "\n\n" << std::endl;
- std::cout << "# ---------------------------------------------------------------------------- #" << std::endl;
- std::cout << "# #" << std::endl;
- std::cout << "# ///, //// #" << std::endl;
- std::cout << "# Paul Dechamps \\ /, / >. #" << std::endl;
- std::cout << "# ULiege 2022 \\ /, _/ /. #" << std::endl;
- std::cout << "# v0.1.0 - Memphis \\_ /_/ /. #" << std::endl;
- std::cout << "# \\__/_ < #" << std::endl;
- std::cout << "# /<<< \\_\\_ #" << std::endl;
- std::cout << "# ██╗ ██╗ ██████╗ ██████╗ ██╗ ██╗███████╗ /,)^>>_._ \\ #" << std::endl;
- std::cout << "# ██║ ██║██╔â•â•â•â–ˆâ–ˆâ•—██╔â•â•â–ˆâ–ˆâ•—██║ ██║██╔â•â•â•â•â• (/ \\ /\\\\ #" << std::endl;
- std::cout << "# ███████║██║ ██║██████╔â•â–ˆâ–ˆâ•‘ ██║███████╗ // ```` #" << std::endl;
- std::cout << "# ██╔â•â•â–ˆâ–ˆâ•‘██║ ██║██╔â•â•â–ˆâ–ˆâ•—██║ ██║╚â•â•â•â•â–ˆâ–ˆâ•‘ ======((`======= #" << std::endl;
- std::cout << "# ██║ ██║╚██████╔â•â–ˆâ–ˆâ•‘ ██║╚██████╔â•â–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ•‘ ___ ___ _ #" << std::endl;
- std::cout << "# â•šâ•â• â•šâ•â• â•šâ•â•â•â•â•â• â•šâ•â• â•šâ•â• â•šâ•â•â•â•â•â• â•šâ•â•â•â•â•â•â• / __| | _ ) | | #" << std::endl;
- std::cout << "# | (__ | _ \\ | |__ #" << std::endl;
- std::cout << "# High-Reynolds Operations for Robust Unsteady \\___| |___/ |____| #" << std::endl;
- std::cout << "# and Separating Coupled Boundary layers #" << std::endl;
- std::cout << "# #" << std::endl;
- std::cout << "# ---------------------------------------------------------------------------- #" << std::endl;
- std::cout << "\n\n" << std::endl;
-
+ std::cout << "\n\n"
+ << std::endl;
+ std::cout << "# ---------------------------------------------------------------------------- #" << std::endl;
+ std::cout << "# #" << std::endl;
+ std::cout << "# ///, //// #" << std::endl;
+ std::cout << "# Paul Dechamps \\ /, / >. #" << std::endl;
+ std::cout << "# ULiege 2022 \\ /, _/ /. #" << std::endl;
+ std::cout << "# v0.1.0 - Memphis \\_ /_/ /. #" << std::endl;
+ std::cout << "# \\__/_ < #" << std::endl;
+ std::cout << "# /<<< \\_\\_ #" << std::endl;
+ std::cout << "# ██╗ ██╗ ██████╗ ██████╗ ██╗ ██╗███████╗ /,)^>>_._ \\ #" << std::endl;
+ std::cout << "# ██║ ██║██╔â•â•â•â–ˆâ–ˆâ•—██╔â•â•â–ˆâ–ˆâ•—██║ ██║██╔â•â•â•â•â• (/ \\ /\\\\ #" << std::endl;
+ std::cout << "# ███████║██║ ██║██████╔â•â–ˆâ–ˆâ•‘ ██║███████╗ // ```` #" << std::endl;
+ std::cout << "# ██╔â•â•â–ˆâ–ˆâ•‘██║ ██║██╔â•â•â–ˆâ–ˆâ•—██║ ██║╚â•â•â•â•â–ˆâ–ˆâ•‘ ======((`======= #" << std::endl;
+ std::cout << "# ██║ ██║╚██████╔â•â–ˆâ–ˆâ•‘ ██║╚██████╔â•â–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ–ˆâ•‘ ___ ___ _ #" << std::endl;
+ std::cout << "# â•šâ•â• â•šâ•â• â•šâ•â•â•â•â•â• â•šâ•â• â•šâ•â• â•šâ•â•â•â•â•â• â•šâ•â•â•â•â•â•â• / __| | _ ) | | #" << std::endl;
+ std::cout << "# | (__ | _ \\ | |__ #" << std::endl;
+ std::cout << "# High-Reynolds Operations for Robust Unsteady \\___| |___/ |____| #" << std::endl;
+ std::cout << "# and Separating Coupled Boundary layers #" << std::endl;
+ std::cout << "# #" << std::endl;
+ std::cout << "# ---------------------------------------------------------------------------- #" << std::endl;
+ std::cout << "\n\n"
+ << std::endl;
}
\ No newline at end of file
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index e5c72884c7712fdd382eb35cb238aa2e5e4eeaf0..42a326ab048ecb15484ec22766dc710ad865b25c 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -10,49 +10,48 @@ namespace dart
class DART_API ViscSolver
{
public:
- double Cdt;
- double Cdf;
- double Cdp;
- double Re;
+ double Cdt,
+ Cdf,
+ Cdp;
+
+ bool printOn = true;
-private:
- double Minf;
- double CFL0;
- size_t nbElmUpper;
- bool stagPtMvmt;
- bool mapMesh=false;
-
- TimeSolver *tSolver;
- SolutionInterp *meshConverter;
-
-public:
- std::vector<BLRegion *> bl;
+ std::vector<std::vector<BLRegion *>> Sections;
- bool printOn = true;
+private:
+ double Re,
+ CFL0;
+ bool mapMesh = false;
- ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor);
+ std::vector<bool> stagPtMvmt;
+ std::vector<unsigned int> nbElmUpper; /* Nb of elements on the upper side. Used for stagnation point movement control. */
- ~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);
+ TimeSolver *tSolver;
+ SolutionInterp *meshConverter;
- 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();
+public:
+ ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor);
+ ~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();
private:
- void CheckWaves(size_t iPoint, BLRegion *bl);
- void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
- void ComputeDrag();
- void ComputeBlwVel();
- void PrintBanner();
-
+ void CheckWaves(size_t iPoint, BLRegion *bl);
+ void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
+ void ComputeDrag();
+ void ComputeBlwVel();
+ void PrintBanner();
};
} // namespace dart
#endif //WVISCSOLVER_H
\ No newline at end of file
diff --git a/dart/tests/bli_Polar.py b/dart/tests/bli_Polar.py
index 45330d4eb301b60cbafb132352e7a1570f1949a8..f5907b4567ed1d0728cecbe988382511c68ce379 100644
--- a/dart/tests/bli_Polar.py
+++ b/dart/tests/bli_Polar.py
@@ -58,10 +58,10 @@ def main():
# define flow variables
Re = 1e7
alpha=0
- alphaSeq = np.arange(-5,5.5,0.5)
+ alphaSeq = np.arange(-5,5.5,5)
M_inf = 0.
- meshFactor = 2
+ meshFactor = 1
CFL0 = 1
# ========================================================================================
@@ -105,7 +105,7 @@ def main():
print('SOLVERS statistics')
print(coupler.tms)
- viscU.PlotPolar(polar)
+ #viscU.PlotPolar(polar)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
diff --git a/dart/tests/bli_Sep.py b/dart/tests/bli_Sep.py
index 4477a8b355316d49c9ffc7abdadf88f90a7b3b67..85b9092ae60916604d5c7af7442d60b5399021ce 100644
--- a/dart/tests/bli_Sep.py
+++ b/dart/tests/bli_Sep.py
@@ -134,7 +134,7 @@ def main():
tests.add(CTest('Cd', vsolver.Cdt, 0.01452, 1e-3, forceabs=True))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0105, 1e-3, forceabs=True))
tests.add(CTest('xtr_top', xtr[0], 0.01, 0.03, forceabs=True))
- tests.add(CTest('xtr_bot', xtr[1], 0.99, 0.03, forceabs=True))
+ tests.add(CTest('xtr_bot', xtr[1], 1.00, 0.03, forceabs=True))
else:
raise Exception('Test not defined for this flow')
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 01084014b245d91a42d24910eeb3eaf51e6f48f5..b250e08a9a9cfb29130ce6fa5e36eb8c3cc0bdc7 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -37,6 +37,7 @@
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 numpy as np
@@ -62,7 +63,7 @@ def main():
M_inf = 0.
meshFactor = 1
- CFL0 = 5
+ CFL0 = 1
plotVar = 'H'
@@ -70,7 +71,12 @@ def main():
c_ref = 1
dim = 2
-
+
+ if dim == 2:
+ nSection = 1
+ else:
+ nSection = 10
+
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
@@ -88,10 +94,9 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- isolver = floD.newton(pbl)
vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
-
- coupler = floVC.Coupler(isolver, vsolver)
+ isolverAPI = viscAPI.dartAPI(floD.newton(pbl), blw[0], blw[1])
+ coupler = floVC.Coupler(isolverAPI, vsolver)
coupler.CreateProblem(blw[0], blw[1])
diff --git a/dart/tests/bli_python.py b/dart/tests/blipython.py
similarity index 97%
rename from dart/tests/bli_python.py
rename to dart/tests/blipython.py
index ca729b16b449ff9e27ee27bac861781c2e107435..99e0816e43936b69b96f65126798551e9a59836e 100755
--- a/dart/tests/bli_python.py
+++ b/dart/tests/blipython.py
@@ -60,8 +60,8 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5.*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]
@@ -86,9 +86,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.001}
+ 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
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 9e721714c9adad9b9173aa138c9d9e1f67912a8c..bb82b95736333ad634de34e6eef44a091f1da1c7 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -29,12 +29,11 @@ import math
import fwk
import dart.viscUtils as viscU
import numpy as np
-from matplotlib import pyplot as plt
class Coupler:
- def __init__(self, iSolver, vSolver) -> None:
- self.iSolver=iSolver
- self.vSolver=vSolver
+ def __init__(self, _iSolverAPI, vSolver) -> None:
+ self.iSolverAPI = _iSolverAPI
+ self.vSolver = vSolver
self.maxCouplIter = 150
self.couplTol = 1e-4
@@ -43,14 +42,6 @@ class Coupler:
self.resetInv = False
pass
- def CreateProblem(self, _boundaryAirfoil, _boundaryWake):
-
- self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
- self.blwVel = [None, None, None]
- self.deltaStar = [None, None, None]
- self.xx = [None, None, None]
- pass
-
def Run(self):
# Initilize meshes in c++ objects
@@ -67,8 +58,8 @@ class Coupler:
# Run inviscid solver
self.tms['inviscid'].start()
if self.resetInv:
- self.iSolver.reset()
- self.iSolver.run()
+ self.iSolverAPI.ResetSolver()
+ self.iSolverAPI.RunSolver()
self.tms['inviscid'].stop()
# Extract inviscid boundary
diff --git a/dart/viiIPS.py b/dart/viiIPS.py
new file mode 100644
index 0000000000000000000000000000000000000000..61f5d957fe6f1011ae109f723ea91e238580a2d4
--- /dev/null
+++ b/dart/viiIPS.py
@@ -0,0 +1,191 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Information passing structures (IPS) for communication between solvers.
+
+import numpy as np
+
+class viscousIPS:
+ def __init__(self, _name) -> None:
+ self.name = _name
+
+ 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]))
+
+
+
+class inviscidIPS:
+ def __init__(self, _boundary):
+ self.boundary = _boundary # boundary of interest
+ self.nN = len(self.boundary.nodes) # number of nodes
+ self.nE = len(self.boundary.tag.elems) # number of elements
+ self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
+ self.u = np.zeros(self.nE) # blowing velocity
+ self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
+ self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
+ 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
+
+class Airfoil(inviscidIPS):
+ def __init__(self, _boundary):
+ inviscidIPS.__init__(self, _boundary)
+
+ def connectList(self):
+ ''' Sort the value read by the viscous solver/ Create list of connectivity
+ '''
+ N1 = np.zeros(self.nN, dtype=int) # Node number
+ connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+ connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+ data = np.zeros((self.boundary.nodes.size(), 10))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos[0]
+ data[i,2] = n.pos[1]
+ data[i,3] = n.pos[2]
+ data[i,4] = self.v[i,0]
+ data[i,5] = self.v[i,1]
+ data[i,6] = self.Me[i]
+ data[i,7] = self.rhoe[i]
+ data[i,8] = self.deltaStar[i]
+ data[i,9] = self.xx[i]
+ i += 1
+ # Table containing the element and its nodes
+ eData = np.zeros((self.nE,3), dtype=int)
+ for i in range(0, self.nE):
+ eData[i,0] = self.boundary.tag.elems[i].no
+ 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
+ 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
+ idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
+ idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
+ ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
+ ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
+ if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
+ upperGlobTE = int(data[idxTE[0],0])
+ lowerGlobTE = int(data[idxTE[1],0])
+ else:
+ upperGlobTE = int(data[idxTE[1],0])
+ lowerGlobTE = int(data[idxTE[0],0])
+ # Connectivity
+ connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+ N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
+ connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+ i = 1
+ upperTE = 0
+ lowerTE = 0
+ # Sort the suction part
+ while upperTE == 0:
+ N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+ connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
+ connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+ if eData[connectListElems[i],2] == upperGlobTE:
+ upperTE = 1
+ i += 1
+ # Sort the pressure side
+ connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+ connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+ N1[i] = eData[connectListElems[i],2]
+ while lowerTE == 0:
+ N1[i+1] = eData[connectListElems[i],1] # First node of the element
+ connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
+ connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+ if eData[connectListElems[i+1],1] == lowerGlobTE:
+ lowerTE = 1
+ i += 1
+ connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
+ data[:,0] = data[connectListNodes,0]
+ data[:,1] = data[connectListNodes,1]
+ data[:,2] = data[connectListNodes,2]
+ data[:,3] = data[connectListNodes,3]
+ data[:,4] = data[connectListNodes,4]
+ data[:,5] = data[connectListNodes,5]
+ data[:,6] = data[connectListNodes,6]
+ data[:,7] = data[connectListNodes,7]
+ data[:,8] = data[connectListNodes,8]
+ data[:,9] = data[connectListNodes,9]
+ # Separated upper and lower part
+ data = np.delete(data,0,1)
+ uData = data[0:np.argmax(data[:,0])+1]
+ lData = data[np.argmax(data[:,0])+1:None]
+ lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+ return connectListNodes, connectListElems,[uData, lData]
+
+class Wake(inviscidIPS):
+ def __init__(self, _boundary):
+ inviscidIPS.__init__(self, _boundary)
+
+ def connectList(self):
+ ''' Sort the value read by the viscous solver/ Create list of connectivity
+ '''
+ N1 = np.zeros(self.nN, dtype=int) # Node number
+ connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+ connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+ data = np.zeros((self.boundary.nodes.size(), 10))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos[0]
+ data[i,2] = n.pos[1]
+ data[i,3] = n.pos[2]
+ data[i,4] = self.v[i,0]
+ data[i,5] = self.v[i,1]
+ data[i,6] = self.Me[i]
+ data[i,7] = self.rhoe[i]
+ data[i,8] = self.deltaStar[i]
+ data[i,9] = self.xx[i]
+ i += 1
+ # Table containing the element and its nodes
+ eData = np.zeros((self.nE,3), dtype=int)
+ for i in range(0, self.nE):
+ eData[i,0] = self.boundary.tag.elems[i].no
+ eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+ eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+ connectListNodes = data[:,1].argsort()
+ connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
+ for i in range(1, len(eData[:,0])):
+ connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
+ data[:,0] = data[connectListNodes,0]
+ data[:,1] = data[connectListNodes,1]
+ data[:,2] = data[connectListNodes,2]
+ data[:,3] = data[connectListNodes,3]
+ data[:,4] = data[connectListNodes,4]
+ data[:,5] = data[connectListNodes,5]
+ data[:,6] = data[connectListNodes,6]
+ data[:,7] = data[connectListNodes,7]
+ data[:,8] = data[connectListNodes,8]
+ data[:,9] = data[connectListNodes,9]
+ # Separated upper and lower part
+ data = np.delete(data,0,1)
+ return connectListNodes, connectListElems, data
\ No newline at end of file
diff --git a/dart/viscAPI.py b/dart/viscAPI.py
new file mode 100644
index 0000000000000000000000000000000000000000..12a99ea795d4a2c186dbcb95f21ee8c1de84a41b
--- /dev/null
+++ b/dart/viscAPI.py
@@ -0,0 +1,66 @@
+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
diff --git a/dart/viscUtils.py b/dart/viscUtils.py
index 763036be821231b502f18cf9f64cdbfc25608207..74b523eff13fd159fbcb8f0559c63af6861d7552 100644
--- a/dart/viscUtils.py
+++ b/dart/viscUtils.py
@@ -1,4 +1,8 @@
+import dart
+def initBL(Re, Minf, CFL0, meshFactor):
+ solver = dart.ViscSolver(Re, Minf, CFL0, meshFactor)
+ return solver
def Dictionary(blScalars, blVectors, xtr):