diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index b7acfff4d3aa5a69bc1fda57fe636e31ca81a21f..ffe14f8058309d6f5b2462101488514518308589 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -31,3 +31,4 @@
#include "wTimeSolver.h"
#include "wViscFluxes.h"
#include "wClosures.h"
+#include "wSolutionInterp.h"
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index 4e12b24f3c3e837d9e707aa53e39885a56f8d2dc..29154d99b35476438d3686513978d8055ee9d180 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -63,6 +63,7 @@ threads="1"
%include "wTimeSolver.h"
%include "wViscFluxes.h"
%include "wClosures.h"
+%include "wSolutionInterp.h"
%include "wFluid.h"
%include "wAssign.h"
diff --git a/dart/default.py b/dart/default.py
index 51b3757950f3a815b34009c09fd726d09e76179b..81e5153ec78fac11692a4e555ca09c384e32442f 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -23,8 +23,8 @@ from fwk.wutils import parseargs
import dart
import tbox
-def initBL(Re, Minf, CFL0):
- solver = dart.ViscSolver(Re, Minf, CFL0)
+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):
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 3a2af7174074169a14ca4270371d7a08eaabbf2d..21a46657c2baf959cd560e78dded1dda357e26bb 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -69,6 +69,7 @@ class BLRegion;
class TimeSolver;
class ViscSolver;
class Closures;
+class SolutionInterp;
// General
class F0Ps;
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index d9c89a69c04e4a1d6fbc2e2e5f5010e618f5a974..8d001e8b4c5c4357a039c6840727aad942786380 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -96,6 +96,7 @@ void BLRegion::SetStagBC(double Re)
std::vector<double> ClosParam = closSolver->LaminarClosures(U[0], U[1], Ret[0], Me[0], name);
+ std::cout << ClosParam[3] << std::endl;
Hstar[0] = ClosParam[0];
Hstar2[0] = ClosParam[1];
Hk[0] = ClosParam[2];
@@ -141,7 +142,6 @@ void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<d
{
Regime[k] = 1;
}
- PrintSolution(0);
}
void BLRegion::ComputeBlwVel()
@@ -175,4 +175,4 @@ void BLRegion::PrintSolution(size_t iPoint)
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/wSolutionInterp.cpp b/dart/src/wSolutionInterp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ca13de009747b8501ca8d2d696b80dda27d2abb7
--- /dev/null
+++ b/dart/src/wSolutionInterp.cpp
@@ -0,0 +1,97 @@
+#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)
+{
+ if (meshFactor<=1)
+ {
+ std::cout << "Can not create finer mesh while meshFactor is " << meshFactor << std::endl;
+ return cVector;
+ }
+ else
+ {
+ std::vector<double> fVector;
+ 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.push_back(cVector[cPoint] + fPoint * dv / meshFactor);
+ }
+ }
+ std::cout << cVector.back() << std::endl;
+ fVector.push_back(cVector.back());
+ return fVector;
+ }
+}
+
+std::vector<double> SolutionInterp::CreateFineMesh(std::vector<double> cMesh)
+{
+ if (meshFactor<1)
+ {
+ std::cout << "SolutionInterp::CreateFineMesh: Mesh factor should be > 1." << std::endl;
+ /*std::stringstream err;
+ err << "SolutionInterp::CreateFineMesh: Mesh factor should be > 1. Value: " << meshFactor << "\n";
+ throw std::runtime_error(err.str());*/
+ }
+
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..d023d57fed99e71dbceba2d6c625af526a514bd9
--- /dev/null
+++ b/dart/src/wSolutionInterp.h
@@ -0,0 +1,33 @@
+#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> CreateFineMesh(std::vector<double> cMesh);
+ void ResetmeshFactor(unsigned int newmeshFactor);
+ std::vector<double> InterpolateToFineMesh(std::vector<double> cVector);
+ 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 77ea0594c9f323cdcec41204f0d0ef201a223559..e0d17f706a6bb296dbd435f86d01478dc1d61a4e 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -120,14 +120,8 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
bl->U[iPoint*nVar + 1] = std::max(bl->U[iPoint*nVar + 1], 1.0005);
}
-
SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
normSysRes = (SysRes + slnIncr/timeStep).norm();
- if (bl->name == "wake" && iPoint == 20)
- {
- std::cout << normSysRes/normSysRes0 << std::scientific;
- bl->PrintSolution(iPoint);
- }
if (normSysRes <= tol * normSysRes0)
{
@@ -192,7 +186,6 @@ double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<d
{
if (std::isnan(bl->U[iPoint*nVar + k]))
{
- //std::cout << "Point " << iPoint << " : Solution reset; NaN identified at position " << k << std::endl;
ResetSolution(iPoint, bl, temp_U);
itFrozenJac0 = 1; // Impose continuous Jacobian evaluation.
nErrors+=1;
@@ -204,7 +197,6 @@ double TimeSolver::ControlIntegration(size_t iPoint, BLRegion *bl, std::vector<d
void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double> temp_U){
- //std::cout << "Reseting point" << iPoint << std::endl;
unsigned int nVar = bl->GetnVar();
/* Reset solution */
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 9d0e577eeda60896b732c8c2c9ed2151b28cc9a9..bce4235b212fcd20640fa96ed78ed263cf7a00e3 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -1,23 +1,45 @@
+/* ----------------------------------------------------------------------------------- */
+/* */
+/* Coupled integral boundary layer solver */
+/* */
+/* Université de Liège */
+/* February 2022 */
+/* Author: Paul Dechamps */
+/* v0.1 */
+/* ----------------------------------------------------------------------------------- */
+
+/* --- Project includes -------------------------------------------------------------- */
+
#include "wViscSolver.h"
#include "wBLRegion.h"
#include "wTimeSolver.h"
+#include "wSolutionInterp.h"
#include <iostream>
#include <vector>
+/* --- Namespaces -------------------------------------------------------------------- */
using namespace tbox;
using namespace dart;
+/* ----------------------------------------------------------------------------------- */
-ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
+/**
+ * @brief Viscous solver and boundary layer(s) initialization
+ */
+ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor)
{
+ PrintBanner();
+
/* Freestream parameters */
- Re = _Re;
- Minf = _Minf;
+
+ Re = _Re; /* Freestream Reynolds number */
+ Minf = _Minf; /* Freestream Mach number */
/* User input numerical parameters */
- CFL0 = _CFL0;
+
+ CFL0 = _CFL0; /* Starting CFL number */
std::cout << "--- Viscous solver setup ---\n"
<< "Reynolds number: " << Re << " \n"
@@ -33,20 +55,34 @@ ViscSolver::ViscSolver(double _Re, double _Minf, double _CFL0)
bl.push_back(region);
} // end for
- tSolver = new TimeSolver(CFL0, Minf, Re);
-
bl[0]->name = "upper";
bl[1]->name = "lower";
bl[2]->name = "wake";
+ /* Temporal solver initialization */
+
+ tSolver = new TimeSolver(CFL0, Minf, Re);
+
+ /* Mesh conversion setup */
+
+ if (meshFactor > 1)
+ {
+ mapMesh = true;
+ }
+ meshConverter = new SolutionInterp(meshFactor);
+
+ /* Stagnation point movements monitoring and control */
+
nbElmUpper = 0;
stagPtMvmt = false;
- std::cout << "Initializing boundary layer" << std::endl;
+ std::cout << "Boundary layer initialized" << std::endl;
} // end ViscSolver
-
+/**
+ * @brief Viscous solver and subsequent dependancies destruction
+ */
ViscSolver::~ViscSolver()
{
for (size_t i=0; i<bl.size(); ++i){
@@ -54,30 +90,77 @@ ViscSolver::~ViscSolver()
} // end for
delete tSolver;
+ delete meshConverter;
- std::cout << "~ViscSolver()\n";
-
+ PrintBanner();
+ std::cout << "~ViscSolver()\n";
} // end ~ViscSolver
+/**
+ * @brief Initialize 1D mesh of specified region.
+ */
void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM)
{
- bl[region]->SetMesh(locM, globM);
+ if (mapMesh)
+ {
+ std::vector<double> fLocMesh = meshConverter->CreateFineMesh(locM);
+ std::vector<double> fGlobMesh = meshConverter->CreateFineMesh(globM);
+ 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);
+ 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)
{
- bl[region]->SetInvBC(_Ue, _Me, _Rhoe);
+
+ if (mapMesh)
+ {
+ std::vector<double> UeFine = meshConverter->InterpolateToFineMesh(_Ue);
+ std::vector<double> MeFine = meshConverter->InterpolateToFineMesh(_Me);
+ std::vector<double> RhoeFine = meshConverter->InterpolateToFineMesh(_Rhoe);
+
+ bl[region]->SetInvBC(UeFine, MeFine, RhoeFine);
+ std::cout << "- " << bl[region]->name << ": Inviscid boundary interpolated." << std::endl;
+ }
+ else
+ {
+ bl[region]->SetInvBC(_Ue, _Me, _Rhoe);
+ }
}
+/**
+ * @brief Set mesh and displacement thickness at the edge of the boundary layer.
+ * Used for the interaction law.
+ */
void ViscSolver::SendExtVariables(size_t region, std::vector<double> _xxExt, std::vector<double> _deltaStarExt)
{
- bl[region]->SetExtVariables(_xxExt, _deltaStarExt);
+ if (mapMesh)
+ {
+ std::vector<double> xxExtFine = meshConverter->InterpolateToFineMesh(_xxExt);
+ std::vector<double> deltaStarExtFine = meshConverter->InterpolateToFineMesh(_deltaStarExt);
+ 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 */
@@ -98,7 +181,7 @@ int ViscSolver::Run(unsigned int couplIter){
else{
tSolver->SetCFL0(0.5);
- tSolver->SetitFrozenJac(1);
+ tSolver->SetitFrozenJac(5);
tSolver->SetinitSln(false);
stagPtMvmt = false;
std::cout << "Updating solution" << std::endl;
@@ -190,6 +273,9 @@ int ViscSolver::Run(unsigned int couplIter){
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. */
@@ -205,6 +291,9 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
}
}
+/**
+ * @brief Computes turbulent solution @ the transition point and averages solutions.
+ */
void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
{
/* Averages solution on transition marker */
@@ -276,14 +365,28 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->Us[iPoint] = turbParam[7];
}
+/**
+ * @brief Extracts blowing velocity distribution of the corresponding region.
+ */
std::vector<double> ViscSolver::ExtractBlowingVel(size_t iRegion)
{
return bl[iRegion]->Blowing;
}
+
+/**
+ * @brief Extracts displacement thickness distribution of the corresponding region.
+ */
std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
{
return bl[iRegion]->DeltaStar;
}
+
+/**
+ * @brief Extracts mesh distribution of the corresponding region.
+ *
+ * @param iRegion Region ID.
+ * @return LocMarkersOut Local markers of the BL mesh.
+ */
std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
{
std::vector<double> LocMarkersOut;
@@ -293,7 +396,13 @@ std::vector<double> ViscSolver::ExtractXx(size_t iRegion)
return LocMarkersOut;
}
-
+/**
+ * @brief Friction, pressure and total drag computation of the configuration.
+ *
+ * Cdtot = theta * Ue^((H+5)/2).
+ * Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
+ * Cdp = Cdtot - Cdf.
+ */
void ViscSolver::ComputeDrag()
{
@@ -332,4 +441,70 @@ void ViscSolver::ComputeDrag()
/* Compute pressure drag coefficient. */
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(17);
+ 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)
+ {
+ 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
+ }
+ }
+ return blVectors;
+}
+
+std::vector<double> ViscSolver::Getxtr()
+{
+ 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;
+
}
\ No newline at end of file
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index c1e6b07ffa77e2bc916dd3e728a61b6484dcec08..a8286b46153e204b0b5658c40e950e55537fe3f1 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -20,13 +20,15 @@ private:
double CFL0;
size_t nbElmUpper;
bool stagPtMvmt;
+ bool mapMesh=false;
TimeSolver *tSolver;
+ SolutionInterp *meshConverter;
public:
std::vector<BLRegion *> bl;
- ViscSolver(double _Re, double _Minf, double _CFL0);
+ ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor);
~ViscSolver();
void InitMeshes(size_t region, std::vector<double> locM, std::vector<double> globM);
@@ -39,13 +41,15 @@ public:
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();
};
} // namespace dart
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 551f771a908724a19d8cfe5b123c92d8f0630fe7..24c59b0b175398c2cf3b67b13faf676295bdc0ae 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -60,7 +60,7 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5.*math.pi/180
+ alpha = 2.*math.pi/180
M_inf = 0.
#user_xtr=[0.41,0.41]
@@ -86,7 +86,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
#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()
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index e5c52372fbcf8fbf97a7ea24fe2e81f28de53b35..0046b19aaae268d2fa76610b12c312b66da85145 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -36,6 +36,7 @@
import dart.utils as floU
import dart.default as floD
+import dart.viscUtils as viscU
import dart.vCoupler as floVC
@@ -56,33 +57,22 @@ def main():
# define flow variables
Re = 1e7
alpha = 2*math.pi/180
- M_inf = 0.
+ M_inf = 0.2
- #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
+ meshFactor = 1
CFL0 = 1
+ plotVar = 'Hk'
+
# ========================================================================================
- 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}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -97,7 +87,7 @@ def main():
tms['solver'].start()
isolver = floD.newton(pbl)
- vsolver = floD.initBL(Re, M_inf, CFL0)
+ vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
coupler = floVC.Coupler(isolver, vsolver)
@@ -126,8 +116,10 @@ def main():
#floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
- #val=validation.Validation(case)
- #val.main(isolver.Cl,vsolver.wData)
+ blScalars, blVectors = viscU.ExtractVars(vsolver)
+ xtr=vsolver.Getxtr()
+ wData = viscU.Dictionary(blScalars, blVectors, xtr)
+ viscU.PlotVars(wData, plotVar)
# check results
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index 21f29243ac35f7ffafea0dc08b832d61d6899613..0884c4cda7679d426998c48f022d663160369d62 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -25,6 +25,7 @@ from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
from fwk.coloring import ccolors
+import fwk
import numpy as np
from matplotlib import pyplot as plt
@@ -35,6 +36,8 @@ class Coupler:
self.maxCouplIter = 100
self.couplTol = 1e-4
+
+ self.tms=fwk.Timers()
pass
def CreateProblem(self, _boundaryAirfoil, _boundaryWake):
@@ -56,7 +59,9 @@ class Coupler:
# Run inviscid solver
+ self.tms['inv'].start()
self.iSolver.run()
+ self.tms['inv'].stop()
# Extract inviscid boundary
@@ -99,12 +104,16 @@ class Coupler:
print('Inviscid boundary imposed')
# Run viscous solver
+ self.tms['visc'].start()
exitCode = self.vSolver.Run(couplIter)
+ self.tms['visc'].stop()
print('Viscous ops terminated with exit code ', exitCode)
error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
if error <= self.couplTol and exitCode==0:
+
+ print(self.tms)
print('')
print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
.format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
diff --git a/dart/viscUtils.py b/dart/viscUtils.py
new file mode 100644
index 0000000000000000000000000000000000000000..7cb75d80d32bcc20beaf49ca660da4ff3de1e227
--- /dev/null
+++ b/dart/viscUtils.py
@@ -0,0 +1,92 @@
+
+
+
+def Dictionary(blScalars, blVectors, xtr):
+ """Create dictionary with the coordinates and the boundary layer parameters
+ """
+
+ wData = {}
+ # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
+ wData['xx'] = blVectors[0]
+ wData['x/c'] = blVectors[1]
+ wData['deltaStar'] = blVectors[2]
+ wData['theta'] = blVectors[3]
+ wData['H'] = blVectors[4]
+ wData['N'] = blVectors[5]
+ wData['Ue'] = blVectors[6]
+ wData['Ct'] = blVectors[7]
+ wData['Hk'] = blVectors[8]
+ wData['Hstar'] = blVectors[9]
+ wData['HStar2'] = blVectors[10]
+ wData['Cf'] = blVectors[11]
+ wData['Cd'] = blVectors[12]
+ wData['delta'] = blVectors[13]
+ wData['VtExt'] = blVectors[14]
+ wData['Me'] = blVectors[15]
+ wData['Rhoe'] = blVectors[16]
+ wData['xtrTop'] = xtr[0]
+ wData['xtrBot'] = xtr[1]
+ return wData
+
+def WriteFile(wData, Re, toW=['delta*', 'H', 'Cf', 'Ctau']):
+ """Write the results in airfoil_viscous.dat
+ """
+
+ # Write
+ print('Writing file: airfoil_viscous.dat...', end = '')
+ f = open('airfoil_viscous.dat', 'w+')
+
+ f.write('$Aerodynamic coefficients\n')
+ f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(Re, wData['Cd'], wData['Cdp'],
+ wData['Cdf'], wData['xtrTop'],
+ wData['xtrBot']))
+ f.write('$Boundary layer variables\n')
+ f.write(' x y z x/c')
+
+ for s in toW:
+ f.write(' {0:>15s}'.format(s))
+ f.write('\n')
+
+ for i in range(len(self.data[:,0])):
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
+ for s in toW:
+ f.write(' {0:15.6f}'.format(wData[s][i]))
+ f.write('\n')
+
+ f.close()
+ print('done.')
+
+def ExtractVars(viscSolver):
+ import numpy as np
+
+ blScalars = [viscSolver.Cdt, viscSolver.Cdf, viscSolver.Cdp]
+
+ # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
+
+ blVectors = viscSolver.ExtractSolution()
+
+
+ return blScalars, blVectors
+
+def PlotVars(wData, var='H'):
+ from matplotlib import pyplot as plt
+
+ nUpper = len(wData['x/c'])
+ for i in range(len(wData['x/c'])):
+ if wData['x/c'][i] == 1.0:
+ nUpper = i+1
+ break
+
+ # [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]
+ # [xx, x, deltaStar, theta, H, N, Ue, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta, VtExt, Me, Rhoe]
+
+ plt.plot(wData['x/c'][:nUpper], wData[var][:nUpper], color='#0072BD', linewidth=3)
+ plt.plot(wData['x/c'][nUpper:], wData[var][nUpper:], color='#D95319', linewidth=3)
+ plt.axvline(x=wData['xtrTop'], linestyle="--", color='#0072BD')
+ plt.axvline(x=wData['xtrBot'], linestyle="--", color='#D95319')
+ plt.xlabel('$x/c$')
+ plt.ylabel(var)
+ plt.xlim([0,1])
+ plt.legend(['Upper side', 'Lower side'], frameon=False)
+ plt.show()
\ No newline at end of file
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index 4d83869c771a6eec9dbaf8dca1f3a9adbbe7eefa..b64b9ee563b5a3fce91e4e06209f117c35c3ac1b 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -23,6 +23,7 @@ from dart.viscous.Master.DAirfoil import Airfoil
from dart.viscous.Master.DWake import Wake
from matplotlib import pyplot as plt
+import fwk
import numpy as np
import tbox.utils as tboxU
@@ -39,6 +40,8 @@ class Coupler:
self.tol = _tol # tolerance of the VII
self.writer = _writer
+ self.tms=fwk.Timers()
+
def run(self):
"""Viscous inviscid coupling.
"""
@@ -58,8 +61,10 @@ class Coupler:
# Run inviscid solver
print('+----------------------------- Inviscid solver --------------------------------+')
+ self.tms['inv'].start()
self.isolver.reset()
self.isolver.run()
+ self.tms['inv'].stop()
for n in range(0, len(self.group)):
@@ -79,7 +84,9 @@ class Coupler:
# Run viscous solver
print('+------------------------------ Viscous solver --------------------------------+')
+ self.tms['visc'].start()
self.vsolver.run(self.group)
+ self.tms['visc'].stop()
for n in range(0, len(self.group)):
for i in range(0, self.group[n].nE):
@@ -91,6 +98,7 @@ class Coupler:
# Check convergence and output coupling iteration.
if error <= self.tol:
+ print(self.tms)
print('')
print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
.format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)