diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index 54d3a31e17266bd06bf7acd18aac25e805049383..5d2bae91d01c90df1fd03e8233753424a8794503 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -26,3 +26,7 @@
#include "wProblem.h"
#include "wSolver.h"
#include "wWake.h"
+#include "wViscSolver.h"
+#include "wBLRegion.h"
+#include "wTimeSolver.h"
+#include "wViscFluxes.h"
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index e16e83975d1e51bdb49efa367eb14f935deb0e43..b1eba25b6a021bf5d6c1d3a84b4ea45112ff5420 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -58,6 +58,11 @@ threads="1"
%shared_ptr(dart::Newton);
%shared_ptr(dart::Adjoint);
+%include "wViscSolver.h"
+%include "wBLRegion.h"
+%include "wTimeSolver.h"
+%include "wViscFluxes.h"
+
%include "wFluid.h"
%include "wAssign.h"
%include "wFreestream.h"
diff --git a/dart/default.py b/dart/default.py
index 3a03df5626170c58ef3be646579d518954072a69..e12c4ab8c0118f8dfcf4b2061cc41f840f94da4b 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,6 +23,12 @@ from fwk.wutils import parseargs
import dart
import tbox
+def initBL(Re, Minf, CFL0):
+ print('Re default.py', Re, Minf, CFL0)
+ solver = dart.ViscSolver(Re, Minf, CFL0)
+ print('ok')
+ return solver
+
def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
"""Initialize mesh and mesh writer
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 2afc45904a47f3012a57caf58b553233bfe721cc..72c71064973cf3954e5cdb8e85577529417d6c1b 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -64,6 +64,10 @@ class Picard;
class Newton;
class FpeLSFunction;
class Adjoint;
+class ViscSolver;
+class BLRegion;
+class TimeSolver;
+class ViscSolver;
// General
class F0Ps;
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..56115e9bc067e56ae18a2a153a11ce26f0c33449
--- /dev/null
+++ b/dart/src/wBLRegion.cpp
@@ -0,0 +1,168 @@
+#include "wBLRegion.h"
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+using namespace tbox;
+using namespace dart;
+
+
+BLRegion::BLRegion()
+{
+ std::cout << "Coucou de BLRegion\n" << std::endl;
+
+
+} // end BLRegion
+
+
+BLRegion::~BLRegion()
+{
+ 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);
+} // End SetMesh
+
+void BLRegion::SetInvBC(std::vector<double> _Ue,
+ std::vector<double> _Me,
+ std::vector<double> _Rhoe)
+{
+ VtExt = _Ue;
+ Me = _Me;
+ Rhoe = _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] = 0;
+ U[3] = Ret[0]/(U[0]*Re);
+ } // End if
+
+ LaminarClosures(0);
+} // End SetStagBC
+
+unsigned int BLRegion::GetnVar(){return nVar;}
+size_t BLRegion::GetnMarkers(){return nMarkers;}
+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 << "Point " << iPoint << std::endl;
+ std::cout << "Theta " << U[iPoint * nVar + 0] << " \n" << std::endl;
+ std::cout << "H " << U[iPoint * nVar + 1] << " \n" << std::endl;
+ std::cout << "N " << U[iPoint * nVar + 2] << " \n" << std::endl;
+ std::cout << "ue " << U[iPoint * nVar + 3] << " \n" << std::endl;
+ std::cout << "Ct " << U[iPoint * nVar + 4] << " \n" << std::endl;
+ }
+
+void BLRegion::LaminarClosures(size_t iPoint)
+{
+ double theta = U[iPoint*nVar+0];
+ double H = U[iPoint*nVar+1];
+ double ret = Ret[iPoint];
+ double me = Me[iPoint];
+ double hk = (H - 0.29*me*me)/(1+0.113*me*me); // Kinematic shape factor
+ if(name == "airfoil"){
+ hk = std::max(hk, 1.02);
+ }
+ else{
+ hk = std::max(hk, 1.00005);
+ }
+ hk = std::min(hk,6.0);
+ double hstar2 = (0.064/(hk-0.8)+0.251)*me*me; // Density shape parameter
+ double delta = std::min((3.15+ H + (1.72/(hk-1)))*theta, 12*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);
+ }
+ 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
+
+ /* Cf */
+
+ double cf;
+ if(hk < 5.5){
+ double tmp = (5.5-hk)*(5.5-hk)*(5.5-hk) / (hk+1);
+ cf = (-0.07 + 0.0727*tmp)/ret;
+ }
+ else if(hk >= 5.5){
+ double tmp = 1 - 1/(hk-4.5);
+ cf = (-0.07 + 0.015*tmp*tmp) /ret;
+ }
+
+ /* Cd */
+
+ double cd;
+ 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;
+ }
+ Cd[iPoint] = cd;
+ Cf[iPoint] = cf;
+ Delta[iPoint] = delta;
+ Hstar[iPoint] = hstar;
+ Hstar2[iPoint] = hstar2;
+ Hk[iPoint] = hk;
+ Cteq[iPoint] = 0;
+ Us[iPoint] = 0;
+}
+
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f0e2ac2e4a8ad6418bbb1f958aba5ebd80bd94c
--- /dev/null
+++ b/dart/src/wBLRegion.h
@@ -0,0 +1,88 @@
+#ifndef WBLREGION_H
+#define WBLREGION_H
+
+#include "dart.h"
+#include <vector>
+#include <string>
+
+namespace dart
+{
+
+class DART_API BLRegion
+{
+
+private:
+
+ unsigned int nVar = 5;
+
+ /* Mesh parameters */
+
+ 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;
+
+public:
+
+ std::string name;
+
+ /* 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<double> Regime;
+
+ 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 SetStagBC(double Re);
+
+ unsigned int GetnVar();
+ size_t GetnMarkers();
+ 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 LaminarClosures(size_t iPoint);
+};
+} // namespace dart
+#endif //WBLREGION_H
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..47c09d27f7fcddd6a5f894669116d849db06217b
--- /dev/null
+++ b/dart/src/wTimeSolver.cpp
@@ -0,0 +1,150 @@
+#include "wTimeSolver.h"
+#include "wBLRegion.h"
+#include "wViscFluxes.h"
+#include <Eigen/Dense>
+#include <iostream>
+#include <vector>
+#include <cmath>
+
+using namespace Eigen;
+using namespace tbox;
+using namespace dart;
+
+
+TimeSolver::TimeSolver(double _CFL0, double _Minf, unsigned long _maxIt, double _tol, unsigned int _itFrozenJac)
+{
+ std::cout << "Coucou de TimeSolver\n" << std::endl;
+
+ CFL0 = _CFL0;
+ maxIt = _maxIt;
+ tol = _tol;
+ itFrozenJac0 = _itFrozenJac;
+ Minf = _Minf;
+
+ SysMatrix = new ViscFluxes();
+
+} // end TimeSolver
+
+
+TimeSolver::~TimeSolver()
+{
+ 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
+} // End InitialCondition
+
+int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
+{
+
+ // Save initial condition
+
+ unsigned int nVar = bl->GetnVar();
+
+ std::vector<double> temp_U;
+ for (size_t i=0; i<nVar; ++i){
+ temp_U.push_back(bl->U[iPoint*nVar+i]);
+ } // End for
+
+ // Initialize time integration variables
+
+ double dx = bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1);
+
+ 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;
+
+ Vector<double, 5> SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
+ double normSysRes0 = SysRes.norm();
+ double normSysRes;
+
+ Matrix<double, 5, 5> JacMatrix;
+ Vector<double, 5> slnIncr;
+
+ unsigned int nErrors = 0; // Number of errors encountered
+ unsigned int innerIters = 0; // Inner (non-linear) iterations
+ while (innerIters < maxIt){
+
+ if(innerIters % itFrozenJac){
+ JacMatrix = SysMatrix->ComputeJacobian(iPoint, bl, SysRes, 1e-8);
+ }
+
+ // Ici on peut faire autrement : .asDiagonal()
+ slnIncr = (JacMatrix + IMatrix).partialPivLu().solve(SysRes);
+
+ for(size_t k=0; k<slnIncr.size(); ++k){
+ bl->U[iPoint*nVar+k] += slnIncr(k);
+ }
+
+ SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
+ normSysRes = (SysRes + slnIncr/timeStep).norm();
+
+ // CFL adaptation.
+
+ CFL = CFL0* pow(normSysRes0/normSysRes, 0.7);
+
+ CFL = std::max(CFL, 1.0);
+ timeStep = SetTimeStep(CFL, dx, bl->GetVtExt(iPoint));
+ IMatrix = Matrix<double, 5, 5>::Identity() / timeStep;
+
+ innerIters++;
+ } // End while
+ return innerIters;
+ return 0;
+} // End Integration
+
+double TimeSolver::SetTimeStep(double CFL, double dx, double invVel)
+{
+ /* Speed of sound. */
+ double gradU2 = (invVel*invVel);
+ double soundSpeed = ComputeSoundSpeed(gradU2);
+
+ /* Velocity of the fastest travelling wave */
+ double advVel = soundSpeed + invVel;
+
+ /* Time step computation */
+ return CFL*dx/advVel;
+}
+
+double TimeSolver::ComputeSoundSpeed(double gradU2)
+{
+ return sqrt(1 / (Minf * Minf) + 0.5 * (gamma - 1) * (1 - gradU2));
+}
+
+
+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
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..b85ac58145e51425189fc0ca6f421c09f2f50cbf
--- /dev/null
+++ b/dart/src/wTimeSolver.h
@@ -0,0 +1,45 @@
+#ifndef WTIMESOLVER_H
+#define WTIMESOLVER_H
+
+#include "dart.h"
+#include "wViscFluxes.h"
+#include <vector>
+#include <string>
+
+namespace dart
+{
+
+class DART_API TimeSolver
+{
+
+private:
+
+ double CFL0;
+ unsigned long maxIt;
+ double tol;
+ unsigned int itFrozenJac0;
+ bool initSln = true;
+ double gamma = 1.4;
+ double Minf;
+
+ ViscFluxes *SysMatrix;
+
+public:
+ TimeSolver(double _CFL0, double _Minf, unsigned long _maxIt = 1000, 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);
+};
+} // namespace dart
+#endif //WTimeSolver_H
\ No newline at end of file
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0374f60d19eea90895555a9e856291bb2f69c4c6
--- /dev/null
+++ b/dart/src/wViscFluxes.cpp
@@ -0,0 +1,399 @@
+#define _USE_MATH_DEFINES
+
+#include "wViscFluxes.h"
+#include "wBLRegion.h"
+#include <iostream>
+#include <Eigen/Dense>
+#include <cmath>
+#include <algorithm>
+#include <string>
+
+using namespace dart;
+using namespace Eigen;
+
+
+ViscFluxes::ViscFluxes()
+{
+ std::cout << "Coucou de ViscFluxes\n" << std::endl;
+
+
+} // end ViscFluxes
+
+
+ViscFluxes::~ViscFluxes()
+{
+ 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);
+}
+
+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
+
+ for (size_t iVar=0; iVar<nVar; ++iVar){
+ double 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"){
+ double dissipRatio = 0.9;
+ } // End if
+ else{
+ double dissiRatio = 1;
+ } // End else
+
+ bl->Ret[iPoint] = std::max(u[3] * u[0], 1.0);
+ bl->DeltaStar[iPoint] = u[0] * u[1];
+
+ /* Boundary layer closure relations */
+
+ if (bl->Regime[iPoint] == 0){
+ LaminarClosures(iPoint, bl);
+ } // End if
+ else if (bl->Regime[iPoint] == 1){
+ TurbulentClosures(iPoint, bl);
+ } // 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 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))/dx;
+ double ddeltaStarExt_dx = (bl->GetDeltaStarExt(iPoint) - bl->GetDeltaStarExt(iPoint))/dx;
+
+ double c=4/(M_PI*dx);
+ double cExt=4/(M_PI*(bl->GetXxExt(iPoint) - bl->GetXxExt(iPoint-1)));
+
+ /* Amplification routine */
+ double dN_dx;
+ double Ax;
+
+ if (bl->xtrF!=NULL && bl->Regime[iPoint] == 0){
+ double Ax = AmplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
+ double dN_dx = (bl->U[iPoint * nVar + 2] - bl->U[(iPoint-1)*nVar + 2])/(bl->GetLocMarkers(iPoint) - bl->GetLocMarkers(iPoint-1));
+
+ } // End if
+ else{
+ double Ax = 0;
+ double dN_dx = 0;
+ } // End else
+
+ double Mea = bl->GetMe(iPoint);
+
+ double Reta = bl->Ret[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 Usa = bl->Us[iPoint];
+
+ /* Space part */
+
+ spaceVector(0) = dT_dx + (2 + u[1] - Mea * Mea) * (u[0]/u[3]) * due_dx - Cfa;
+ spaceVector(1) = u[1] * 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*(bl->Cteq[iPoint] - 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;
+ }
+
+ /* Time part */
+
+ timeMatrix(0,0) = u[1]/u[3];
+ timeMatrix(0,1) = u[0]/u[3];
+ timeMatrix(0,2) = 0;
+ timeMatrix(0,3) = u[0]*u[1]/(u[3]*u[3]);
+ timeMatrix(0,4) = 0;
+
+ timeMatrix(1,0) = (1+u[1]*(1-Hstara))/u[3];
+ timeMatrix(1,1) = (1-Hstara)*u[1]/u[3];
+ timeMatrix(1,2) = 0;
+ timeMatrix(1,3) = (2-Hstara*u[1])*u[0]/(u[3]*u[3]);
+ timeMatrix(1,4) = 0;
+
+ timeMatrix(2,0) = 0;
+ timeMatrix(2,1) = 0;
+ timeMatrix(2,2) = 1;
+ timeMatrix(2,3) = 0;
+ timeMatrix(2,4) = 0;
+
+ timeMatrix(3,0) = -c*u[1];
+ timeMatrix(3,1) = -c*u[0];
+ timeMatrix(3,2) = 0;
+ timeMatrix(3,3) = 1;
+ timeMatrix(3,4) = 0;
+
+ timeMatrix(4,0) = 0;
+ timeMatrix(4,1) = 0;
+ timeMatrix(4,2) = 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;
+ timeMatrix(4,4) = 0;
+ }
+
+ Vector<double, 5> sysRes = timeMatrix.partialPivLu().solve(spaceVector);
+ return sysRes;
+}
+
+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;
+}
+
+void ViscFluxes::LaminarClosures(size_t iPoint, BLRegion *bl)
+{
+ unsigned int nVar = bl->GetnVar();
+ double theta = bl->U[iPoint*nVar+0];
+ double H = bl->U[iPoint*nVar+1];
+ double Ret = bl->Ret[iPoint];
+ double Me = bl->GetMe(iPoint);
+ double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me); // Kinematic shape factor
+ if(bl->name == "airfoil"){
+ Hk = std::max(Hk, 1.02);
+ }
+ else{
+ Hk = std::max(Hk, 1.00005);
+ }
+ Hk = std::min(Hk,6.0);
+ double Hstar2 = (0.064/(Hk-0.8)+0.251)*Me*Me; // Density shape parameter
+ double delta = std::min((3.15+ H + (1.72/(Hk-1)))*theta, 12*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);
+ }
+ 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
+
+ /* Cf */
+
+ double Cf;
+ if(Hk < 5.5){
+ double tmp = (5.5-Hk)*(5.5-Hk)*(5.5-Hk) / (Hk+1);
+ Cf = (-0.07 + 0.0727*tmp)/Ret;
+ }
+ else if(Hk >= 5.5){
+ double tmp = 1 - 1/(Hk-4.5);
+ Cf = (-0.07 + 0.015*tmp*tmp) /Ret;
+ }
+
+ /* Cd */
+
+ double Cd;
+ 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 (bl->name == "wake"){
+ Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
+ Cf = 0;
+ }
+ bl->Cd[iPoint] = Cd;
+ bl->Cf[iPoint] = Cf;
+ bl->Delta[iPoint] = delta;
+ bl->Hstar[iPoint] = Hstar;
+ bl->Hstar2[iPoint] = Hstar2;
+ bl->Hk[iPoint] = Hk;
+ bl->Cteq[iPoint] = 0;
+ bl->Us[iPoint] = 0;
+}
+
+void ViscFluxes::TurbulentClosures(size_t iPoint, BLRegion *bl)
+{
+ unsigned int nVar = bl->GetnVar();
+ double theta = bl->U[iPoint*nVar+0];
+ double H = bl->U[iPoint*nVar+1];
+ double Ct = bl->U[iPoint*nVar+4];
+ double Ret = bl->Ret[iPoint];
+ double Me = bl->GetMe(iPoint);
+
+ H = std::max(H, 1.0005);
+ /* Eliminate absurd transients */
+ Ct = std::max(std::min(Ct, 0.30), 0.0000001);
+ double Hk = (H - 0.29*Me*Me)/(1+0.113*Me*Me);
+ Hk = std::max(std::min(Hk, 10.0), 1.00005);
+ double Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me*Me;
+ //gam = 1.4 - 1
+ //Fc = np.sqrt(1+0.5*gam*Me**2)
+ double Fc = sqrt(1+0.2*Me*Me);
+
+ double H0;
+ if(Ret <= 400){
+ H0 = 4;
+ }
+ else if(Ret > 400){
+ H0 = 3 + (400/Ret);
+ }
+
+ Ret = std::max(Ret, 200.0);
+ 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 if (Hk > H0){
+ 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
+ //logRt = np.log(Ret/Fc)
+ //logRt = max(logRt,3)
+ //arg = max(-1.33*Hk, -20)
+ 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);
+
+ int n;
+ double Cf;
+ double Hkc;
+ double Cdw;
+ double Cdd;
+ double Cdl;
+ double Cteq;
+ double Hmin;
+ double Fl;
+ double Dfac;
+
+ if (bl->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;
+ } // End if
+ n = 1;
+ // Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-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
+ Hmin = 1 + 2.1/log(Ret);
+ Fl = (Hk-1)/(Hmin-1);
+ 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
+ } // End else
+ double Cd = n*(Cdw+ Cdd + Cdl);
+ // Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk)*((Hk-1)/(6.432*Hk)));
+ bl->Cd[iPoint] = Cd;
+ bl->Cf[iPoint] = Cf;
+ bl->Delta[iPoint] = delta;
+ bl->Hstar[iPoint] = Hstar;
+ bl->Hstar2[iPoint] = Hstar2;
+ bl->Hk[iPoint] = Hk;
+ bl->Cteq[iPoint] = Cteq;
+ bl->Us[iPoint] = Us;
+}
\ No newline at end of file
diff --git a/dart/src/wViscFluxes.h b/dart/src/wViscFluxes.h
new file mode 100644
index 0000000000000000000000000000000000000000..8e9ac8ef9fe59096aa9b93e8c3cd5d4521f83b9a
--- /dev/null
+++ b/dart/src/wViscFluxes.h
@@ -0,0 +1,28 @@
+#ifndef WViscFluxes_H
+#define WViscFluxes_H
+
+#include "dart.h"
+#include "wBLRegion.h"
+#include <Eigen/Dense>
+#include <vector>
+#include <string>
+
+namespace dart
+{
+
+class DART_API ViscFluxes
+{
+public:
+ ViscFluxes();
+ ~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);
+ void LaminarClosures(size_t iPoint, BLRegion *bl);
+ void TurbulentClosures(size_t iPoint, BLRegion *bl);
+
+};
+} // namespace dart
+#endif //WViscFluxes_H
\ No newline at end of file
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ba78fd28f634469afb37ffa462b03d446e81151
--- /dev/null
+++ b/dart/src/wViscSolver.cpp
@@ -0,0 +1,110 @@
+#include "wViscSolver.h"
+#include "wBLRegion.h"
+#include "wTimeSolver.h"
+#include <iostream>
+#include <vector>
+
+
+using namespace tbox;
+using namespace dart;
+
+
+ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
+{
+
+ /* Freestream parameters */
+ Re = _Re;
+ Minf = _Minf;
+
+ /* User input numerical parameters */
+ CFL0 = _CFL0;
+
+ /* Parameters */
+ double cd = 0.0;
+
+ std::cout << "--- Viscous solver setup ---\n"
+ << "Reynolds number: " << Re
+ << "Freestream Mach number: " << Minf
+ << "Initial CFL number: " << CFL0
+ << std::endl;
+
+ /* Initialzing boundary layer */
+
+ //std::vector<dart::BLRegion> bl;
+
+ std::cout << "--- Jusqu'ici ok ---\n" << std::endl;
+
+ for (size_t i=0; i<3; ++i)
+ {
+ std::cout << "yop" << std::endl;
+ BLRegion *region = new BLRegion();
+ std::cout << "yop2" << std::endl;
+ bl.push_back(region);
+ } // end for
+
+ tSolver = new TimeSolver(CFL0, Minf);
+
+ bl[0]->name = "upper";
+ bl[1]->name = "lower";
+ bl[2]->name = "wake";
+
+ std::cout << " Boundary layer initialized.\n" << std::endl;
+
+} // end ViscSolver
+
+
+ViscSolver::~ViscSolver()
+{
+ for (size_t i=0; i<bl.size(); ++i){
+ delete bl[i];
+ } // end for
+
+ delete tSolver;
+
+ std::cout << "~ViscSolver()\n";
+
+
+} // end ~ViscSolver
+
+void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM)
+{
+ bl[region]->SetMesh(locM, globM);
+}
+
+void ViscSolver::SetInvBC(size_t region, std::vector<double> _Ue,
+ std::vector<double> _Me,
+ std::vector<double> _Rhoe)
+{
+ bl[region]->SetInvBC(_Ue, _Me, _Rhoe);
+}
+
+
+int ViscSolver::Run(){
+ std::cout << "Entered run"
+ << std::endl;
+
+ /* Set boundary condition */
+
+ bl[0]->SetStagBC(Re);
+ bl[1]->SetStagBC(Re);
+
+ /* Loop over the regions */
+
+ for (size_t iRegion = 0; iRegion < bl.size(); ++iRegion){
+
+ /* Loop over points */
+
+ for (size_t iPoint = 1; iPoint < bl[iRegion]->GetnMarkers(); ++iPoint){
+
+ tSolver->InitialCondition(iPoint, bl[iRegion]);
+
+ std::cout << "Point " << iPoint << " H " << bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 3] << std::endl;
+
+ //tSolver->Integration(iPoint, bl[iRegion]);
+
+ }
+
+
+ }
+ return 0; // Successfull exit
+};
\ No newline at end of file
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..5a369f480d4b07620b19963dbe2fd9aeb6bc8290
--- /dev/null
+++ b/dart/src/wViscSolver.h
@@ -0,0 +1,38 @@
+#ifndef WVISCSOLVER_H
+#define WVISCSOLVER_H
+
+#include "dart.h"
+#include "wBLRegion.h"
+
+namespace dart
+{
+
+class DART_API ViscSolver
+{
+public:
+ double cd;
+
+private:
+ double Re;
+ double Minf;
+ double CFL0;
+
+ TimeSolver *tSolver;
+
+public:
+ std::vector<BLRegion *> bl;
+
+ ViscSolver(double _Re, double _Minf, double _CFL0);
+
+ ~ViscSolver();
+ void InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM);
+ void SetInvBC(size_t region, std::vector<double> _Ue,
+ std::vector<double> _Me,
+ std::vector<double> _Rhoe);
+
+
+ int Run();
+
+};
+} // namespace dart
+#endif //WVISCSOLVER_H
\ No newline at end of file
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
index c29ee906ffb6d1829ea00f41e0ff26c0e5e64264..f65132bf5140085e838e3bd3c6c985e7cb43c8d5 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -62,9 +62,9 @@ def main():
# define flow variables
Re = 1e7
alpha = 2*math.pi/180
- M_inf = 0.73
+ M_inf = 0.715
- user_xtr=[0, 0]
+ user_xtr=[0.43, None]
user_Ti=None
mapMesh = 1
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
new file mode 100644
index 0000000000000000000000000000000000000000..8391d1f52c94c7d5bca46afb408361614aaf1d98
--- /dev/null
+++ b/dart/tests/blicpp.py
@@ -0,0 +1,152 @@
+#!/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.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+# -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+
+import dart.vCoupler as floVC
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def main():
+ print('Start')
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ Re = 1e7
+ alpha = 0.*math.pi/180
+ M_inf = 0.
+
+ #user_xtr=[0.41,0.41]
+ #user_xtr=[0,None]
+ user_xtr=[None,None]
+ user_Ti=None
+
+ mapMesh = 0
+ nFactor = 3
+
+ # Time solver parameters
+ Time_Domain = True # True for unsteady solver, False for steady solver
+ CFL0 = 5
+
+ # ========================================================================================
+
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ c_ref = 1
+ dim = 2
+ tol = 1e-4 # tolerance for the VII (usually one drag count)
+ case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ tms['msh'].stop()
+
+ # set the problem and add medium, initial/boundary conditions
+ tms['pre'].start()
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+ tms['pre'].stop()
+
+ # solve viscous problem
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+
+ isolver = floD.newton(pbl)
+ vsolver = floD.initBL(Re, M_inf, CFL0)
+
+ coupler = floVC.Coupler(isolver, vsolver)
+
+ coupler.CreateProblem(blw[0], blw[1])
+
+ coupler.Run()
+ tms['solver'].stop()
+
+ # extract Cp
+ Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+
+
+
+ # visualize solution and plot results
+ floD.initViewer(pbl)
+ tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
+
+ val=validation.Validation(case)
+ val.main(isolver.Cl,vsolver.wData)
+
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
+ tests.add(CTest('Cd', vsolver.Cd, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
+ else:
+ raise Exception('Test not defined for this flow')
+
+ tests.run()
+
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a575f0732bb7db28b8293c55acd9b8d1406e996
--- /dev/null
+++ b/dart/vCoupler.py
@@ -0,0 +1,114 @@
+#!/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.
+
+
+## VII Coupler
+# Paul Dechamps
+
+from dart.viscous.Master.DAirfoil import Airfoil
+from dart.viscous.Master.DWake import Wake
+from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
+
+from fwk.coloring import ccolors
+
+import numpy as np
+
+class Coupler:
+ def __init__(self, iSolver, vSolver) -> None:
+ self.iSolver=iSolver
+ self.vSolver=vSolver
+
+ self.maxCouplIter = 2
+ self.couplTol = 1e-4
+ pass
+
+ def CreateProblem(self, _boundaryAirfoil, _boundaryWake):
+
+ self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+ pass
+
+ def Run(self):
+
+ # Initilize meshes in c++ objects
+
+ couplIter = 0
+ while couplIter < self.maxCouplIter:
+
+ cdPrev = 0
+
+ # Run inviscid solver
+
+ self.iSolver.run()
+
+ # Extract inviscid boundary
+
+ for n in range(0, len(self.group)):
+
+ for i in range(0, len(self.group[n].boundary.nodes)):
+
+ self.group[n].v[i,0] = self.iSolver.U[self.group[n].boundary.nodes[i].row][0]
+ self.group[n].v[i,1] = self.iSolver.U[self.group[n].boundary.nodes[i].row][1]
+ self.group[n].v[i,2] = self.iSolver.U[self.group[n].boundary.nodes[i].row][2]
+ self.group[n].Me[i] = self.iSolver.M[self.group[n].boundary.nodes[i].row]
+ self.group[n].rhoe[i] = self.iSolver.rho[self.group[n].boundary.nodes[i].row]
+
+ dataIn=[None,None]
+ for n in range(0,len(self.group)):
+ self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList()
+ fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, 0, 1)
+ if couplIter == 0:
+ # Initialize mesh
+ self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
+ self.vSolver.InitMeshes(1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ self.vSolver.InitMeshes(2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:], fMeshDict['globCoord'][fMeshDict['bNodes'][2]:])
+ print('fro mpy : Mesh ok')
+ self.vSolver.SetInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
+ self.vSolver.SetInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ self.vSolver.SetInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
+ print('from py: InvBC ok')
+ # Run viscous solver
+
+ exitCode = self.vSolver.Run()
+ print('From python exitcode: ', exitCode)
+
+ error = (self.vSolver.cd - cdPrev)/self.vSolver.cd if self.vSolver.cd != 0 else 1
+
+ if error <= self.couplTol:
+ print('')
+ print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
+ .format(couplIter, self.vsolver.cd, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
+ print('')
+ return 0
+
+ print('')
+ print(ccolors.ANSI_BLUE, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} > {:<2.3f}'
+ .format(couplIter, self.vSolver.cd, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
+ print('')
+
+ # Set inviscid boundary condition
+
+ for n in range(0, len(self.group)):
+
+ for i in range(0, self.group[n].nE):
+
+ self.group[n].boundary.setU(i, self.group[n].u[i])
+
+ couplIter += 1
+ else:
+ print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
+ return couplIter
+
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index 44877c7d62e24180c383d8319156baa755c36ed9..cb6d02c925a2736334e5110b99580da1b919a2a2 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -42,38 +42,38 @@ class GroupConstructor():
if mapMesh:
# Save initial mesh.
- cMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
- cMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
+ cMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
+ cMeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][:,0])]
cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
- cxxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][1:,8],dataIn[1][:,8]))
+ cxxExt = np.concatenate((dataIn[0][0][:,8], dataIn[0][1][:,8],dataIn[1][:,8]))
# Create fine mesh.
fMeshUp = self.createFineMesh(dataIn[0][0][:,0], nFactor)
fMeshLw = self.createFineMesh(dataIn[0][1][:,0], nFactor)
fMeshWk = self.createFineMesh(dataIn[1][:,0] , nFactor)
- fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
- fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
+ fMesh = np.concatenate((fMeshUp, fMeshLw, fMeshWk))
+ fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw)]
fxx_up = self.interpolateToFineMesh(xx_up, fMeshUp, nFactor)
fxx_lw = self.interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
fxx_wk = self.interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
- fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
+ fxx = np.concatenate((fxx_up, fxx_lw, fxx_wk))
fvtExt_up = self.interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
fvtExt_lw = self.interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
fvtExt_wk = self.interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
- fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
+ fvtExt = np.concatenate((fvtExt_up, fvtExt_lw, fvtExt_wk))
fMe_up = self.interpolateToFineMesh(dataIn[0][0][:,5], fMeshUp, nFactor)
fMe_lw = self.interpolateToFineMesh(dataIn[0][1][:,5], fMeshLw, nFactor)
fMe_wk = self.interpolateToFineMesh(dataIn[1][:,5] , fMeshWk, nFactor)
- fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
+ fMe = np.concatenate((fMe_up, fMe_lw, fMe_wk))
frho_up = self.interpolateToFineMesh(dataIn[0][0][:,6], fMeshUp, nFactor)
frho_lw = self.interpolateToFineMesh(dataIn[0][1][:,6], fMeshLw, nFactor)
frho_wk = self.interpolateToFineMesh(dataIn[1][:,6] , fMeshWk, nFactor)
- frho = np.concatenate((frho_up, frho_lw[1:], frho_wk))
+ frho = np.concatenate((frho_up, frho_lw, frho_wk))
fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
@@ -83,17 +83,17 @@ class GroupConstructor():
fxxExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
fxxExt_wk = self.interpolateToFineMesh(dataIn[1][:,8] , fMeshWk, nFactor)
- fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
- fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
+ fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw, fdStarExt_wk))
+ fxxext = np.concatenate((fxxExt_up, fxxExt_lw, fxxExt_wk))
- fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
+ fdv = np.concatenate((dv_up, dv_lw, dv_wk))
fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
- cMe = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][1:,5], dataIn[1][:,5]))
- cRho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][1:,6], dataIn[1][:,6]))
- cdStarExt = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][1:,7], dataIn[1][:,7]))
+ cMe = np.concatenate((dataIn[0][0][:,5], dataIn[0][1][:,5], dataIn[1][:,5]))
+ cRho = np.concatenate((dataIn[0][0][:,6], dataIn[0][1][:,6], dataIn[1][:,6]))
+ cdStarExt = np.concatenate((dataIn[0][0][:,7], dataIn[0][1][:,7], dataIn[1][:,7]))
cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
@@ -112,7 +112,7 @@ class GroupConstructor():
else:
- Mesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ Mesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][:,0], dataIn[1][:,0]))
MeshbNodes = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0]) + len(dataIn[0][1][1:,0])]
xx = np.concatenate((xx_up, xx_lw[1:], xx_wk))