Hybrid parallel boundary conditions for turbulence and compressible flow solvers#975
Hybrid parallel boundary conditions for turbulence and compressible flow solvers#975pcarruscag merged 42 commits intodevelopfrom
Conversation
…d changing the signature of all functions)
| for (iVar = 0; iVar < nVar; iVar++) | ||
| Solution[iVar] = 0.5* (Solution_i[iVar] + Solution_j[iVar]); |
There was a problem hiding this comment.
This function (CSolver::SetRotatingFrame_GCL) was averaging the solution but then updating the iPoint / jPoint residuals with Solution_i / j instead, if anyone knows of a good reason for that let me know and I'll put it back, otherwise "it's fixed".
There was a problem hiding this comment.
Actually this apparent fix makes one of the turbomachinery cases fail, so I've put the apparent bug back.
There was a problem hiding this comment.
Can you make an issue for this please with a link to the spot in the code? I think this could be a bug too but need to look through it...
…it and CHT interface
| */ | ||
| void SetSpline(vector<su2double> &x, vector<su2double> &y, unsigned long n, su2double yp1, su2double ypn, vector<su2double> &y2); | ||
|
|
||
| /*! | ||
| * \brief Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai’s in order), | ||
| and given the array y2a[1..n], which is the output from spline above, and given a value of | ||
| x, this routine returns a cubic-spline interpolated value y. | ||
| Numerical Recipes: The Art of Scientific Computing, Third Edition in C++. | ||
| * \return The interpolated value of for x. | ||
| */ | ||
| su2double GetSpline(vector<su2double> &xa, vector<su2double> &ya, vector<su2double> &y2a, unsigned long n, su2double x); |
There was a problem hiding this comment.
This code was unused and duplicated in CGeometry.
| SU2_OMP_MASTER | ||
| { | ||
| CNumerics* conv_bound_numerics = numerics[CONV_BOUND_TERM + omp_get_thread_num()*MAX_TERMS]; | ||
| CNumerics* visc_bound_numerics = numerics[VISC_BOUND_TERM + omp_get_thread_num()*MAX_TERMS]; |
There was a problem hiding this comment.
Usual change no 1: One numerics object per thread.
| Residual_i = new su2double[nVar](); | ||
| Residual_j = new su2double[nVar](); | ||
| Res_Conv = new su2double[nVar](); | ||
| Res_Visc = new su2double[nVar](); | ||
| Res_Sour = new su2double[nVar](); |
There was a problem hiding this comment.
Usual change no 2: Avoid using the small arrays of CSolver.
| SU2_OMP_FOR_DYN(OMP_MIN_SIZE) | ||
| for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { |
There was a problem hiding this comment.
Usual change no 3: Add worksharing directives to for loops.
| AddDynamicGridResidualContribution(iPoint, Point_Normal, geometry, UnitNormal, | ||
| Area, geometry->nodes->GetGridVel(iPoint), | ||
| Jacobian_i, Res_Conv, Res_Visc); |
There was a problem hiding this comment.
Significant change no 1: Find duplicated code and move it to something that can be reused.
| void CNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, | ||
| CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { | ||
|
|
||
| LinSysRes.SubtractBlock(iPoint, Res_Visc); | ||
| BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); | ||
| } | ||
|
|
||
| /*--- Enforce the no-slip boundary condition in a strong way by | ||
| modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ | ||
| void CNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, | ||
| CConfig *config, unsigned short val_marker) { | ||
|
|
||
| if (implicit) { | ||
| for (iVar = 1; iVar <= nDim; iVar++) { | ||
| total_index = iPoint*nVar+iVar; | ||
| Jacobian.DeleteValsRowi(total_index); | ||
| } | ||
| } | ||
| } | ||
| } | ||
| BC_Isothermal_Wall_Generic(geometry, solver_container, conv_numerics, nullptr, config, val_marker, true); | ||
| } |
There was a problem hiding this comment.
Significant change no 2: The Isothermal and ConjugateHeat BC's are basically the same thing but the latter gets its wall temperature from a different place, so I moved the difference to a small private function, half the code, half the bugs.
| su2double Vel[3] = {0.0, 0.0, 0.0}, VelNormal, VelTang[3], VelTangMod, VelInfMod, WallDist[3], WallDistMod; | ||
| su2double T_Normal, P_Normal; | ||
| su2double Density_Wall, T_Wall, P_Wall, Lam_Visc_Wall, Tau_Wall = 0.0, Tau_Wall_Old = 0.0; | ||
| su2double *Coord, *Coord_Normal; | ||
| su2double diff, Delta; | ||
| su2double U_Tau, U_Plus, Gam, Beta, Phi, Q, Y_Plus_White, Y_Plus; | ||
| su2double TauElem[3], TauNormal, TauTangent[3], WallShearStress; | ||
| su2double Gas_Constant = config->GetGas_ConstantND(); | ||
| su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; |
There was a problem hiding this comment.
Usual change no 5: Define everything as it is used.
| void CTurbSSTSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, | ||
| CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { | ||
|
|
||
| unsigned long iPoint, jPoint, iVertex, total_index; | ||
| unsigned short iDim, iVar; | ||
| su2double distance, density = 0.0, laminar_viscosity = 0.0, beta_1; | ||
|
|
||
| for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { | ||
| iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); | ||
|
|
||
| /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ | ||
| if (geometry->nodes->GetDomain(iPoint)) { | ||
|
|
||
| /*--- distance to closest neighbor ---*/ | ||
| jPoint = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); | ||
| distance = 0.0; | ||
| for (iDim = 0; iDim < nDim; iDim++) { | ||
| distance += (geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim))* | ||
| (geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim)); | ||
| } | ||
| distance = sqrt(distance); | ||
|
|
||
| /*--- Set wall values ---*/ | ||
|
|
||
| density = solver_container[FLOW_SOL]->GetNodes()->GetDensity(jPoint); | ||
| laminar_viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(jPoint); | ||
|
|
||
| beta_1 = constants[4]; | ||
|
|
||
| Solution[0] = 0.0; | ||
| Solution[1] = 60.0*laminar_viscosity/(density*beta_1*distance*distance); | ||
|
|
||
| /*--- Set the solution values and zero the residual ---*/ | ||
| nodes->SetSolution_Old(iPoint,Solution); | ||
| nodes->SetSolution(iPoint,Solution); | ||
| LinSysRes.SetBlock_Zero(iPoint); | ||
|
|
||
| /*--- Change rows of the Jacobian (includes 1 in the diagonal) ---*/ | ||
| for (iVar = 0; iVar < nVar; iVar++) { | ||
| total_index = iPoint*nVar+iVar; | ||
| Jacobian.DeleteValsRowi(total_index); | ||
| } | ||
|
|
||
| } | ||
| } | ||
| BC_HeatFlux_Wall(geometry, solver_container, conv_numerics, visc_numerics, config, val_marker); | ||
|
|
||
| } |
|
@talbring files have been moved |
talbring
left a comment
There was a problem hiding this comment.
Very nice. You also did a good clean-up of some stuff :) Thanks!
Proposed Changes
So far BC's were done by a single thread, this PR addresses that which required making all that code thread-safe. There is also some duplication-avoiding cleanup.
A batch of regression tests for hybrid parallel (compilation + execution) was added. Boundary conditions are perhaps the easiest places to introduce non-thread-safe code, due to the common pattern of writing into temporary arrays of CSolver (e.g. Solution_i), thus it becomes important to have some tests in place.
Related Work
Continuation of #789, so far as compressible flow goes we now have "full coverage", the only major areas left out (that matter to scalability) are MPI and parts of the initial preprocessing.
PR Checklist