Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hybrid parallel boundary conditions for turbulence and compressible flow solvers #975

Merged
merged 42 commits into from
May 28, 2020

Conversation

pcarruscag
Copy link
Member

@pcarruscag pcarruscag commented May 9, 2020

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

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags).
  • My contribution is commented and consistent with SU2 style.

Comment on lines -2747 to -2748
for (iVar = 0; iVar < nVar; iVar++)
Solution[iVar] = 0.5* (Solution_i[iVar] + Solution_j[iVar]);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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".

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually this apparent fix makes one of the turbomachinery cases fail, so I've put the apparent bug back.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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...

@pcarruscag pcarruscag changed the title Hybrid parallel boundary conditions for turbulence and compressible flow solvers [WIP] Hybrid parallel boundary conditions for turbulence and compressible flow solvers May 11, 2020
Comment on lines -8351 to -8361
*/
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);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code was unused and duplicated in CGeometry.

Comment on lines -81 to +82
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];
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Usual change no 1: One numerics object per thread.

Comment on lines -175 to -179
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]();
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Usual change no 2: Avoid using the small arrays of CSolver.

Comment on lines +7286 to 7287
SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) {
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Usual change no 3: Add worksharing directives to for loops.

Comment on lines +1156 to +1158
AddDynamicGridResidualContribution(iPoint, Point_Normal, geometry, UnitNormal,
Area, geometry->nodes->GetGridVel(iPoint),
Jacobian_i, Res_Conv, Res_Visc);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Significant change no 1: Find duplicated code and move it to something that can be reused.

Comment on lines +1385 to 1395
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);
}
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Comment on lines -1840 to -1848
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;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Usual change no 5: Define everything as it is used.

Comment on lines 432 to 437
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);

}
Copy link
Member Author

@pcarruscag pcarruscag May 13, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why...

@pcarruscag pcarruscag mentioned this pull request May 27, 2020
@pcarruscag
Copy link
Member Author

@talbring files have been moved

Copy link
Member

@talbring talbring left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice. You also did a good clean-up of some stuff :) Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants