From 7b0a6935239583048879fcb0021cc881a385bd9d Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paulzer@Pauls-MacBook-Pro.local>
Date: Wed, 20 Apr 2022 03:17:53 +0200
Subject: [PATCH] Refined to find error
---
dart/pyVII/vCoupler.py | 22 ++--
dart/pyVII/vCoupler2.py | 2 +-
dart/src/wBLRegion.cpp | 49 +++++----
dart/src/wBLRegion.h | 2 +-
dart/src/wClosures.cpp | 17 +--
dart/src/wTimeSolver.cpp | 33 +++---
dart/src/wTimeSolver.h | 18 ++--
dart/src/wViscFluxes.cpp | 35 +++---
dart/src/wViscFluxes.h | 2 +-
dart/src/wViscMesh.cpp | 13 ++-
dart/src/wViscMesh.h | 3 +-
dart/src/wViscSolver.cpp | 120 ++++++++++-----------
dart/src/wViscSolver.h | 2 +-
dart/tests/lift.py | 17 ++-
dart/viscous/Drivers/DGroupConstructor.py | 8 +-
dart/viscous/Drivers/DGroupConstructor2.py | 8 +-
dart/viscous/Physics/DDataStructures.py | 4 +-
17 files changed, 188 insertions(+), 167 deletions(-)
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 9de6006..8876da9 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -36,7 +36,7 @@ class Coupler:
self.iSolver=iSolver
self.vSolver=vSolver
- self.maxCouplIter = 100
+ self.maxCouplIter = 150
self.couplTol = 1e-4
self.tms=fwk.Timers()
@@ -88,7 +88,7 @@ class Coupler:
# print('Viscous ops terminated with exit code ', exitCode)
self.tms['processing'].start()
- self.__ImposeBlwVel()
+ self.__ImposeBlwVel(couplIter)
self.tms['processing'].stop()
error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
@@ -107,7 +107,7 @@ class Coupler:
xtr = [1, 1]
xtr[0]=self.vSolver.Getxtr(0,0)
xtr[1]=self.vSolver.Getxtr(0,1)
- print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter,
+ print(' {:>4.0f}| {:>10.10f} {:>10.10f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter,
self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
couplIter += 1
@@ -124,7 +124,7 @@ class Coupler:
self.iSolver.printOn = False
self.vSolver.printOn = False
- coeffTable = np.empty((3,len(alphaSeq)))
+ coeffTable = np.zeros((3,len(alphaSeq)))
coeffTable[0,:] = alphaSeq
for ind, alpha in enumerate(alphaSeq):
@@ -139,10 +139,9 @@ class Coupler:
return coeffTable
- def __ImposeBlwVel(self):
+ def __ImposeBlwVel(self, couplIter):
# Retreive and set blowing velocities.
-
for iRegion in range(3):
self.blwVel[iRegion] = np.asarray(self.vSolver.ExtractBlowingVel(0, iRegion))
self.deltaStar[iRegion] = np.asarray(self.vSolver.ExtractDeltaStar(0, iRegion))
@@ -163,9 +162,6 @@ class Coupler:
self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
- """
- plt.plot(self.group[0].u)
- plt.show()"""
for n in range(0, len(self.group)):
for i in range(0, self.group[n].nE):
@@ -190,9 +186,10 @@ class Coupler:
if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter==0):
- # print('STAG POINT MOVED')
+ print('STAG POINT MOVED')
# Initialize mesh
+
self.vSolver.SetGlobMesh(0, 0, fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['yGlobCoord'][:fMeshDict['bNodes'][1]], fMeshDict['zGlobCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.SetGlobMesh(0, 1, fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.SetGlobMesh(0, 2, fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][2]:], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][2]:])
@@ -241,8 +238,9 @@ class Coupler:
plt.ylim([-0.3, 0.2])
plt.title('Uy = f(x)')
plt.show()"""
-
-
+ """plt.plot(fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['vx'][:fMeshDict['bNodes'][1]])
+ plt.plot(fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
+ plt.show()"""
self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.SetVelocities(0, 0, fMeshDict['vx'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]])
self.vSolver.SetVelocities(0, 1, fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
diff --git a/dart/pyVII/vCoupler2.py b/dart/pyVII/vCoupler2.py
index b388662..6292c25 100644
--- a/dart/pyVII/vCoupler2.py
+++ b/dart/pyVII/vCoupler2.py
@@ -114,7 +114,7 @@ class Coupler:
self.iSolver.printOn = False
self.vSolver.printOn = False
- coeffTable = np.empty((3,len(alphaSeq)))
+ coeffTable = np.zeros((3,len(alphaSeq)))
coeffTable[0,:] = alphaSeq
for ind, alpha in enumerate(alphaSeq):
diff --git a/dart/src/wBLRegion.cpp b/dart/src/wBLRegion.cpp
index 43a776e..0998854 100644
--- a/dart/src/wBLRegion.cpp
+++ b/dart/src/wBLRegion.cpp
@@ -6,13 +6,16 @@
#include <vector>
#include <cmath>
+#include <iomanip>
+#include <deque>
+
using namespace dart;
BLRegion::BLRegion(double _xtrF)
{
if (_xtrF < 0. && _xtrF != -1.)
{
- std::cout << "Unrecognized transition location " << _xtrF << std::endl;
+ std::cout << "Fatal error: Invalid transition location " << _xtrF << std::endl;
throw;
}
xtrF = _xtrF;
@@ -20,7 +23,7 @@ BLRegion::BLRegion(double _xtrF)
mesh = new ViscMesh();
invBnd = new InvBnd();
- transMarker = mesh->GetnMarkers();
+ transMarker = 0;
xtr = 1.0;
} // end BLRegion
@@ -56,7 +59,7 @@ void BLRegion::SetMesh(std::vector<double> x,
Us.resize(nMarkers,0.);
DeltaStar.resize(nMarkers,0.);
Delta.resize(nMarkers,0.);
- Regime.resize(nMarkers,0.);
+ Regime.resize(nMarkers, 0);
}
} // End SetMesh
@@ -73,15 +76,15 @@ void BLRegion::SetStagBC(double Re)
U[0] = 1e-6; // Theta
} // end else
U[1] = 2.23; // H
- U[2] = 0; // N
+ U[2] = 0.; // N
U[3] = invBnd->GetVt(0); // ue
- U[4] = 0; // Ctau
+ U[4] = 0.; // Ctau
Ret[0] = U[3] * U[0] * Re;
if (Ret[0] < 1)
{
- Ret[0] = 1;
+ Ret[0] = 1.;
U[3] = Ret[0] / (U[0] * Re);
} // End if
@@ -104,10 +107,13 @@ void BLRegion::SetStagBC(double Re)
void BLRegion::SetWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
{
-
+ if (name!="wake")
+ {
+ std::cout << "Warning: Imposing wake boundary condition on " << name << std::endl;
+ }
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[2] = 9.;
U[3] = invBnd->GetVt(0); /* Inviscid velocity. */
U[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / U[0]; /* ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low). */
@@ -154,18 +160,17 @@ void BLRegion::ComputeBlwVel()
void BLRegion::PrintSolution(size_t iPoint) const
{
- 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;
+ std::cout << "Pt " << iPoint << "Reg " << Regime[iPoint] << std::endl;
+ std::cout << std::scientific << std::setprecision(15);
+ std::cout << std::setw(10) << "T " << U[iPoint*nVar + 0] << "\n"
+ << std::setw(10) << "H " << U[iPoint*nVar + 1] << "\n"
+ << std::setw(10) << "N " << U[iPoint*nVar + 2] << "\n"
+ << std::setw(10) << "ue " << U[iPoint*nVar + 3] << "\n"
+ << std::setw(10) << "Ct " << U[iPoint*nVar + 4] << "\n"
+ << std::setw(10) << "dStar (BL) " << DeltaStar[iPoint] << "\n"
+ << std::setw(10) << "dStar (old) " << invBnd->GetDeltaStar(iPoint) << "\n"
+ << std::setw(10) << "Vt " << invBnd->GetVt(iPoint) << "\n"
+ << std::setw(10) << "Me " << invBnd->GetMe(iPoint) << "\n"
+ << std::setw(10) << "Rhoe " <<invBnd->GetRhoe(iPoint) << "\n"
+ << "\n";
}
\ No newline at end of file
diff --git a/dart/src/wBLRegion.h b/dart/src/wBLRegion.h
index 0e2cc5e..00e111c 100644
--- a/dart/src/wBLRegion.h
+++ b/dart/src/wBLRegion.h
@@ -15,7 +15,7 @@ class DART_API BLRegion
{
private:
- unsigned int nVar = 5;
+ unsigned int const nVar = 5;
/* Transition related variables */
diff --git a/dart/src/wClosures.cpp b/dart/src/wClosures.cpp
index 669c38e..f9874c4 100644
--- a/dart/src/wClosures.cpp
+++ b/dart/src/wClosures.cpp
@@ -86,7 +86,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
if (name == "wake")
{
Cd = 2*(1.10*(1.0-1.0/Hk)*(1.0-1.0/Hk)/Hk)* (Hstar/(2*Ret));
- Cf = 0;
+ Cf = 0.;
}
double Us = 0.;
@@ -108,7 +108,7 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
if (std::isnan(ClosureVars[i]))
{
std::cout << "lam : nan detected in closureVars at position " << i << std::endl;
- throw;
+ //throw;
}
}
return ClosureVars;
@@ -116,7 +116,10 @@ std::vector<double> Closures::LaminarClosures(double theta, double H, double Ret
std::vector<double> Closures::TurbulentClosures(double theta, double H, double Ct, double Ret, double Me, std::string name) const
{
-
+ if (std::isnan(Ret))
+ {
+ std::cout << "Ret is nan. H:" << H << std::endl;
+ }
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
//Ct = std::max(std::min(Ct, 0.30), 0.0000001);
@@ -165,8 +168,8 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
double logRt = std::log(Ret/Fc);
if (std::isnan(logRt))
{
- std::cout << "logRt is nan " << std::endl;
- throw;
+ std::cout << "logRt is nan. Ret= " << Ret << std::endl;
+ //throw;
}
logRt = std::max(logRt, 3.0);
double arg = std::max(-1.33*Hk, -20.0);
@@ -226,7 +229,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
if (n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H)<0)
{
std::cout << "Negative sqrt encountered " << std::endl;
- throw;
+ //throw;
}
double Cteq = std::sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc*Hkc)/((1-Us)*(Hk*Hk)*H));
Cd = n*(Cdw+ Cdd + Cdl);
@@ -248,7 +251,7 @@ std::vector<double> Closures::TurbulentClosures(double theta, double H, double C
if (std::isnan(ClosureVars[i]))
{
std::cout << "turb : nan detected in closureVars at position " << i << std::endl;
- throw;
+ //throw;
}
}
return ClosureVars;
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index ae69b24..42870da 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -19,6 +19,7 @@ TimeSolver::TimeSolver(double _CFL0, double _Minf, double Re, unsigned long _max
tol = _tol;
itFrozenJac0 = _itFrozenJac;
Minf = std::max(_Minf, 0.1);
+ initSln = true;
SysMatrix = new ViscFluxes(Re);
@@ -95,16 +96,20 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
VectorXd slnIncr = VectorXd::Zero(5);
- nErrors = 0; // Number of errors encountered
+ unsigned int nErrors = 0; // Number of errors encountered
unsigned int innerIters = 0; // Inner (non-linear) iterations
+ /* for (int i=0; i<1e10; ++i)
+ {
+ double a = 1+1;
+ } */
while (innerIters < maxIt)
{
- if (bl->name=="upper" && iPoint == 1)
+ /* if (bl->name=="upper" && iPoint == 1 && innerIters == 0)
{
std::cout << "it " << innerIters << " res " << std::log10(normSysRes/normSysRes0) << std::endl;
- bl->PrintSolution(iPoint);
- }
+ bl->PrintSolution(0);
+ } */
if (innerIters % itFrozenJac == 0)
{
@@ -119,6 +124,10 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
}
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);
+ if (std::isnan(bl->U[iPoint * nVar + 1]))
+ {
+ std::cout << iPoint << " nan at it " << innerIters << std::endl;
+ }
SysRes = SysMatrix->ComputeFluxes(iPoint, bl);
normSysRes = (SysRes + slnIncr / timeStep).norm();
@@ -132,7 +141,7 @@ int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
normSysRes = (SysRes + slnIncr / timeStep).norm();
itFrozenJac = 1;
std::cout << "normSysRes is nan " << std::endl;
- throw;
+ //throw;
}
if (normSysRes <= tol * normSysRes0)
@@ -197,6 +206,9 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
/* Reset closures */
+ bl->Ret[iPoint] = std::max( bl->U[iPoint * nVar + 3] * bl->U[iPoint * nVar + 0] * SysMatrix->Re, 1.0); /* Reynolds number based on theta */
+ bl->DeltaStar[iPoint] = bl->U[iPoint * nVar + 0] * bl->U[iPoint * nVar + 1]; /* Displacement thickness */
+
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->invBnd->GetMe(iPoint), bl->name);
@@ -222,13 +234,4 @@ void TimeSolver::ResetSolution(size_t iPoint, BLRegion *bl, std::vector<double>
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; }
-
-void TimeSolver::SetinitSln(bool _initSln) { initSln = _initSln; }
\ No newline at end of file
+} // End ResetSolution
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.h b/dart/src/wTimeSolver.h
index 5866372..6ab52df 100644
--- a/dart/src/wTimeSolver.h
+++ b/dart/src/wTimeSolver.h
@@ -18,9 +18,8 @@ private:
double tol;
unsigned int itFrozenJac0;
bool initSln;
- double gamma = 1.4;
+ double const gamma = 1.4;
double Minf;
- unsigned int nErrors;
ViscFluxes *SysMatrix;
@@ -29,13 +28,14 @@ public:
~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);
+
+ double GetCFL0() const {return CFL0;};
+ void SetCFL0(double _CFL0) {if(_CFL0>0){CFL0 = _CFL0;}};
+ unsigned long GetmaxIt() const {return maxIt;};
+ void SetmaxIt(unsigned long _maxIt) {maxIt = _maxIt;};
+ unsigned int GetitFrozenJac() const {return itFrozenJac0;};
+ void SetitFrozenJac(unsigned int _itFrozenJac) {itFrozenJac0 = _itFrozenJac;};
+ void SetinitSln(bool _initSln) {initSln = _initSln;};
private:
double SetTimeStep(double CFL, double dx, double invVel) const;
diff --git a/dart/src/wViscFluxes.cpp b/dart/src/wViscFluxes.cpp
index 5d039e5..34334a5 100644
--- a/dart/src/wViscFluxes.cpp
+++ b/dart/src/wViscFluxes.cpp
@@ -72,9 +72,14 @@ VectorXd ViscFluxes::BLlaws(size_t iPoint, BLRegion *bl, std::vector<double> u)
} // End if
else
{
- dissipRatio = 1;
+ dissipRatio = 1.;
} // End else
+ if (std::isnan(u[1]))
+ {
+ std::cout << "H is nan. Point " << iPoint << bl->name << std::endl;
+ }
+
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 */
@@ -245,22 +250,22 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
double Dgr = 0.08;
double Hk1 = 3.5;
- double Hk2 = 4;
- double Hmi = (1 / (Hk - 1));
+ 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)
+ if (Ret <= 0.)
{
- Ret = 1;
+ Ret = 1.;
}
- if (log10(Ret) < logRet_crit - Dgr)
+ if (std::log10(Ret) < logRet_crit - Dgr)
{
- Ax = 0;
+ Ax = 0.;
}
else
{
- double Rnorm = (log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
+ double Rnorm = (std::log10(Ret) - (logRet_crit - Dgr)) / (2 * Dgr);
double Rfac = 0.;
- if (Rnorm <= 0)
+ if (Rnorm <= 0.)
{
Rfac = 0;
}
@@ -273,16 +278,16 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
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 m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3 * Hmi * Hmi * Hmi + 0.1 * std::exp(-20. * Hmi);
double arg = 3.87 * Hmi - 2.52;
- double dndRet = 0.028 * (Hk - 1) - 0.0345 * exp(-(arg * arg));
+ double dndRet = 0.028 * (Hk - 1) - 0.0345 * std::exp(-(arg * arg));
Ax = (m * dndRet / theta) * Rfac;
// Set correction for separated profiles
if (Hk > Hk1)
{
double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
double Hfac = 0.;
- if (Hnorm >= 1)
+ if (Hnorm >= 1.)
{
Hfac = 1;
}
@@ -291,9 +296,9 @@ double ViscFluxes::AmplificationRoutine(double Hk, double Ret, double theta) con
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;
+ double Gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5));
+ double Tnr = std::tanh(1.2 * (std::log10(Ret) - Gro));
+ double Ax2 = (0.086 * Tnr - 0.25 / std::pow(Hk - 1., 1.5)) / theta;
if (Ax2 < 0)
{
Ax2 = 0;
diff --git a/dart/src/wViscFluxes.h b/dart/src/wViscFluxes.h
index e6881f7..c295c98 100644
--- a/dart/src/wViscFluxes.h
+++ b/dart/src/wViscFluxes.h
@@ -13,9 +13,9 @@ namespace dart
class DART_API ViscFluxes
{
private:
- double Re;
public:
+ double Re;
ViscFluxes(double _Re);
~ViscFluxes();
Eigen::VectorXd ComputeFluxes(size_t iPoint, BLRegion *bl);
diff --git a/dart/src/wViscMesh.cpp b/dart/src/wViscMesh.cpp
index 0a36a37..aafc154 100644
--- a/dart/src/wViscMesh.cpp
+++ b/dart/src/wViscMesh.cpp
@@ -21,18 +21,21 @@ using namespace dart;
/* ----------------------------------------------------------------------------------- */
-ViscMesh::ViscMesh(){prevnMarkers = 0;}
+ViscMesh::ViscMesh(){nMarkers = 0; nElm = 0; prevnMarkers = 0;}
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())
+ if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size())
{
- throw std::runtime_error("dart::ViscMesh Wrong mesh size.");
+ throw std::runtime_error("dart::ViscMesh Wrong mesh sizes.");
}
nMarkers = _x.size();
+ x.resize(nMarkers, 0.);
+ y.resize(nMarkers, 0.);
+ z.resize(nMarkers, 0.);
nElm = nMarkers - 1;
x = _x;
y = _y;
@@ -44,12 +47,12 @@ void ViscMesh::SetGlob(std::vector<double> _x, std::vector<double> _y, std::vect
void ViscMesh::ComputeBLcoord()
{
Loc.resize(nMarkers, 0.);
- std::vector<double> dx(nElm,0.);
alpha.resize(nElm, 0.);
+ dx.resize(nElm, 0.);
for (size_t iPoint=0; iPoint<nElm; ++iPoint)
{
- alpha[iPoint] = std::atan2(y[iPoint+1]-y[iPoint], x[iPoint+1]-x[iPoint]-1);
+ alpha[iPoint] = std::atan2(y[iPoint+1]-y[iPoint], x[iPoint+1]-x[iPoint]);
/* dx = sqrt(Delta x ^2 + Delta y ^2) */
dx[iPoint] = std::sqrt((x[iPoint+1]-x[iPoint])*(x[iPoint+1]-x[iPoint]) + (y[iPoint+1]-y[iPoint])*(y[iPoint+1]-y[iPoint]));
Loc[iPoint + 1] = Loc[iPoint] + dx[iPoint];
diff --git a/dart/src/wViscMesh.h b/dart/src/wViscMesh.h
index 21d481f..ff1a53e 100644
--- a/dart/src/wViscMesh.h
+++ b/dart/src/wViscMesh.h
@@ -36,11 +36,12 @@ class DART_API ViscMesh
double Gety(size_t iPoint) const {return y[iPoint];};
double Getz(size_t iPoint) const {return z[iPoint];};
double GetExt(size_t iPoint) const {return Ext[iPoint];};
+ double GetAlpha(size_t iElm) const {return alpha[iElm];};
int GetDiffnElm() const {return nMarkers-prevnMarkers;};
void SetGlob(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
- void SetExt(std::vector<double> ExtM) {Ext = ExtM;};
+ void SetExt(std::vector<double> ExtM) {Ext.resize(ExtM.size(), 0.); Ext = ExtM;};
void ComputeBLcoord();
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 7297dfc..6e06caa 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -90,7 +90,7 @@ ViscSolver::~ViscSolver()
}
delete tSolver;
- PrintBanner();
+ //PrintBanner();
std::cout << "~ViscSolver()\n";
} // end ~ViscSolver
@@ -101,12 +101,12 @@ ViscSolver::~ViscSolver()
*/
int ViscSolver::Run(unsigned int couplIter)
{
- bool lockTrans;
- int pointExitCode;
+ bool lockTrans=false;
+ int pointExitCode = 0;
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
{
- std::cout << "\n - Sec " << iSec << " stagPtMvmt " << stagPtMvmt[iSec] << Sections[iSec][0]->mesh->GetDiffnElm();
+ //std::cout << "\n - Sec " << iSec << " stagPtMvmt " << stagPtMvmt[iSec] << Sections[iSec][0]->mesh->GetDiffnElm();
/* Check for stagnation point movement */
@@ -153,15 +153,15 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->SetitFrozenJac(1);
tSolver->SetinitSln(true);
- std::vector<double> xxExtCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
- for (size_t iPoint=0; iPoint<xxExtCurr.size(); ++iPoint)
+ if (Sections[iSec][iRegion]->mesh->GetDiffnElm())
{
- xxExtCurr[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
- }
- Sections[iSec][iRegion]->mesh->SetExt(xxExtCurr);
+ std::vector<double> xxExtCurr(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.);
+ for (size_t iPoint=0; iPoint<xxExtCurr.size(); ++iPoint)
+ {
+ xxExtCurr[iPoint] = Sections[iSec][iRegion]->mesh->GetLoc(iPoint);
+ }
+ Sections[iSec][iRegion]->mesh->SetExt(xxExtCurr);
- if (couplIter == 0)
- {
std::vector<double> DeltaStarZeros(Sections[iSec][iRegion]->mesh->GetnMarkers(), 0.0);
Sections[iSec][iRegion]->invBnd->SetDeltaStar(DeltaStarZeros);
@@ -205,10 +205,13 @@ int ViscSolver::Run(unsigned int couplIter)
}
/* Impose wake boundary condition */
- std::vector<double> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0);
- std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0);
- 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 */
+ std::vector<double> UpperCond(Sections[iSec][iRegion]->GetnVar(), 0.);
+ std::vector<double> LowerCond(Sections[iSec][iRegion]->GetnVar(), 0.);
+ for (size_t k=0; k<UpperCond.size(); ++k)
+ {
+ UpperCond[k] = Sections[iSec][0]->U[(Sections[iSec][0]->mesh->GetnMarkers() - 1) * Sections[iSec][0]->GetnVar() + k];
+ LowerCond[k] = Sections[iSec][0]->U[(Sections[iSec][1]->mesh->GetnMarkers() - 1) * Sections[iSec][1]->GetnVar() + k];
+ }
Sections[iSec][iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
lockTrans = true;
}
@@ -221,42 +224,15 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
- /* Time integration */
- /* if (iRegion == 0 && iPoint == 1)
- {
- std::cout << " iSec " << iSec
- << " iPoint " << iPoint
- << " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
- << " Vt " << Sections[iSec][iRegion]->invBnd->GetVt(iPoint)
- << " Me " << Sections[iSec][iRegion]->invBnd->GetMe(iPoint)
- << " xxExt " << Sections[iSec][iRegion]->mesh->GetExt(iPoint)
- << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
- << std::endl;
- Sections[iSec][iRegion]->PrintSolution(iPoint);
- } */
-
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
- /* if (iRegion == 1 && iPoint == 265)
- {
- std::cout << " iSec " << iSec
- << " iPoint " << iPoint
- << " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
- << " Vt " << Sections[iSec][iRegion]->invBnd->GetVt(iPoint)
- << " Me " << Sections[iSec][iRegion]->invBnd->GetMe(iPoint)
- << " xxExt " << Sections[iSec][iRegion]->mesh->GetExt(iPoint)
- << " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
- << std::scientific;
- Sections[iSec][iRegion]->PrintSolution(iPoint);
- } */
-
/* Unsucessfull convergence output */
if (pointExitCode != 0)
{
if (printOn){std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;}
solverExitCode = -1; // Convergence failed
- std::cout << " NC " << iPoint << " ";
+ //std::cout << " NC " << iPoint << " ";
}
@@ -303,11 +279,11 @@ void ViscSolver::CheckWaves(size_t iPoint, BLRegion *bl)
{
/* 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 * std::pow(1 / (bl->Hk[iPoint] - 1), 0.43) + 0.7 * (std::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)
+ if (std::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.
}
@@ -334,9 +310,9 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
bl->transMarker = iPoint; /* Save transition marker */
/* Loop over remaining points in the region */
- for (size_t k = 0; k < bl->mesh->GetnMarkers() - iPoint; ++k)
+ for (size_t k = iPoint; k < bl->mesh->GetnMarkers(); ++k)
{
- bl->Regime[iPoint + k] = 1; // Set Turbulent regime.
+ bl->Regime[k] = 1; // Set Turbulent regime.
}
/* Compute transition location */
@@ -345,14 +321,14 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
/* Percentage of the subsection corresponding to a turbulent flow */
double avTurb = (bl->mesh->Getx(iPoint) - bl->xtr) / (bl->mesh->Getx(iPoint) - bl->mesh->Getx(iPoint - 1));
- double avLam = 1 - avTurb;
+ 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->invBnd->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) */
+ /* Avoid starting with ill conditioned IC for turbulent computation. (Regime was switched above) */
/* 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.
@@ -368,7 +344,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
/* 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 + 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 */
@@ -451,45 +427,59 @@ void ViscSolver::SetGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x
void ViscSolver::SetVelocities(size_t iSec, size_t iRegion, std::vector<double> vx, std::vector<double> vy, std::vector<double> vz)
{
+ if (vx.size() != vy.size() || vy.size() != vz.size() || vx.size() != vz.size())
+ {
+ throw std::runtime_error("dart::ViscSolver Wrong velocity vector sizes.");
+ }
+ if (vx.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver Velocity vector is not coherent with the mesh.");
+ }
+
std::vector<double> Vt(vx.size(), 0.0);
- double alpha=0.0;
for (size_t iPoint=0; iPoint<vx.size()-1; ++iPoint)
{
- alpha = std::atan2((Sections[iSec][iRegion]->mesh->Gety(iPoint+1) - Sections[iSec][iRegion]->mesh->Gety(iPoint)),
- (Sections[iSec][iRegion]->mesh->Getx(iPoint+1) - Sections[iSec][iRegion]->mesh->Getx(iPoint)));
-
- Vt[iPoint] = vx[iPoint]*std::cos(alpha) + vy[iPoint]*std::sin(alpha);
- Vt[iPoint+1] = vx[iPoint+1]*std::cos(alpha) + vy[iPoint+1]*std::sin(alpha);
- /*if (Vt[iPoint] < 0)
- {
- std::cout << "Negative velocity at point " << iPoint << std::endl;
- //throw;
- }*/
- // std::cout << iSec << iRegion << iPoint << "vx " << vx[iPoint] << "vy" << vy[iPoint] << "vz" << vz[iPoint] << "vt" << Vt[iPoint] << std::endl;
+ Vt[iPoint] = vx[iPoint]*std::cos(Sections[iSec][iRegion]->mesh->GetAlpha(iPoint)) + vy[iPoint]*std::sin(Sections[iSec][iRegion]->mesh->GetAlpha(iPoint));
+ Vt[iPoint+1] = vx[iPoint+1]*std::cos(Sections[iSec][iRegion]->mesh->GetAlpha(iPoint)) + vy[iPoint+1]*std::sin(Sections[iSec][iRegion]->mesh->GetAlpha(iPoint));
}
- //Vt.back() = vx.back()*cos(alpha) + vy.back()*std::sin(alpha);
Sections[iSec][iRegion]->invBnd->SetVt(Vt);
}
void ViscSolver::SetMe(size_t iSec, size_t iRegion, std::vector<double> _Me)
{
+ if (_Me.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver Mach number vector is not coherent with the mesh.");
+ }
Sections[iSec][iRegion]->invBnd->SetMe(_Me);
}
void ViscSolver::SetRhoe(size_t iSec, size_t iRegion, std::vector<double> _Rhoe)
{
+ if (_Rhoe.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver Density vector is not coherent with the mesh.");
+ }
Sections[iSec][iRegion]->invBnd->SetRhoe(_Rhoe);
}
void ViscSolver::SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt)
{
+ if (_xxExt.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver External mesh vector is not coherent with the mesh.");
+ }
Sections[iSec][iRegion]->mesh->SetExt(_xxExt);
}
void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt)
{
+ if (_DeltaStarExt.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver External delta star vector is not coherent with the mesh.");
+ }
Sections[iSec][iRegion]->invBnd->SetDeltaStar(_DeltaStarExt);
}
@@ -533,6 +523,10 @@ std::vector<double> ViscSolver::ExtractUe(size_t iSec, size_t iRegion) const
void ViscSolver::SetUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
{
+ if (_ue.size() != Sections[iSec][iRegion]->mesh->GetnMarkers())
+ {
+ throw std::runtime_error("dart::ViscSolver External velocity vector is not coherent with the mesh.");
+ }
Sections[iSec][iRegion]->invBnd->SetVtOld(_ue);
}
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index eafbb9a..85f3292 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -21,7 +21,7 @@ private:
double Re,
CFL0;
- int solverExitCode = 0;
+ int solverExitCode;
std::vector<bool> stagPtMvmt;
diff --git a/dart/tests/lift.py b/dart/tests/lift.py
index 623c495..9f598bb 100644
--- a/dart/tests/lift.py
+++ b/dart/tests/lift.py
@@ -32,6 +32,7 @@ import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+from matplotlib import pyplot as plt
def main():
# timer
@@ -39,7 +40,7 @@ def main():
tms['total'].start()
# define flow variables
- alpha = 2*math.pi/180
+ alpha = 12*math.pi/180
M_inf = 0.
c_ref = 1
dim = 2
@@ -47,7 +48,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.0075, 'msLe' : 0.0075}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.001, 'msLe' : 0.001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -60,7 +61,16 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
- solver.run()
+ solver.printOn = False
+ for i in range(10):
+ solver.run()
+ print('{:.15f}'.format(solver.Cl))
+ Ue = floU.extract(bnd.groups[0].tag.elems, solver.U)
+ solver.reset()
+ #plt.plot(Ue[:,0], Ue[:,3])
+ plt.plot(Ue[:,0], Ue[:,4], 'x')
+ plt.show()
+
solver.save(gmshWriter)
tms['solver'].stop()
@@ -101,7 +111,6 @@ def main():
else:
raise Exception('Test not defined for this flow')
tests.run()"""
- print('{:.14f}'.format(solver.Cl))
# eof
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index 2d91256..036d604 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -81,9 +81,9 @@ class GroupConstructor():
cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
- dataUpper = np.empty((len(fMeshUp), 0))
- dataLower = np.empty((len(fMeshLw), 0))
- dataWake = np.empty((len(fMeshWk), 0))
+ dataUpper = np.zeros((len(fMeshUp), 0))
+ dataLower = np.zeros((len(fMeshLw), 0))
+ dataWake = np.zeros((len(fMeshWk), 0))
for iParam in range(len(dataIn[0][0][0,:])):
interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
dataUpper = np.column_stack((dataUpper, interpParamUp))
@@ -154,7 +154,7 @@ class GroupConstructor():
def interpolateToFineMesh(self, cVector, fMesh, nFactor):
- fVector = np.empty(len(fMesh)-1)
+ fVector = np.zeros(len(fMesh)-1)
for cPoint in range(len(cVector) - 1):
dv = cVector[cPoint + 1] - cVector[cPoint]
diff --git a/dart/viscous/Drivers/DGroupConstructor2.py b/dart/viscous/Drivers/DGroupConstructor2.py
index 44877c7..57cf76a 100755
--- a/dart/viscous/Drivers/DGroupConstructor2.py
+++ b/dart/viscous/Drivers/DGroupConstructor2.py
@@ -97,9 +97,9 @@ class GroupConstructor():
cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
- dataUpper = np.empty((len(fMeshUp), 0))
- dataLower = np.empty((len(fMeshLw), 0))
- dataWake = np.empty((len(fMeshWk), 0))
+ dataUpper = np.zeros((len(fMeshUp), 0))
+ dataLower = np.zeros((len(fMeshLw), 0))
+ dataWake = np.zeros((len(fMeshWk), 0))
for iParam in range(len(dataIn[0][0][0,:])):
interpParamUp = self.interpolateToFineMesh(dataIn[0][0][:,iParam], fMeshUp, nFactor)
dataUpper = np.column_stack((dataUpper, interpParamUp))
@@ -164,7 +164,7 @@ class GroupConstructor():
def interpolateToFineMesh(self, cVector, fMesh, nFactor):
- fVector = np.empty(len(fMesh)-1)
+ fVector = np.zeros(len(fMesh)-1)
for cPoint in range(len(cVector) - 1):
dv = cVector[cPoint + 1] - cVector[cPoint]
diff --git a/dart/viscous/Physics/DDataStructures.py b/dart/viscous/Physics/DDataStructures.py
index ca867fb..1404b36 100755
--- a/dart/viscous/Physics/DDataStructures.py
+++ b/dart/viscous/Physics/DDataStructures.py
@@ -80,8 +80,8 @@ class PSolution:
def initializeSolution(self, PGeometry):
nMarkers = PGeometry.GetnMarkers()
- self.u = np.empty(self.__nVar * nMarkers)
- self.blPar = np.empty(self.__nVar * nMarkers)
+ self.u = np.zeros(self.__nVar * nMarkers)
+ self.blPar = np.zeros(self.__nVar * nMarkers)
return self.u, self.blPar
--
GitLab