@@ -1985,180 +1985,120 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes
19851985void CEulerSolver::SetInitialCondition (CGeometry **geometry, CSolver ***solver_container, CConfig *config, unsigned long TimeIter) {
19861986
19871987 const bool restart = (config->GetRestart () || config->GetRestart_Flow ());
1988- const bool rans = (config->GetKind_Turb_Model () != NONE);
1989- const bool dual_time = ((config->GetTime_Marching () == DT_STEPPING_1ST) ||
1990- (config->GetTime_Marching () == DT_STEPPING_2ND));
19911988 const bool SubsonicEngine = config->GetSubsonicEngine ();
19921989
1993- /* --- Start OpenMP parallel region. ---*/
1994-
1995- SU2_OMP_PARALLEL {
1996-
1997- unsigned long iPoint;
1998- unsigned short iMesh, iDim;
1999- su2double X0[MAXNDIM] = {0.0 }, X1[MAXNDIM] = {0.0 }, X2[MAXNDIM] = {0.0 },
2000- X1_X0[MAXNDIM] = {0.0 }, X2_X0[MAXNDIM] = {0.0 }, X2_X1[MAXNDIM] = {0.0 },
2001- CP[MAXNDIM] = {0.0 }, Distance, DotCheck, Radius;
2002-
2003- /* --- Check if a verification solution is to be computed. ---*/
2004- if ((VerificationSolution) && (TimeIter == 0 ) && !restart) {
1990+ /* --- Use default implementation, then add solver-specifics. ---*/
20051991
2006- /* --- Loop over the multigrid levels. ---*/
2007- for (iMesh = 0 ; iMesh <= config->GetnMGLevels (); iMesh++) {
1992+ BaseClass::SetInitialCondition (geometry, solver_container, config, TimeIter);
20081993
2009- /* --- Loop over all grid points. ---*/
2010- SU2_OMP_FOR_STAT (omp_chunk_size)
2011- for (iPoint = 0 ; iPoint < geometry[iMesh]->GetnPoint (); iPoint++) {
1994+ /* --- Set subsonic initial condition for engine intakes at iteration 0 ---*/
20121995
2013- /* Set the pointers to the coordinates and solution of this DOF. */
2014- const su2double *coor = geometry[iMesh]->nodes ->GetCoord (iPoint);
2015- su2double *solDOF = solver_container[iMesh][FLOW_SOL]->GetNodes ()->GetSolution (iPoint);
1996+ if (!SubsonicEngine || (TimeIter != 0 ) || restart) return ;
20161997
2017- /* Set the solution in this DOF to the initial condition provided by
2018- the verification solution class. This can be the exact solution,
2019- but this is not necessary. */
2020- VerificationSolution->GetInitialCondition (coor, solDOF);
2021-
2022- }
2023- }
2024- }
2025-
2026- /* --- Set subsonic initial condition for engine intakes ---*/
2027-
2028- if (SubsonicEngine) {
2029-
2030- /* --- Set initial boundary condition at iteration 0 ---*/
1998+ /* --- Start OpenMP parallel region. ---*/
20311999
2032- if ((TimeIter == 0 ) && (!restart)) {
2000+ SU2_OMP_PARALLEL {
20332001
2034- su2double Velocity_Cyl[3 ] = {0.0 , 0.0 , 0.0 }, Velocity_CylND[3 ] = {0.0 , 0.0 , 0.0 }, Viscosity_Cyl,
2035- Density_Cyl, Density_CylND, Pressure_CylND, ModVel_Cyl, ModVel_CylND, Energy_CylND,
2036- T_ref = 0.0 , S = 0.0 , Mu_ref = 0.0 ;
2037- const su2double *Coord, *SubsonicEngine_Cyl, *SubsonicEngine_Values;
2002+ unsigned long iPoint;
2003+ unsigned short iMesh, iDim;
2004+ su2double X0[MAXNDIM] = {0.0 }, X1[MAXNDIM] = {0.0 }, X2[MAXNDIM] = {0.0 },
2005+ X1_X0[MAXNDIM] = {0.0 }, X2_X0[MAXNDIM] = {0.0 }, X2_X1[MAXNDIM] = {0.0 },
2006+ CP[MAXNDIM] = {0.0 }, Distance, DotCheck, Radius;
20382007
2039- SubsonicEngine_Values = config->GetSubsonicEngine_Values ();
2040- su2double Mach_Cyl = SubsonicEngine_Values[0 ];
2041- su2double Alpha_Cyl = SubsonicEngine_Values[1 ];
2042- su2double Beta_Cyl = SubsonicEngine_Values[2 ];
2043- su2double Pressure_Cyl = SubsonicEngine_Values[3 ];
2044- su2double Temperature_Cyl = SubsonicEngine_Values[4 ];
2008+ su2double Velocity_Cyl[MAXNDIM] = {0.0 }, Velocity_CylND[MAXNDIM] = {0.0 }, Viscosity_Cyl,
2009+ Density_Cyl, Density_CylND, Pressure_CylND, ModVel_Cyl, ModVel_CylND, Energy_CylND,
2010+ T_ref = 0.0 , S = 0.0 , Mu_ref = 0.0 ;
2011+ const su2double *Coord, *SubsonicEngine_Cyl, *SubsonicEngine_Values;
20452012
2046- su2double Alpha = Alpha_Cyl*PI_NUMBER/180.0 ;
2047- su2double Beta = Beta_Cyl*PI_NUMBER/180.0 ;
2013+ SubsonicEngine_Values = config->GetSubsonicEngine_Values ();
2014+ su2double Mach_Cyl = SubsonicEngine_Values[0 ];
2015+ su2double Alpha_Cyl = SubsonicEngine_Values[1 ];
2016+ su2double Beta_Cyl = SubsonicEngine_Values[2 ];
2017+ su2double Pressure_Cyl = SubsonicEngine_Values[3 ];
2018+ su2double Temperature_Cyl = SubsonicEngine_Values[4 ];
20482019
2049- su2double Gamma_Minus_One = Gamma - 1 .0 ;
2050- su2double Gas_Constant = config-> GetGas_Constant () ;
2020+ su2double Alpha = Alpha_Cyl*PI_NUMBER/ 180 .0 ;
2021+ su2double Beta = Beta_Cyl*PI_NUMBER/ 180.0 ;
20512022
2052- su2double Mach2Vel_Cyl = sqrt (Gamma*Gas_Constant*Temperature_Cyl);
2023+ su2double Gamma_Minus_One = Gamma - 1.0 ;
2024+ su2double Gas_Constant = config->GetGas_Constant ();
20532025
2054- for (iMesh = 0 ; iMesh <= config-> GetnMGLevels (); iMesh++) {
2026+ su2double Mach2Vel_Cyl = sqrt (Gamma*Gas_Constant*Temperature_Cyl);
20552027
2056- SU2_OMP_FOR_STAT (omp_chunk_size)
2057- for (iPoint = 0 ; iPoint < geometry[iMesh]->GetnPoint (); iPoint++) {
2028+ for (iMesh = 0 ; iMesh <= config->GetnMGLevels (); iMesh++) {
20582029
2059- Velocity_Cyl[0 ] = cos (Alpha)*cos (Beta)*Mach_Cyl*Mach2Vel_Cyl;
2060- Velocity_Cyl[1 ] = sin (Beta)*Mach_Cyl*Mach2Vel_Cyl;
2061- Velocity_Cyl[2 ] = sin (Alpha)*cos (Beta)*Mach_Cyl*Mach2Vel_Cyl;
2030+ auto FlowNodes = solver_container[iMesh][FLOW_SOL]->GetNodes ();
20622031
2063- ModVel_Cyl = 0.0 ;
2064- for (iDim = 0 ; iDim < nDim; iDim++) {
2065- ModVel_Cyl += Velocity_Cyl[iDim]*Velocity_Cyl[iDim];
2066- }
2067- ModVel_Cyl = sqrt (ModVel_Cyl);
2068-
2069- if (config->GetViscous ()) {
2070- if (config->GetSystemMeasurements () == SI) { T_ref = 273.15 ; S = 110.4 ; Mu_ref = 1.716E-5 ; }
2071- if (config->GetSystemMeasurements () == US) {
2072- T_ref = (273.15 - 273.15 ) * 1.8 + 491.67 ;
2073- S = (110.4 - 273.15 ) * 1.8 + 491.67 ;
2074- Mu_ref = 1.716E-5 /47.88025898 ;
2075- }
2076- Viscosity_Cyl = Mu_ref*(pow (Temperature_Cyl/T_ref, 1.5 ) * (T_ref+S)/(Temperature_Cyl+S));
2077- Density_Cyl = config->GetReynolds ()*Viscosity_Cyl/(ModVel_Cyl*config->GetLength_Reynolds ());
2078- Pressure_Cyl = Density_Cyl*Gas_Constant*Temperature_Cyl;
2079- }
2080- else {
2081- Density_Cyl = Pressure_Cyl/(Gas_Constant*Temperature_Cyl);
2082- }
2032+ SU2_OMP_FOR_STAT (omp_chunk_size)
2033+ for (iPoint = 0 ; iPoint < geometry[iMesh]->GetnPoint (); iPoint++) {
20832034
2084- Density_CylND = Density_Cyl/config->GetDensity_Ref ();
2085- Pressure_CylND = Pressure_Cyl/config->GetPressure_Ref ();
2035+ Velocity_Cyl[0 ] = cos (Alpha)*cos (Beta)*Mach_Cyl*Mach2Vel_Cyl;
2036+ Velocity_Cyl[1 ] = sin (Beta)*Mach_Cyl*Mach2Vel_Cyl;
2037+ Velocity_Cyl[2 ] = sin (Alpha)*cos (Beta)*Mach_Cyl*Mach2Vel_Cyl;
20862038
2087- for (iDim = 0 ; iDim < nDim; iDim++) {
2088- Velocity_CylND[iDim] = Velocity_Cyl[iDim]/config->GetVelocity_Ref ();
2089- }
2039+ ModVel_Cyl = GeometryToolbox::Norm (nDim, Velocity_Cyl);
20902040
2091- ModVel_CylND = 0.0 ;
2092- for (iDim = 0 ; iDim < nDim; iDim++) {
2093- ModVel_CylND += Velocity_CylND[iDim]*Velocity_CylND[iDim];
2041+ if (config->GetViscous ()) {
2042+ if (config->GetSystemMeasurements () == SI) { T_ref = 273.15 ; S = 110.4 ; Mu_ref = 1.716E-5 ; }
2043+ if (config->GetSystemMeasurements () == US) {
2044+ T_ref = (273.15 - 273.15 ) * 1.8 + 491.67 ;
2045+ S = (110.4 - 273.15 ) * 1.8 + 491.67 ;
2046+ Mu_ref = 1.716E-5 /47.88025898 ;
20942047 }
2095- ModVel_CylND = sqrt (ModVel_CylND);
2096-
2097- Energy_CylND = Pressure_CylND/(Density_CylND*Gamma_Minus_One)+0.5 *ModVel_CylND*ModVel_CylND;
2048+ Viscosity_Cyl = Mu_ref*(pow (Temperature_Cyl/T_ref, 1.5 ) * (T_ref+S)/(Temperature_Cyl+S));
2049+ Density_Cyl = config->GetReynolds ()*Viscosity_Cyl/(ModVel_Cyl*config->GetLength_Reynolds ());
2050+ Pressure_Cyl = Density_Cyl*Gas_Constant*Temperature_Cyl;
2051+ }
2052+ else {
2053+ Density_Cyl = Pressure_Cyl/(Gas_Constant*Temperature_Cyl);
2054+ }
20982055
2099- Coord = geometry[iMesh]->nodes ->GetCoord (iPoint);
2056+ Density_CylND = Density_Cyl/config->GetDensity_Ref ();
2057+ Pressure_CylND = Pressure_Cyl/config->GetPressure_Ref ();
21002058
2101- SubsonicEngine_Cyl = config->GetSubsonicEngine_Cyl ();
2059+ for (iDim = 0 ; iDim < nDim; iDim++) {
2060+ Velocity_CylND[iDim] = Velocity_Cyl[iDim]/config->GetVelocity_Ref ();
2061+ }
21022062
2103- X0[0 ] = Coord[0 ]; X0[1 ] = Coord[1 ]; X0[2 ] = Coord[2 ];
2104- X1[0 ] = SubsonicEngine_Cyl[0 ]; X1[1 ] = SubsonicEngine_Cyl[1 ]; X1[2 ] = SubsonicEngine_Cyl[2 ];
2105- X2[0 ] = SubsonicEngine_Cyl[3 ]; X2[1 ] = SubsonicEngine_Cyl[4 ]; X2[2 ] = SubsonicEngine_Cyl[5 ];
2106- Radius = SubsonicEngine_Cyl[6 ];
2063+ ModVel_CylND = GeometryToolbox::Norm (nDim, Velocity_CylND);
21072064
2108- for (iDim = 0 ; iDim < nDim; iDim++) {
2109- X2_X1[iDim]= X1[iDim] - X2[iDim];
2110- X1_X0[iDim]= X0[iDim] - X1[iDim];
2111- X2_X0[iDim]= X0[iDim] - X2[iDim];
2112- }
2065+ Energy_CylND = Pressure_CylND/(Density_CylND*Gamma_Minus_One)+0.5 *ModVel_CylND*ModVel_CylND;
21132066
2114- CP[0 ] = (X2_X1[1 ]*X1_X0[2 ] - X2_X1[2 ]*X1_X0[1 ]);
2115- CP[1 ] = (X2_X1[2 ]*X1_X0[0 ] - X2_X1[0 ]*X1_X0[2 ]);
2116- CP[2 ] = (X2_X1[0 ]*X1_X0[1 ] - X2_X1[1 ]*X1_X0[0 ]);
2067+ Coord = geometry[iMesh]->nodes ->GetCoord (iPoint);
21172068
2118- Distance = sqrt ((CP[ 0 ]*CP[ 0 ]+CP[ 1 ]*CP[ 1 ]+CP[ 2 ]*CP[ 2 ])/(X2_X1[ 0 ]*X2_X1[ 0 ]+X2_X1[ 1 ]*X2_X1[ 1 ]+X2_X1[ 2 ]*X2_X1[ 2 ]) );
2069+ SubsonicEngine_Cyl = config-> GetSubsonicEngine_Cyl ( );
21192070
2120- DotCheck = -(X1_X0[0 ]*X2_X1[0 ]+X1_X0[1 ]*X2_X1[1 ]+X1_X0[2 ]*X2_X1[2 ]);
2121- if (DotCheck < 0.0 ) Distance = sqrt (X1_X0[0 ]*X1_X0[0 ]+X1_X0[1 ]*X1_X0[1 ]+X1_X0[2 ]*X1_X0[2 ]);
2071+ X0[0 ] = Coord[0 ]; X0[1 ] = Coord[1 ]; if (nDim==3 ) X0[2 ] = Coord[2 ];
2072+ X1[0 ] = SubsonicEngine_Cyl[0 ]; X1[1 ] = SubsonicEngine_Cyl[1 ]; X1[2 ] = SubsonicEngine_Cyl[2 ];
2073+ X2[0 ] = SubsonicEngine_Cyl[3 ]; X2[1 ] = SubsonicEngine_Cyl[4 ]; X2[2 ] = SubsonicEngine_Cyl[5 ];
2074+ Radius = SubsonicEngine_Cyl[6 ];
21222075
2123- DotCheck = (X2_X0[0 ]*X2_X1[0 ]+X2_X0[1 ]*X2_X1[1 ]+X2_X0[2 ]*X2_X1[2 ]);
2124- if (DotCheck < 0.0 ) Distance = sqrt (X2_X0[0 ]*X2_X0[0 ]+X2_X0[1 ]*X2_X0[1 ]+X2_X0[2 ]*X2_X0[2 ]);
2076+ GeometryToolbox::Distance (3 , X1, X2, X2_X1);
2077+ GeometryToolbox::Distance (3 , X0, X1, X1_X0);
2078+ GeometryToolbox::Distance (3 , X0, X2, X2_X0);
21252079
2126- if (Distance < Radius) {
2080+ GeometryToolbox::CrossProduct (X2_X1, X1_X0, CP);
21272081
2128- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution (iPoint, 0 , Density_CylND);
2129- for (iDim = 0 ; iDim < nDim; iDim++)
2130- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution (iPoint, iDim+1 , Density_CylND*Velocity_CylND[iDim]);
2131- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution (iPoint, nVar-1 , Density_CylND*Energy_CylND);
2082+ Distance = sqrt (GeometryToolbox::SquaredNorm (3 ,CP) / GeometryToolbox::SquaredNorm (3 ,X2_X1));
21322083
2133- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution_Old (iPoint, 0 , Density_CylND);
2134- for (iDim = 0 ; iDim < nDim; iDim++)
2135- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution_Old (iPoint, iDim+1 , Density_CylND*Velocity_CylND[iDim]);
2136- solver_container[iMesh][FLOW_SOL]->GetNodes ()->SetSolution_Old (iPoint, nVar-1 , Density_CylND*Energy_CylND);
2084+ DotCheck = -GeometryToolbox::DotProduct (3 , X1_X0, X2_X1);
2085+ if (DotCheck < 0.0 ) Distance = GeometryToolbox::Norm (3 , X1_X0);
21372086
2138- }
2087+ DotCheck = GeometryToolbox::DotProduct (3 , X2_X0, X2_X1);
2088+ if (DotCheck < 0.0 ) Distance = GeometryToolbox::Norm (3 , X2_X0);
21392089
2090+ if (Distance < Radius) {
2091+ FlowNodes->SetSolution (iPoint, 0 , Density_CylND);
2092+ for (iDim = 0 ; iDim < nDim; iDim++)
2093+ FlowNodes->SetSolution (iPoint, iDim+1 , Density_CylND*Velocity_CylND[iDim]);
2094+ FlowNodes->SetSolution (iPoint, nVar-1 , Density_CylND*Energy_CylND);
21402095 }
21412096
2142- /* --- Set the MPI communication ---*/
2143-
2144- solver_container[iMesh][FLOW_SOL]->InitiateComms (geometry[iMesh], config, SOLUTION);
2145- solver_container[iMesh][FLOW_SOL]->CompleteComms (geometry[iMesh], config, SOLUTION);
2146-
2147- solver_container[iMesh][FLOW_SOL]->InitiateComms (geometry[iMesh], config, SOLUTION_OLD);
2148- solver_container[iMesh][FLOW_SOL]->CompleteComms (geometry[iMesh], config, SOLUTION_OLD);
2149-
21502097 }
21512098
2152- }
2153-
2154- }
2099+ FlowNodes->Set_OldSolution ();
21552100
2156- /* --- Make sure that the solution is well initialized for unsteady
2157- calculations with dual time-stepping (load additional restarts for 2nd-order). ---*/
2158-
2159- if (dual_time && ((TimeIter == 0 ) || (restart && (TimeIter == config->GetRestart_Iter ()))) ) {
2160- PushSolutionBackInTime (TimeIter, restart, rans, solver_container, geometry, config);
2161- }
2101+ }
21622102
21632103 } // end SU2_OMP_PARALLEL
21642104
0 commit comments