Hybrid parallel (OpenMP) support for incompressible solvers by pcarruscag · Pull Request #1178 · su2code/SU2
/*--- Read the restart data from either an ASCII or binary SU2 file. ---*/
if (config->GetRead_Binary_Restart()) { Read_SU2_Restart_Binary(geometry[MESH_0], config, restart_filename); } else { Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); }
/*--- Load data from the restart into correct containers. ---*/
counter = 0; for (iPoint_Global = 0; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++ ) {
/*--- Retrieve local index. If this node from the restart file lives on the current processor, we will load and instantiate the vars. ---*/
iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
if (iPoint_Local > -1) {
/*--- We need to store this point's data, so jump to the correct offset in the buffer of data from the restart file and load it. ---*/
index = counter*Restart_Vars[1] + skipVars; for (iVar = 0; iVar < nVar_Restart; iVar++) Solution[iVar] = Restart_Data[index+iVar]; nodes->SetSolution(iPoint_Local,Solution); iPoint_Global_Local++;
/*--- For dynamic meshes, read in and store the grid coordinates and grid velocities for each node. ---*/
if (dynamic_grid && val_update_geo) {
/*--- Read in the next 2 or 3 variables which are the grid velocities ---*/ /*--- If we are restarting the solution from a previously computed static calculation (no grid movement) ---*/ /*--- the grid velocities are set to 0. This is useful for FSI computations ---*/
/*--- Rewind the index to retrieve the Coords. ---*/ index = counter*Restart_Vars[1]; for (iDim = 0; iDim < nDim; iDim++) { Coord[iDim] = Restart_Data[index+iDim]; }
su2double GridVel[MAXNDIM] = {0.0}; if (!steady_restart) { /*--- Move the index forward to get the grid velocities. ---*/ index = counter*Restart_Vars[1] + skipVars + nVar_Restart + turbVars; for (iDim = 0; iDim < nDim; iDim++) { GridVel[iDim] = Restart_Data[index+iDim]; } }
for (iDim = 0; iDim < nDim; iDim++) { geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]); } }
/*--- For static FSI problems, grid_movement is 0 but we need to read in and store the grid coordinates for each node (but not the grid velocities, as there are none). ---*/
if (static_fsi && val_update_geo) { /*--- Rewind the index to retrieve the Coords. ---*/ index = counter*Restart_Vars[1]; for (iDim = 0; iDim < nDim; iDim++) { Coord[iDim] = Restart_Data[index+iDim];}
for (iDim = 0; iDim < nDim; iDim++) { geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]); } }
/*--- Increment the overall counter for how many points have been loaded. ---*/ counter++;
} }
/*--- Detect a wrong solution file ---*/
if (iPoint_Global_Local < nPointDomain) { SU2_MPI::Error(string("The solution file ") + restart_filename + string(" doesn't match with the mesh file!\n") + string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); }
/*--- Update the geometry for flows on deforming meshes ---*/
if ((dynamic_grid || static_fsi) && val_update_geo) {
/*--- Communicate the new coordinates and grid velocities at the halos ---*/
geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES); geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES);
if (dynamic_grid) { geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); }
/*--- Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/
geometry[MESH_0]->SetControlVolume(config, UPDATE); geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry[MESH_0]->SetMaxLength(config);
/*--- Update the multigrid structure after setting up the finest grid, including computing the grid velocities on the coarser levels. ---*/
for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { iMeshFine = iMesh-1; geometry[iMesh]->SetControlVolume(config, geometry[iMeshFine], UPDATE); geometry[iMesh]->SetBoundControlVolume(config, geometry[iMeshFine],UPDATE); geometry[iMesh]->SetCoord(geometry[iMeshFine]); if (dynamic_grid) { geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMeshFine], config); } geometry[iMesh]->SetMaxLength(config); } }
/*--- Communicate the loaded solution on the fine grid before we transfer it down to the coarse levels. We alo call the preprocessing routine on the fine level in order to have all necessary quantities updated, especially if this is a turbulent simulation (eddy viscosity). ---*/
solver[MESH_0][FLOW_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); solver[MESH_0][FLOW_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION);
/*--- For turbulent simulations the flow preprocessing is done by the turbulence solver * after it loads its variables (they are needed to compute flow primitives). ---*/ if (!turbulent) { solver[MESH_0][FLOW_SOL]->Preprocessing(geometry[MESH_0], solver[MESH_0], config, MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, false); }
/*--- Interpolate the solution down to the coarse multigrid levels ---*/
for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { for (iPoint = 0; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) { Area_Parent = geometry[iMesh]->nodes->GetVolume(iPoint); for (iVar = 0; iVar < nVar; iVar++) Solution[iVar] = 0.0; for (iChildren = 0; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) { Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren); Area_Children = geometry[iMesh-1]->nodes->GetVolume(Point_Fine); Solution_Fine = solver[iMesh-1][FLOW_SOL]->GetNodes()->GetSolution(Point_Fine); for (iVar = 0; iVar < nVar; iVar++) { Solution[iVar] += Solution_Fine[iVar]*Area_Children/Area_Parent; } } solver[iMesh][FLOW_SOL]->GetNodes()->SetSolution(iPoint,Solution); } solver[iMesh][FLOW_SOL]->InitiateComms(geometry[iMesh], config, SOLUTION); solver[iMesh][FLOW_SOL]->CompleteComms(geometry[iMesh], config, SOLUTION);
if (!turbulent) { solver[iMesh][FLOW_SOL]->Preprocessing(geometry[iMesh], solver[iMesh], config, iMesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); } }
/*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy ---*/ if (dual_time && config->GetGrid_Movement() && !config->GetDeform_Mesh() && (config->GetKind_GridMovement() != RIGID_MOTION)) { Restart_OldGeometry(geometry[MESH_0], config); }
/*--- Delete the class memory that is used to load the restart. ---*/
delete [] Restart_Vars; Restart_Vars = nullptr; delete [] Restart_Data; Restart_Data = nullptr; LoadRestart_impl(geometry, solver, config, val_iter, val_update_geo, Solution, nVar_Restart);