(fix) Fixed bug in transition

deltaStar was not recomputed for the transition point after averaging solutions, leading to oscillations in the transition region. Test cases have been updated accordingly and deltaStar computation is now tested in the 2D subsonic test case

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;

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()
 

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')