diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 8a4fb03d44511956c40cce0a98e8bfcc4cfca715..050c6dc0610fcbf86dbdd3b378af857d8476202d 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -257,13 +257,9 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
for (size_t k = iPoint; k < bl->nNodes; ++k)
bl->regime[k] = 1; // Set Turbulent regime.
- // Compute transition location.
- // if (forced)
- // bl->xtr = bl->xtrF;
- // else
- bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->x[iPoint] - bl->x[iPoint-1]) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->x[iPoint-1];
+ bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->xoc[iPoint] - bl->xoc[iPoint-1]) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->xoc[iPoint-1];
- double avTurb = (bl->x[iPoint] - bl->xtr) / (bl->x[iPoint] - bl->x[iPoint-1]);
+ double avTurb = (bl->xoc[iPoint] - bl->xtr) / (bl->xoc[iPoint] - bl->xoc[iPoint-1]);
double avLam = 1. - avTurb;
// std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
@@ -293,16 +289,19 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
// Solve point in turbulent configuration.
int exitCode = tSolver->integration(iPoint, bl);
- if (exitCode != 0 && verbose > 1)
+ if (exitCode != 0)
std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
// Average both solutions (Now solution @ iPoint is the turbulent solution).
bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
- bl->u[iPoint * nVar + 2] = 9.; // N.
+ bl->u[iPoint * nVar + 2] = bl->getnCrit(); // 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.
+ bl->deltaStar[iPoint-1] = bl->u[(iPoint-1)*nVar + 0] * bl->u[(iPoint-1)*nVar + 1];
+ bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1];
+
// Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
std::vector<double> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->Me[iPoint-1], bl->name);
@@ -328,7 +327,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
}
/**
- * @brief Friction, pressure and total drag computation qt each station of the configuration.
+ * @brief Friction, pressure and total drag computation at each station of the configuration.
*
* @tparam Cdt = theta * Ue^((H+5)/2).
* @tparam Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
@@ -357,11 +356,9 @@ void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, size_
double Cdf_lower = 0.0;
for (size_t k = 1; k < sec[1]->nNodes; ++k)
Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (sec[1]->x[k] - sec[1]->x[k-1]) / 2;
-
- // Friction drag coefficient.
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
- // Total drag coefficient. Compute from upper and lower TE if there is no wake.
+ // Total drag coefficient. Computed from upper and lower TE if there is no wake.
if (sec.size() > 2 && sec[2]->getName() == "wake")
{
size_t lastWkPt = (sec[2]->nNodes - 1) * nVar;
diff --git a/blast/tests/dart/naca0012_2D.py b/blast/tests/dart/naca0012_2D.py
index 04ab7904d42a3722cdce12bd87553759d8099649..a5cd8f1a9849fad2ea56d07ee15a968dbde53124 100644
--- a/blast/tests/dart/naca0012_2D.py
+++ b/blast/tests/dart/naca0012_2D.py
@@ -143,6 +143,7 @@ def main():
tests.add(CTest('Cdf', vsol.Cdf, 0.0047, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
tests.add(CTest('xtr_top', vsol.getAverageTransition(0), 0.196, 0.03, forceabs=True)) # XFOIL 0.1877
tests.add(CTest('xtr_bot', vsol.getAverageTransition(1), 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+ tests.add(CTest('error dStar', np.linalg.norm(vSolution[0]['deltaStar'] - vSolution[0]['theta']*vSolution[0]['H']), 0., 0., forceabs=True))
tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 19, 0, forceabs=True))
tests.run()
diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
index ef462b89c1a17f7710be828104fba0cd4dd64599..19a1e1e9dffa1abdebe11f595d104afbb8a2ca83 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/raeValidation.py
@@ -29,6 +29,7 @@ import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import os.path
+import tbox
import math
@@ -42,7 +43,7 @@ def cfgInviscid(nthrds, verb):
'Verb' : verb, # verbosity
# Model (geometry or mesh)
'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/rae_2.geo', # Input file containing the model
- 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.0075, 'msLe' : 0.001}, # parameters for input file model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -64,13 +65,13 @@ def cfgInviscid(nthrds, verb):
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
- 'LSolver' : 'SparseLu', # inner solver (Pardiso, MUMPS or GMRES)
+ 'LSolver' : 'PARDISO', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
'Rel_tol' : 1e-6, # relative tolerance on solver residual
'Abs_tol' : 1e-8, # absolute tolerance on solver residual
- 'Max_it' : 50 # solver maximum number of iterations
+ 'Max_it' : 100 # solver maximum number of iterations
}
def cfgBlast(verb):
@@ -84,7 +85,7 @@ def cfgBlast(verb):
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0],
- 'xtrF' : [0., 0.],# Forced transition location
+ 'xtrF' : [0.03, 0.03],# Forced transition location
'interpolator' : 'Matching',
'nDim' : 2
}
@@ -98,6 +99,12 @@ def main():
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
+ try:
+ _ = tbox.Pardiso()
+ except:
+ print('PARDISO not found, using SparseLu')
+ icfg['LSolver'] = 'SparseLu'
+
gr = 1.1
if gr == 1.0:
icfg['Pars']['msF'] = icfg['Pars']['msLe']
@@ -132,12 +139,12 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', isol.getCl(), 0.730, 5e-2))
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0095, 1e-3, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.740, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0096, 1e-3, forceabs=True))
tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0132, 1e-3, forceabs=True))
- tests.add(CTest('Cdf', vsol.Cdf, 0.0068, 1e-3, forceabs=True))
+ tests.add(CTest('Cdf', vsol.Cdf, 0.0067, 1e-3, forceabs=True))
if icfg['LSolver'] == 'SparseLu':
- tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 21, 0, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 34, 0, forceabs=True))
tests.run()
expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')