From aea19ac9b9e1ba17949dcdad510e515397f892b8 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Tue, 19 Feb 2019 20:24:42 +0100 Subject: [PATCH 1/8] Bugfix sym plane for compressible and incompressible solver. Working but cleaning necessary. --- SU2_CFD/src/solver_direct_mean.cpp | 226 ++++++++++++++++++++++- SU2_CFD/src/solver_direct_mean_inc.cpp | 236 ++++++++++++++++++++++++- 2 files changed, 454 insertions(+), 8 deletions(-) mode change 100644 => 100755 SU2_CFD/src/solver_direct_mean.cpp mode change 100644 => 100755 SU2_CFD/src/solver_direct_mean_inc.cpp diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp old mode 100644 new mode 100755 index 03f919c07b22..c810175d4bed --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -12481,12 +12481,230 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta } -void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, - CConfig *config, unsigned short val_marker) { +void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) { - /*--- Call the Euler residual ---*/ + unsigned short iDim, iVar; + unsigned long iVertex, iPoint, Point_Normal; + + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + + su2double *Normal = new su2double[nDim]; + + su2double UnitNormal[nDim], NormalArea[nDim], Velocity_i[nDim], Tangential[nDim]; + su2double ProjVelocity_i = 0.0; + + su2double ProjGradient[nDim+2]; + su2double *V_reflected, *V_domain; + su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; + su2double **Grad_Prim = new su2double*[nPrimVarGrad]; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + Grad_Reflected[iVar] = new su2double[nDim]; + Grad_Prim[iVar] = new su2double[nDim]; + } - BC_Euler_Wall(geometry, solver_container, conv_numerics, config, val_marker); + /*--- Loop over all the vertices on this boundary marker ---*/ + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + /*--- Allocate the value at the outlet ---*/ + + V_reflected = GetCharacPrimVar(val_marker, iVertex); //TK no idea why this was done like that in inc solver + + iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ + + if (geometry->node[iPoint]->GetDomain()) { + + /*--- Index of the closest interior node ---*/ + + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + + /*--- Normal vector for this vertex (negate for outward convention) ---*/ + + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); + + if (config->GetGrid_Movement()) + conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), + geometry->node[iPoint]->GetGridVel()); + + /*--- Compute unit Normal, to be used for projected velocity ---*/ + + su2double Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Area += Normal[iDim]*Normal[iDim]; + Area = sqrt (Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; + + /*--- Current solution at this boundary node ---*/ + + V_domain = node[iPoint]->GetPrimitive(); + + /*--- Force the velocity to be tangential to the surface. ---*/ + + ProjVelocity_i = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = node[iPoint]->GetVelocity(iDim); + ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; + } + + for (iDim = 0; iDim < nDim; iDim++) { + V_reflected[iDim+1] = Velocity_i[iDim] - 2.0*ProjVelocity_i * UnitNormal[iDim]; + } + + V_reflected[0] = node[iPoint]->GetTemperature(); + V_reflected[nDim+1] = node[iPoint]->GetPressure(); + V_reflected[nDim+2] = node[iPoint]->GetDensity(); + V_reflected[nDim+3] = node[iPoint]->GetEnthalpy(); + //TK JST also uses nDim+4 speed of sound + V_reflected[nDim+8] = node[iPoint]->GetSpecificHeatCp(); + + + /*--- Set various quantities in the solver class ---*/ + + conv_numerics->SetPrimitive(V_domain, V_reflected); + conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary());//TK What are secondarys: thermo-physical properties (partial derivatives) + + /*--- Compute the residual using an upwind scheme ---*/ + + conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); + + /*--- Update residual value ---*/ + + LinSysRes.AddBlock(iPoint, Residual); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) { //TK removed jacobian entries + Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); + } + + /*--- own approach for 2D: enhanced velocity gradient computation---*/ + if (config->GetViscous()) { + + /*--- Set transport properties at the outlet. ---*/ + + V_reflected[nDim+5] = node[iPoint]->GetLaminarViscosity(); + V_reflected[nDim+6] = node[iPoint]->GetEddyViscosity(); + + /*--- Set the normal vector and the coordinates ---*/ + + visc_numerics->SetNormal(Normal); + visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); //TK + + /*--- Primitive variables, and gradient ---*/ + + visc_numerics->SetPrimitive(V_domain, V_reflected); + visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); + + /*--- Get gradients of primitives of boundary cell ---*/ //TK node[iPoint]->GetGradient_Primitive() gets whole matrix at once.. well it gets the pointer + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + for (iDim = 0; iDim < nDim; iDim++) { + Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); + } + } + + /*--- Set gradients of scalars p,T ---*/ + + /*--- Ensure that the normal gradients are also reflected for ... its not these vars any more [p, v_x, v_y, (v_z), T]. ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { //only v and T are necessary + /*--- Compute projected part of the gradient in a dot product ---*/ + ProjGradient[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjGradient[iVar] += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; + /*--- Compute reflected state of the gradients + TK Gradients of velocity components are set here but overridden later ---*/ + for (iDim = 0; iDim < nDim; iDim++) + Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient[iVar]*UnitNormal[iDim]; //TK ProjGradient could be scalar + } + + /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ + + switch( nDim ) { + case 2: { + Tangential[0] = -UnitNormal[1]; + Tangential[1] = UnitNormal[0]; + break; + } + case 3: { + Tangential[0] = -UnitNormal[1]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); + Tangential[1] = UnitNormal[0]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); + Tangential[2] = 0.0; + break; + } + } + + /*--- Compute gradients of normal and tangential velocity ---*/ + + su2double GradNormVel[nDim]; + su2double GradTangVel[nDim]; + for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components + GradNormVel[iVar] = 0.0; + GradTangVel[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normaal + GradNormVel[iVar] += Grad_Prim[iDim+1][iVar] * UnitNormal[iDim]; + GradTangVel[iVar] += Grad_Prim[iDim+1][iVar] * Tangential[iDim]; + } + } + + /*--- Refelect gradients in tangential and normal direction ---*/ + + su2double ReflGradNormVel[nDim]; + su2double ReflGradTangVel[nDim]; + su2double ProjNormVelGrad = 0.0; + su2double ProjTangVelGrad = 0.0; + + for (iDim = 0; iDim < nDim; iDim++) { + ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; + ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; + } + + for (iDim = 0; iDim < nDim; iDim++) { + ReflGradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; + ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; + } + + /*--- Compute Cartesian reflected gradients ---*/ + + for (iVar = 0; iVar < nDim; iVar++) { // loops over the velocity component gradients + for (iDim = 0; iDim < nDim; iDim++) { // loops over the entries of the above + Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; + } + } + + /*--- End own enhanced ---*/ + visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); + + /*--- Turbulent kinetic energy ---*/ + + if (config->GetKind_Turb_Model() == SST) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), + solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); + + /*--- Compute and update residual. Note that the shear stress tensor is computed in the + * following routine based upon the velocity-component gradients. ---*/ + + visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); + + LinSysRes.SubtractBlock(iPoint, Residual); + + /*--- Jacobian contribution for implicit integration ---*/ + if (implicit) + Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); + + } + } + } } diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp old mode 100644 new mode 100755 index ee6a6184c20f..d54b061879f7 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6050,14 +6050,242 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } -void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, - CConfig *config, unsigned short val_marker) { +void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, + unsigned short val_marker) { - /*--- Call the Euler wall residual method. ---*/ + unsigned short iDim, iVar; + unsigned long iVertex, iPoint, Point_Normal; + + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - BC_Euler_Wall(geometry, solver_container, conv_numerics, config, val_marker); + su2double UnitNormal[nDim], Velocity_i[nDim], Tangential[nDim]; + su2double ProjVelocity_i = 0.0; + su2double ProjGradient[nDim+2]; //TK here nDim+2 but Grad_Reflected is nPrimVarGrad big + su2double *V_reflected, *V_domain; + su2double *Normal = new su2double[nDim]; + su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; + su2double **Grad_Prim = new su2double*[nPrimVarGrad]; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + Grad_Reflected[iVar] = new su2double[nDim]; + Grad_Prim[iVar] = new su2double[nDim]; + } + + /*--- Loop over all the vertices on this boundary marker ---*/ + + for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + /*--- Allocate the value at the outlet ---*/ + + V_reflected = GetCharacPrimVar(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->node[iPoint]->GetDomain()) { + + /*--- Index of the closest interior node ---*/ + + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); //Could be used for SetCoord for viscous jacobian (dist_ij) + + /*--- Normal vector for this vertex (negate for outward convention) ---*/ + + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); + + if (config->GetGrid_Movement()) + conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), + geometry->node[iPoint]->GetGridVel()); + + /*--- Compute unit Normal, to be used for projected velocity ---*/ + + su2double Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Area += Normal[iDim]*Normal[iDim]; + Area = sqrt(Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; + + /*--- Current solution at this boundary node ---*/ + + V_domain = node[iPoint]->GetPrimitive(); + + /*--- Set the reflected state based on the boundary node. Scalars are copied and + the velocity is mirrored along the symmetry boundary, i.e. the velocity in + normal direction is substracted twice. ---*/ + + ProjVelocity_i = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = node[iPoint]->GetVelocity(iDim); + ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; + } + + for (iDim = 0; iDim < nDim; iDim++) { + V_reflected[iDim+1] = Velocity_i[iDim] - 2.0 * ProjVelocity_i*UnitNormal[iDim]; + } + + V_reflected[0] = node[iPoint]->GetPressure(); + V_reflected[nDim+1] = node[iPoint]->GetTemperature(); + V_reflected[nDim+2] = node[iPoint]->GetDensity(); + V_reflected[nDim+3] = node[iPoint]->GetBetaInc2(); + V_reflected[nDim+7] = node[iPoint]->GetSpecificHeatCp(); + + conv_numerics->SetPrimitive(V_domain, V_reflected); + + /*--- Compute the residual using an upwind scheme ---*/ + + conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); + + /*--- Update residual value ---*/ + + LinSysRes.AddBlock(iPoint, Residual); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) { + Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); + } + + /*--- ---*/ + if (config->GetViscous()) { + + /*--- Set the normal vector and the coordinates ---*/ + + visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); + visc_numerics->SetNormal(Normal); + + /*--- Set transport properties at the outlet. ---*/ + + V_reflected[nDim+4] = node[iPoint]->GetLaminarViscosity(); + V_reflected[nDim+5] = node[iPoint]->GetEddyViscosity(); + V_reflected[nDim+6] = node[iPoint]->GetThermalConductivity(); + + /*--- Primitive variables ---*/ + + visc_numerics->SetPrimitive(V_domain, V_reflected); + + /*--- For viscous Fluxes also the gradients of the primitives need to be determined. + 1. The gradients of scalars are mirrored along the sym plane + 2. The gradients of the velocity components need more attention ---*/ + + /*--- Get gradients of primitives of boundary cell ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) + for (iDim = 0; iDim < nDim; iDim++) + Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); + + /*--- Ensure that the gradients are also reflected for [p, v_x, v_y, (v_z), T]. ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + /*--- Compute projected part of the gradient in a dot product ---*/ + ProjGradient[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjGradient[iVar] += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; + /*--- Compute reflected state of the gradients. + TK Gradients of velocity components are set here but overridden later ---*/ + for (iDim = 0; iDim < nDim; iDim++) + Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient[iVar]*UnitNormal[iDim]; //TK ProjGradient could be scalar + } + + /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ + + switch( nDim ) { + case 2: { + Tangential[0] = -UnitNormal[1]; + Tangential[1] = UnitNormal[0]; + break; + } + case 3: { + Tangential[0] = -UnitNormal[1]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); + Tangential[1] = UnitNormal[0]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); + Tangential[2] = 0.0; + break; + } + } + + /*--- Compute gradients of normal and tangential velocity ---*/ + + su2double GradNormVel[nDim]; + su2double GradTangVel[nDim]; + for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components + GradNormVel[iVar] = 0.0; + GradTangVel[iVar] = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normaal + GradNormVel[iVar] += Grad_Prim[iDim+1][iVar] * UnitNormal[iDim]; + GradTangVel[iVar] += Grad_Prim[iDim+1][iVar] * Tangential[iDim]; + } + } + + /*--- Refelect gradients in tangential and normal direction ---*/ + + su2double ReflGradNormVel[nDim]; + su2double ReflGradTangVel[nDim]; + su2double ProjNormVelGrad = 0.0; + su2double ProjTangVelGrad = 0.0; + + for (iDim = 0; iDim < nDim; iDim++) { + ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; + ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; + } + + for (iDim = 0; iDim < nDim; iDim++) { + ReflGradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; + ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; + } + + /*--- Compute Cartesian reflected gradients ---*/ + + for (iVar = 0; iVar < nDim; iVar++) { // loops over the velocity component gradients + for (iDim = 0; iDim < nDim; iDim++) { // loops over the entries of the above + Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; + } + } + + /*--- Gradient and limiters ---*/ + + visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), + Grad_Reflected); + + /*--- Turbulent kinetic energy ---*/ + + if (config->GetKind_Turb_Model() == SST) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), + solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); + + /*--- Compute and update residual. Note that the shear stress tensor is computed in the + following routine based upon the velocity-component gradients. ---*/ + + visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); + + LinSysRes.SubtractBlock(iPoint, Residual); + + /*--- Jacobian contribution for implicit integration ---*/ + + if (implicit) + Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); + + } // Viscous + + } + } + + /*--- Free locally allocated memory ---*/ + + delete [] Normal; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + delete [] Grad_Prim[iVar]; + delete [] Grad_Reflected[iVar]; + } + delete [] Grad_Prim; + delete [] Grad_Reflected; } + void CIncEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config) { From 1f72e5ae1e8e47c15c0ae9dd75fa2563c975f849 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 21 Feb 2019 17:34:11 +0100 Subject: [PATCH 2/8] Added documentation. Cleaning. --- .travis.yml | 4 +- SU2_CFD/src/solver_direct_mean.cpp | 125 +++++++++++++++---------- SU2_CFD/src/solver_direct_mean_inc.cpp | 93 ++++++++++-------- 3 files changed, 134 insertions(+), 88 deletions(-) diff --git a/.travis.yml b/.travis.yml index b974c79b1279..25d208b369af 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,11 +12,11 @@ compiler: notifications: email: recipients: - - su2code-dev@lists.stanford.edu + - tobias.kattmann@de.bosch.com branches: only: - - develop + - bugfix_BC-sym-plane python: - 2.7 diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index d33837e3a55d..cce3b0f28dfb 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -12496,13 +12496,12 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - su2double *Normal = new su2double[nDim]; - su2double UnitNormal[nDim], NormalArea[nDim], Velocity_i[nDim], Tangential[nDim]; su2double ProjVelocity_i = 0.0; - su2double ProjGradient[nDim+2]; - su2double *V_reflected, *V_domain; + su2double ProjGradient; + su2double *Normal = new su2double[nDim]; + su2double *V_reflected, *V_domain;// TK why here no new su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; su2double **Grad_Prim = new su2double*[nPrimVarGrad]; for (iVar = 0; iVar < nPrimVarGrad; iVar++) { @@ -12514,9 +12513,9 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Allocate the value at the outlet ---*/ + /*--- Allocate the reflected state at the symmetry boundary ---*/ - V_reflected = GetCharacPrimVar(val_marker, iVertex); //TK no idea why this was done like that in inc solver + V_reflected = GetCharacPrimVar(val_marker, iVertex); // TK Why this method iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -12534,11 +12533,13 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); + /*--- Grid movement ---*/ + if (config->GetGrid_Movement()) conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); - /*--- Compute unit Normal, to be used for projected velocity ---*/ + /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ su2double Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) @@ -12551,8 +12552,18 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Current solution at this boundary node ---*/ V_domain = node[iPoint]->GetPrimitive(); - - /*--- Force the velocity to be tangential to the surface. ---*/ + + /*--- Set the reflected state based on the boundary node. Scalars are copied and + the velocity is mirrored along the symmetry boundary, i.e. the velocity in + normal direction is substracted twice. ---*/ + + /*--- T, vx, vy, vz, p, rho, enthalpy, speed of sound, lam visco, eddy visco, ?, specific heat cp, ? ---*/ + + for(iVar = 0; iVar < nPrimVar; iVar++) + V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); + + /*--- Compute velocity in normal direction (ProjVelcity_i) und substract twice from + velocity in normal direction. ---*/ ProjVelocity_i = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -12562,20 +12573,12 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) { V_reflected[iDim+1] = Velocity_i[iDim] - 2.0*ProjVelocity_i * UnitNormal[iDim]; - } - - V_reflected[0] = node[iPoint]->GetTemperature(); - V_reflected[nDim+1] = node[iPoint]->GetPressure(); - V_reflected[nDim+2] = node[iPoint]->GetDensity(); - V_reflected[nDim+3] = node[iPoint]->GetEnthalpy(); - //TK JST also uses nDim+4 speed of sound - V_reflected[nDim+8] = node[iPoint]->GetSpecificHeatCp(); - + } /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_reflected); - conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary());//TK What are secondarys: thermo-physical properties (partial derivatives) + conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); //TK What are secondarys: thermo-physical properties (partial derivatives) /*--- Compute the residual using an upwind scheme ---*/ @@ -12587,48 +12590,57 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) { //TK removed jacobian entries + if (implicit) { Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - /*--- own approach for 2D: enhanced velocity gradient computation---*/ - if (config->GetViscous()) { + /*--- The viscous fluxes of the Navier-Stokes equations depend on the Primitive + variables and thier gradients. The viscous numerics container is filled + just as the convective numerics container, but the primitive gradients of + the reflected state have be determined such that symmetry at the boundary is + enforced. ---*/ - /*--- Set transport properties at the outlet. ---*/ - - V_reflected[nDim+5] = node[iPoint]->GetLaminarViscosity(); - V_reflected[nDim+6] = node[iPoint]->GetEddyViscosity(); + if (config->GetViscous()) { /*--- Set the normal vector and the coordinates ---*/ visc_numerics->SetNormal(Normal); - visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); //TK - - /*--- Primitive variables, and gradient ---*/ + visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); + + /*--- Set the primitive and Secondary variables ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); - /*--- Get gradients of primitives of boundary cell ---*/ //TK node[iPoint]->GetGradient_Primitive() gets whole matrix at once.. well it gets the pointer - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - for (iDim = 0; iDim < nDim; iDim++) { + /*--- For viscous Fluxes also the gradients of the primitives need to be determined. + 1. The gradients of scalars are mirrored along the sym plane just as velocity for the primitives + 2. The gradients of the velocity components need more attention, i.e. the gradient of the + normal velocity in tangential direction is mirrored and the gradient of the tangential velocity in + normal direction is mirrored. ---*/ + + /*--- Get gradients of primitives of boundary cell ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) + for (iDim = 0; iDim < nDim; iDim++) Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); - } - } - - /*--- Set gradients of scalars p,T ---*/ - - /*--- Ensure that the normal gradients are also reflected for ... its not these vars any more [p, v_x, v_y, (v_z), T]. ---*/ - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { //only v and T are necessary + + /*--- Reflect the gradients for all scalars including the velocity components. + The gradients of the velocity components are overriden later with the + correct values. ---*/ + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient[iVar] = 0.0; + + ProjGradient = 0.0; for (iDim = 0; iDim < nDim; iDim++) - ProjGradient[iVar] += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; + ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; + /*--- Compute reflected state of the gradients TK Gradients of velocity components are set here but overridden later ---*/ + for (iDim = 0; iDim < nDim; iDim++) - Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient[iVar]*UnitNormal[iDim]; //TK ProjGradient could be scalar - } + Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; + } /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ @@ -12639,9 +12651,16 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, break; } case 3: { - Tangential[0] = -UnitNormal[1]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); - Tangential[1] = UnitNormal[0]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); - Tangential[2] = 0.0; + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that ---*/ + + unsigned short Largest, Arbitrary, Zero; + if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} + else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} + else {Largest=2;Arbitrary=1;Zero=0;} + + Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Zero] = 0.0; break; } } @@ -12689,7 +12708,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Turbulent kinetic energy ---*/ - if (config->GetKind_Turb_Model() == SST) + if (config->GetKind_Turb_Model() == SST) //TK check turb solver in visc residual of direct mean visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); @@ -12703,11 +12722,19 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Jacobian contribution for implicit integration ---*/ if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); - } } } - + + /*--- Free locally allocated memory ---*/ + + delete [] Normal; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + delete [] Grad_Prim[iVar]; + delete [] Grad_Reflected[iVar]; + } + delete [] Grad_Prim; + delete [] Grad_Reflected; } void CEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index cd27129b04de..6fc521709713 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6067,9 +6067,8 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, string Marker_Tag = config->GetMarker_All_TagBound(val_marker); su2double UnitNormal[nDim], Velocity_i[nDim], Tangential[nDim]; - su2double ProjVelocity_i = 0.0; + su2double ProjVelocity_i = 0.0, ProjGradient = 0.0; - su2double ProjGradient[nDim+2]; //TK here nDim+2 but Grad_Reflected is nPrimVarGrad big su2double *V_reflected, *V_domain; su2double *Normal = new su2double[nDim]; su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; @@ -6083,7 +6082,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Allocate the value at the outlet ---*/ + /*--- Allocate the reflected state at the symmetry boundary ---*/ V_reflected = GetCharacPrimVar(val_marker, iVertex); @@ -6095,7 +6094,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Index of the closest interior node ---*/ - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); //Could be used for SetCoord for viscous jacobian (dist_ij) + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); // TK Could be used for SetCoord for viscous jacobian (dist_ij) /*--- Normal vector for this vertex (negate for outward convention) ---*/ @@ -6103,11 +6102,13 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); + /*--- Grid movement ---*/ + if (config->GetGrid_Movement()) conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); - /*--- Compute unit Normal, to be used for projected velocity ---*/ + /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ su2double Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) @@ -6124,7 +6125,15 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Set the reflected state based on the boundary node. Scalars are copied and the velocity is mirrored along the symmetry boundary, i.e. the velocity in normal direction is substracted twice. ---*/ - + + /*--- p, vx, vy, vz, T, rho, betainc2, lam visco, eddy visco, thermal cond, specific heat cp, ?, ? ---*/ + + for(iVar = 0; iVar < nPrimVar; iVar++) + V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); + + /*--- Compute velocity in normal direction (ProjVelcity_i) und substract twice from + velocity in normal direction. ---*/ + ProjVelocity_i = 0.0; for (iDim = 0; iDim < nDim; iDim++) { Velocity_i[iDim] = node[iPoint]->GetVelocity(iDim); @@ -6135,14 +6144,11 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, V_reflected[iDim+1] = Velocity_i[iDim] - 2.0 * ProjVelocity_i*UnitNormal[iDim]; } - V_reflected[0] = node[iPoint]->GetPressure(); - V_reflected[nDim+1] = node[iPoint]->GetTemperature(); - V_reflected[nDim+2] = node[iPoint]->GetDensity(); - V_reflected[nDim+3] = node[iPoint]->GetBetaInc2(); - V_reflected[nDim+7] = node[iPoint]->GetSpecificHeatCp(); + /*--- Set Primitive and Secondary for numerics class ---*/ conv_numerics->SetPrimitive(V_domain, V_reflected); - + conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); + /*--- Compute the residual using an upwind scheme ---*/ conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); @@ -6157,43 +6163,52 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - /*--- ---*/ + /*--- The viscous fluxes of the Navier-Stokes equations depend on the Primitive + variables and thier gradients. The viscous numerics container is filled + just as the convective numerics container, but the primitive gradients of + the reflected state have be determined such that symmetry at the boundary is + enforced. ---*/ + if (config->GetViscous()) { /*--- Set the normal vector and the coordinates ---*/ visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); visc_numerics->SetNormal(Normal); - - /*--- Set transport properties at the outlet. ---*/ - - V_reflected[nDim+4] = node[iPoint]->GetLaminarViscosity(); - V_reflected[nDim+5] = node[iPoint]->GetEddyViscosity(); - V_reflected[nDim+6] = node[iPoint]->GetThermalConductivity(); - /*--- Primitive variables ---*/ + /*--- Set the primitive and Secondary variables ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); + visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); /*--- For viscous Fluxes also the gradients of the primitives need to be determined. - 1. The gradients of scalars are mirrored along the sym plane - 2. The gradients of the velocity components need more attention ---*/ + 1. The gradients of scalars are mirrored along the sym plane just as velocity for the primitives + 2. The gradients of the velocity components need more attention, i.e. the gradient of the + normal velocity in tangential direction is mirrored and the gradient of the tangential velocity in + normal direction is mirrored. ---*/ /*--- Get gradients of primitives of boundary cell ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); - /*--- Ensure that the gradients are also reflected for [p, v_x, v_y, (v_z), T]. ---*/ + /*--- Reflect the gradients for all scalars including the velocity components. + The gradients of the velocity components are overriden later with the + correct values. ---*/ + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient[iVar] = 0.0; + + ProjGradient = 0.0; for (iDim = 0; iDim < nDim; iDim++) - ProjGradient[iVar] += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; - /*--- Compute reflected state of the gradients. - TK Gradients of velocity components are set here but overridden later ---*/ + ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; + + /*--- Compute reflected state of the gradients. Gradients of velocity components are set here but overridden later ---*/ + for (iDim = 0; iDim < nDim; iDim++) - Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient[iVar]*UnitNormal[iDim]; //TK ProjGradient could be scalar + Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; } /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ @@ -6205,9 +6220,16 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, break; } case 3: { - Tangential[0] = -UnitNormal[1]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); - Tangential[1] = UnitNormal[0]/sqrt(pow(UnitNormal[0],2) + pow(UnitNormal[1],2)); - Tangential[2] = 0.0; + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that ---*/ + + unsigned short Largest, Arbitrary, Zero; + if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} + else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} + else {Largest=2;Arbitrary=1;Zero=0;} + + Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Zero] = 0.0; break; } } @@ -6252,8 +6274,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Gradient and limiters ---*/ - visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), - Grad_Reflected); + visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); /*--- Turbulent kinetic energy ---*/ @@ -6272,15 +6293,13 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); - - } // Viscous - + } } } /*--- Free locally allocated memory ---*/ - delete [] Normal; + delete [] Normal; for (iVar = 0; iVar < nPrimVarGrad; iVar++) { delete [] Grad_Prim[iVar]; delete [] Grad_Reflected[iVar]; From 6662b7c12138565d1960f3f62350064a0c18fab8 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 21 Feb 2019 20:11:39 +0100 Subject: [PATCH 3/8] Little style changes. --- SU2_CFD/src/solver_direct_mean.cpp | 21 +++++++++------------ SU2_CFD/src/solver_direct_mean_inc.cpp | 21 +++++++++------------ 2 files changed, 18 insertions(+), 24 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index cce3b0f28dfb..eb2587fb8d7a 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -12483,11 +12483,11 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta } -void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, - CSolver **solver_container, - CNumerics *conv_numerics, - CNumerics *visc_numerics, - CConfig *config, +void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, unsigned short val_marker) { unsigned short iDim, iVar; @@ -12536,8 +12536,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Grid movement ---*/ if (config->GetGrid_Movement()) - conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), - geometry->node[iPoint]->GetGridVel()); + conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ @@ -12647,7 +12646,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, switch( nDim ) { case 2: { Tangential[0] = -UnitNormal[1]; - Tangential[1] = UnitNormal[0]; + Tangential[1] = UnitNormal[0]; break; } case 3: { @@ -12697,11 +12696,9 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Compute Cartesian reflected gradients ---*/ - for (iVar = 0; iVar < nDim; iVar++) { // loops over the velocity component gradients - for (iDim = 0; iDim < nDim; iDim++) { // loops over the entries of the above + for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients + for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; - } - } /*--- End own enhanced ---*/ visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 6fc521709713..95cbcfb75ade 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6053,11 +6053,11 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, } -void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, - CSolver **solver_container, - CNumerics *conv_numerics, - CNumerics *visc_numerics, - CConfig *config, +void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, + CSolver **solver_container, + CNumerics *conv_numerics, + CNumerics *visc_numerics, + CConfig *config, unsigned short val_marker) { unsigned short iDim, iVar; @@ -6105,8 +6105,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Grid movement ---*/ if (config->GetGrid_Movement()) - conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), - geometry->node[iPoint]->GetGridVel()); + conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ @@ -6216,7 +6215,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, switch( nDim ) { case 2: { Tangential[0] = -UnitNormal[1]; - Tangential[1] = UnitNormal[0]; + Tangential[1] = UnitNormal[0]; break; } case 3: { @@ -6266,11 +6265,9 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Compute Cartesian reflected gradients ---*/ - for (iVar = 0; iVar < nDim; iVar++) { // loops over the velocity component gradients - for (iDim = 0; iDim < nDim; iDim++) { // loops over the entries of the above + for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients + for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; - } - } /*--- Gradient and limiters ---*/ From 9d9e556c6026f233fd5ef4e0eec46c645f0246c1 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Fri, 22 Feb 2019 15:03:21 +0100 Subject: [PATCH 4/8] More comments and additional code cleaning. --- SU2_CFD/src/solver_direct_mean.cpp | 164 +++++++++++------------- SU2_CFD/src/solver_direct_mean_inc.cpp | 168 +++++++++++-------------- 2 files changed, 151 insertions(+), 181 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index eb2587fb8d7a..9586370d3e00 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -12491,55 +12491,55 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, unsigned short val_marker) { unsigned short iDim, iVar; - unsigned long iVertex, iPoint, Point_Normal; - - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + unsigned long iVertex, iPoint; - su2double UnitNormal[nDim], NormalArea[nDim], Velocity_i[nDim], Tangential[nDim]; - su2double ProjVelocity_i = 0.0; + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - su2double ProjGradient; - su2double *Normal = new su2double[nDim]; - su2double *V_reflected, *V_domain;// TK why here no new + su2double ProjVelocity_i, ProjGradient; + su2double *V_reflected, *V_domain; + + su2double *Normal = new su2double[nDim]; + su2double *UnitNormal = new su2double[nDim]; + su2double *Tangential = new su2double[nDim]; + + /*--- Allocation of primitive gradient arrays. ---*/ su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; - su2double **Grad_Prim = new su2double*[nPrimVarGrad]; + su2double **Grad_Prim = new su2double*[nPrimVarGrad]; for (iVar = 0; iVar < nPrimVarGrad; iVar++) { Grad_Reflected[iVar] = new su2double[nDim]; - Grad_Prim[iVar] = new su2double[nDim]; + Grad_Prim[iVar] = new su2double[nDim]; } - /*--- Loop over all the vertices on this boundary marker ---*/ - + /*--- Loop over all the vertices on this boundary marker. ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Allocate the reflected state at the symmetry boundary ---*/ - - V_reflected = GetCharacPrimVar(val_marker, iVertex); // TK Why this method - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ if (geometry->node[iPoint]->GetDomain()) { - /*--- Index of the closest interior node ---*/ - - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ + /*-------------------------------------------------------------------------------*/ + /*--- Step 1: For the convective fluxes, create a reflected state of the ---*/ + /*--- Primitive variables by copying all interior values to the ---*/ + /*--- reflected. Only the velocity is mirrored along the symmetry ---*/ + /*--- axis. Based on the Upwind_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ + + /*--- Allocate the reflected state at the symmetry boundary. ---*/ + V_reflected = GetCharacPrimVar(val_marker, iVertex); + + /*--- Grid movement ---*/ + if (config->GetGrid_Movement()) + conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); + /*--- Normal vector for this vertex (negate for outward convention). ---*/ geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + for (iDim = 0; iDim < nDim; iDim++) + Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); - /*--- Grid movement ---*/ - - if (config->GetGrid_Movement()) - conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); - - /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ - + /*--- Compute unit normal, to be used for projected velocity and velocity component gradients. ---*/ su2double Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; @@ -12548,66 +12548,54 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Current solution at this boundary node ---*/ - + /*--- Get current solution at this boundary node ---*/ V_domain = node[iPoint]->GetPrimitive(); /*--- Set the reflected state based on the boundary node. Scalars are copied and the velocity is mirrored along the symmetry boundary, i.e. the velocity in normal direction is substracted twice. ---*/ - - /*--- T, vx, vy, vz, p, rho, enthalpy, speed of sound, lam visco, eddy visco, ?, specific heat cp, ? ---*/ - + ProjVelocity_i = 0.0; for(iVar = 0; iVar < nPrimVar; iVar++) V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); - /*--- Compute velocity in normal direction (ProjVelcity_i) und substract twice from - velocity in normal direction. ---*/ - - ProjVelocity_i = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = node[iPoint]->GetVelocity(iDim); - ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; - } - - for (iDim = 0; iDim < nDim; iDim++) { - V_reflected[iDim+1] = Velocity_i[iDim] - 2.0*ProjVelocity_i * UnitNormal[iDim]; - } + /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from + velocity, in normal direction: v_r = v - 2 (v*n)n ---*/ + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity_i += node[iPoint]->GetVelocity(iDim)*UnitNormal[iDim]; - /*--- Set various quantities in the solver class ---*/ + for (iDim = 0; iDim < nDim; iDim++) + V_reflected[iDim+1] = node[iPoint]->GetVelocity(iDim) - 2.0 * ProjVelocity_i*UnitNormal[iDim]; + /*--- Set Primitive and Secondary for numerics class. ---*/ conv_numerics->SetPrimitive(V_domain, V_reflected); - conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); //TK What are secondarys: thermo-physical properties (partial derivatives) - - /*--- Compute the residual using an upwind scheme ---*/ + conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); + /*--- Compute the residual using an upwind scheme. ---*/ conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); /*--- Update residual value ---*/ - LinSysRes.AddBlock(iPoint, Residual); - /*--- Jacobian contribution for implicit integration ---*/ - + /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) { Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - /*--- The viscous fluxes of the Navier-Stokes equations depend on the Primitive - variables and thier gradients. The viscous numerics container is filled - just as the convective numerics container, but the primitive gradients of - the reflected state have be determined such that symmetry at the boundary is - enforced. ---*/ - + /*-------------------------------------------------------------------------------*/ + /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ + /*--- Primitive variables and their gradients. The viscous numerics ---*/ + /*--- container is filled just as the convective numerics container, ---*/ + /*--- but the primitive gradients of the reflected state have to be ---*/ + /*--- determined additionally such that symmetry at the boundary is ---*/ + /*--- enforced. Based on the Viscous_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ if (config->GetViscous()) { - /*--- Set the normal vector and the coordinates ---*/ - + /*--- Set the normal vector and the coordinates. ---*/ visc_numerics->SetNormal(Normal); visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); - /*--- Set the primitive and Secondary variables ---*/ - + /*--- Set the primitive and Secondary variables. ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); @@ -12624,25 +12612,19 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Reflect the gradients for all scalars including the velocity components. The gradients of the velocity components are overriden later with the - correct values. ---*/ - + correct values: grad(V)_r = grad(V) - 2 [grad(V)*n]n, V beeing any primitive ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient = 0.0; for (iDim = 0; iDim < nDim; iDim++) ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; - /*--- Compute reflected state of the gradients - TK Gradients of velocity components are set here but overridden later ---*/ - for (iDim = 0; iDim < nDim; iDim++) Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; } /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ - switch( nDim ) { case 2: { Tangential[0] = -UnitNormal[1]; @@ -12650,8 +12632,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, break; } case 3: { - /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that ---*/ - + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ unsigned short Largest, Arbitrary, Zero; if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} @@ -12664,8 +12645,9 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, } } - /*--- Compute gradients of normal and tangential velocity ---*/ - + /*--- Compute gradients of normal and tangential velocity: + grad(v*n) = grad(v_x) n_x + grad(v_y) n_y (+ grad(v_z) n_z) + grad(v*t) = grad(v_x) t_x + grad(v_y) t_y (+ grad(v_z) t_z) ---*/ su2double GradNormVel[nDim]; su2double GradTangVel[nDim]; for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components @@ -12677,16 +12659,18 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, } } - /*--- Refelect gradients in tangential and normal direction ---*/ - + /*--- Refelect gradients in tangential and normal direction by substracting the normal/tangential + component twice, just as done with velocity above. + grad(v*n)_r = grad(v*n) - 2 {grad([v*n])*t}t + grad(v*t)_r = grad(v*t) - 2 {grad([v*t])*n}n ---*/ su2double ReflGradNormVel[nDim]; su2double ReflGradTangVel[nDim]; su2double ProjNormVelGrad = 0.0; su2double ProjTangVelGrad = 0.0; for (iDim = 0; iDim < nDim; iDim++) { - ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; - ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; + ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; //grad([v*n])*t + ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; //grad([v*t])*n } for (iDim = 0; iDim < nDim; iDim++) { @@ -12694,29 +12678,29 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; } - /*--- Compute Cartesian reflected gradients ---*/ - + /*--- Transfer reflected gradients back into the Cartesian Coordinate system: + grad(v_x)_r = grad(v*n)_r n_x + grad(v*t)_r t_x + grad(v_y)_r = grad(v*n)_r n_y + grad(v*t)_r t_y + ( grad(v_z)_r = grad(v*n)_r n_z + grad(v*t)_r t_z ) ---*/ for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; - /*--- End own enhanced ---*/ + /*--- Set the primitive gradients of the boundary and reflected state. ---*/ visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); - /*--- Turbulent kinetic energy ---*/ - - if (config->GetKind_Turb_Model() == SST) //TK check turb solver in visc residual of direct mean + /*--- Turbulent kinetic energy. ---*/ + if (config->GetKind_Turb_Model() == SST) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); /*--- Compute and update residual. Note that the shear stress tensor is computed in the - * following routine based upon the velocity-component gradients. ---*/ - + following routine based upon the velocity-component gradients. ---*/ visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); LinSysRes.SubtractBlock(iPoint, Residual); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); } @@ -12724,8 +12708,10 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, } /*--- Free locally allocated memory ---*/ - delete [] Normal; + delete [] UnitNormal; + delete [] Tangential; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { delete [] Grad_Prim[iVar]; delete [] Grad_Reflected[iVar]; diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 95cbcfb75ade..e948f44e1c13 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6061,54 +6061,54 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, unsigned short val_marker) { unsigned short iDim, iVar; - unsigned long iVertex, iPoint, Point_Normal; - - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + unsigned long iVertex, iPoint; - su2double UnitNormal[nDim], Velocity_i[nDim], Tangential[nDim]; - su2double ProjVelocity_i = 0.0, ProjGradient = 0.0; + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + su2double ProjVelocity_i, ProjGradient; su2double *V_reflected, *V_domain; - su2double *Normal = new su2double[nDim]; + + su2double *Normal = new su2double[nDim]; + su2double *UnitNormal = new su2double[nDim]; + su2double *Tangential = new su2double[nDim]; + + /*--- Allocation of primitive gradient arrays. ---*/ su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; - su2double **Grad_Prim = new su2double*[nPrimVarGrad]; + su2double **Grad_Prim = new su2double*[nPrimVarGrad]; for (iVar = 0; iVar < nPrimVarGrad; iVar++) { Grad_Reflected[iVar] = new su2double[nDim]; - Grad_Prim[iVar] = new su2double[nDim]; + Grad_Prim[iVar] = new su2double[nDim]; } - /*--- Loop over all the vertices on this boundary marker ---*/ - + /*--- Loop over all the vertices on this boundary marker. ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - /*--- Allocate the reflected state at the symmetry boundary ---*/ - - V_reflected = GetCharacPrimVar(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->node[iPoint]->GetDomain()) { - - /*--- Index of the closest interior node ---*/ - - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); // TK Could be used for SetCoord for viscous jacobian (dist_ij) - - /*--- Normal vector for this vertex (negate for outward convention) ---*/ - - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); + + /*-------------------------------------------------------------------------------*/ + /*--- Step 1: For the convective fluxes, create a reflected state of the ---*/ + /*--- Primitive variables by copying all interior values to the ---*/ + /*--- reflected. Only the velocity is mirrored along the symmetry ---*/ + /*--- axis. Based on the Upwind_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ + + /*--- Allocate the reflected state at the symmetry boundary. ---*/ + V_reflected = GetCharacPrimVar(val_marker, iVertex); /*--- Grid movement ---*/ - if (config->GetGrid_Movement()) conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); - /*--- Compute unit normal, to be used for projected velocity and velocity component gradients ---*/ + /*--- Normal vector for this vertex (negate for outward convention). ---*/ + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) + Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); + /*--- Compute unit normal, to be used for projected velocity and velocity component gradients. ---*/ su2double Area = 0.0; for (iDim = 0; iDim < nDim; iDim++) Area += Normal[iDim]*Normal[iDim]; @@ -6117,66 +6117,54 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Current solution at this boundary node ---*/ - + /*--- Get current solution at this boundary node ---*/ V_domain = node[iPoint]->GetPrimitive(); /*--- Set the reflected state based on the boundary node. Scalars are copied and the velocity is mirrored along the symmetry boundary, i.e. the velocity in normal direction is substracted twice. ---*/ - - /*--- p, vx, vy, vz, T, rho, betainc2, lam visco, eddy visco, thermal cond, specific heat cp, ?, ? ---*/ - for(iVar = 0; iVar < nPrimVar; iVar++) V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); - /*--- Compute velocity in normal direction (ProjVelcity_i) und substract twice from - velocity in normal direction. ---*/ - + /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from + velocity in normal direction: v_r = v - 2 (v*n)n ---*/ ProjVelocity_i = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = node[iPoint]->GetVelocity(iDim); - ProjVelocity_i += Velocity_i[iDim]*UnitNormal[iDim]; - } + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity_i += node[iPoint]->GetVelocity(iDim)*UnitNormal[iDim]; - for (iDim = 0; iDim < nDim; iDim++) { - V_reflected[iDim+1] = Velocity_i[iDim] - 2.0 * ProjVelocity_i*UnitNormal[iDim]; - } - - /*--- Set Primitive and Secondary for numerics class ---*/ - + for (iDim = 0; iDim < nDim; iDim++) + V_reflected[iDim+1] = node[iPoint]->GetVelocity(iDim) - 2.0 * ProjVelocity_i*UnitNormal[iDim]; + + /*--- Set Primitive and Secondary for numerics class. ---*/ conv_numerics->SetPrimitive(V_domain, V_reflected); conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); - /*--- Compute the residual using an upwind scheme ---*/ - + /*--- Compute the residual using an upwind scheme. ---*/ conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); - /*--- Update residual value ---*/ - + /*--- Update residual value ---*/ LinSysRes.AddBlock(iPoint, Residual); - /*--- Jacobian contribution for implicit integration ---*/ - + /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) { Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - - /*--- The viscous fluxes of the Navier-Stokes equations depend on the Primitive - variables and thier gradients. The viscous numerics container is filled - just as the convective numerics container, but the primitive gradients of - the reflected state have be determined such that symmetry at the boundary is - enforced. ---*/ - + + /*-------------------------------------------------------------------------------*/ + /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ + /*--- Primitive variables and their gradients. The viscous numerics ---*/ + /*--- container is filled just as the convective numerics container, ---*/ + /*--- but the primitive gradients of the reflected state have to be ---*/ + /*--- determined additionally such that symmetry at the boundary is ---*/ + /*--- enforced. Based on the Viscous_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ if (config->GetViscous()) { - /*--- Set the normal vector and the coordinates ---*/ - + /*--- Set the normal vector and the coordinates. ---*/ visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); visc_numerics->SetNormal(Normal); - /*--- Set the primitive and Secondary variables ---*/ - + /*--- Set the primitive and Secondary variables. ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); @@ -6186,32 +6174,26 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, normal velocity in tangential direction is mirrored and the gradient of the tangential velocity in normal direction is mirrored. ---*/ - /*--- Get gradients of primitives of boundary cell ---*/ - + /*--- Get gradients of primitives of boundary cell ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); /*--- Reflect the gradients for all scalars including the velocity components. The gradients of the velocity components are overriden later with the - correct values. ---*/ - + correct values: grad(V)_r = grad(V) - 2 [grad(V)*n]n, V beeing any primitive ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient = 0.0; for (iDim = 0; iDim < nDim; iDim++) ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; - /*--- Compute reflected state of the gradients. Gradients of velocity components are set here but overridden later ---*/ - for (iDim = 0; iDim < nDim; iDim++) Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; } - /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ - + /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0. ---*/ switch( nDim ) { case 2: { Tangential[0] = -UnitNormal[1]; @@ -6219,8 +6201,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, break; } case 3: { - /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that ---*/ - + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ unsigned short Largest, Arbitrary, Zero; if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} @@ -6233,29 +6214,32 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, } } - /*--- Compute gradients of normal and tangential velocity ---*/ - + /*--- Compute gradients of normal and tangential velocity: + grad(v*n) = grad(v_x) n_x + grad(v_y) n_y (+ grad(v_z) n_z) + grad(v*t) = grad(v_x) t_x + grad(v_y) t_y (+ grad(v_z) t_z) ---*/ su2double GradNormVel[nDim]; su2double GradTangVel[nDim]; for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components GradNormVel[iVar] = 0.0; GradTangVel[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normaal + for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normal/tangential GradNormVel[iVar] += Grad_Prim[iDim+1][iVar] * UnitNormal[iDim]; GradTangVel[iVar] += Grad_Prim[iDim+1][iVar] * Tangential[iDim]; } } - /*--- Refelect gradients in tangential and normal direction ---*/ - + /*--- Refelect gradients in tangential and normal direction by substracting the normal/tangential + component twice, just as done with velocity above. + grad(v*n)_r = grad(v*n) - 2 {grad([v*n])*t}t + grad(v*t)_r = grad(v*t) - 2 {grad([v*t])*n}n ---*/ su2double ReflGradNormVel[nDim]; su2double ReflGradTangVel[nDim]; su2double ProjNormVelGrad = 0.0; su2double ProjTangVelGrad = 0.0; for (iDim = 0; iDim < nDim; iDim++) { - ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; - ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; + ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; //grad([v*n])*t + ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; //grad([v*t])*n } for (iDim = 0; iDim < nDim; iDim++) { @@ -6263,31 +6247,29 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; } - /*--- Compute Cartesian reflected gradients ---*/ - + /*--- Transfer reflected gradients back into the Cartesian Coordinate system: + grad(v_x)_r = grad(v*n)_r n_x + grad(v*t)_r t_x + grad(v_y)_r = grad(v*n)_r n_y + grad(v*t)_r t_y + ( grad(v_z)_r = grad(v*n)_r n_z + grad(v*t)_r t_z ) ---*/ for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; - /*--- Gradient and limiters ---*/ - + /*--- Set the primitive gradients of the boundary and reflected state. ---*/ visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); - /*--- Turbulent kinetic energy ---*/ - + /*--- Turbulent kinetic energy. ---*/ if (config->GetKind_Turb_Model() == SST) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); - /*--- Compute and update residual. Note that the shear stress tensor is computed in the + /*--- Compute and update residual. Note that the viscous shear stress tensor is computed in the following routine based upon the velocity-component gradients. ---*/ - visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); LinSysRes.SubtractBlock(iPoint, Residual); - /*--- Jacobian contribution for implicit integration ---*/ - + /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); } @@ -6295,8 +6277,10 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, } /*--- Free locally allocated memory ---*/ - delete [] Normal; + delete [] UnitNormal; + delete [] Tangential; + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { delete [] Grad_Prim[iVar]; delete [] Grad_Reflected[iVar]; From 784df8c6acdca0f34e3539f34d62577321971bdb Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 28 Feb 2019 12:24:18 +0100 Subject: [PATCH 5/8] Performance improvements. Moved allocations and unnecessary computations out of the vertex loop. --- SU2_CFD/src/solver_direct_mean.cpp | 175 ++++++++++++++----------- SU2_CFD/src/solver_direct_mean_inc.cpp | 175 ++++++++++++++----------- 2 files changed, 190 insertions(+), 160 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index 9586370d3e00..f4150d2da867 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -12495,19 +12495,68 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - su2double ProjVelocity_i, ProjGradient; + /*--- Allocation of variables necessary for convective fluxes. ---*/ + su2double Area, ProjVelocity_i; su2double *V_reflected, *V_domain; - su2double *Normal = new su2double[nDim]; su2double *UnitNormal = new su2double[nDim]; - su2double *Tangential = new su2double[nDim]; - /*--- Allocation of primitive gradient arrays. ---*/ + /*--- Allocation of variables necessary for viscous fluxes. ---*/ + su2double ProjGradient, ProjNormVelGrad, ProjTangVelGrad; + su2double *Tangential = new su2double[nDim]; + su2double *GradNormVel = new su2double[nDim]; + su2double *GradTangVel = new su2double[nDim]; + + /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/ su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; - su2double **Grad_Prim = new su2double*[nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + for (iVar = 0; iVar < nPrimVarGrad; iVar++) Grad_Reflected[iVar] = new su2double[nDim]; - Grad_Prim[iVar] = new su2double[nDim]; + + /*---------------------------------------------------------------------------------------------*/ + /*--- Preprocessing: On a symmetry-plane, the Unit-Normal is constant. Therefore a constant ---*/ + /*--- Unit-Tangential to that Unit-Normal can be prescribed. The computation ---*/ + /*--- of these vectors is done outside the loop (over all Marker-vertices). ---*/ + /*--- The "Normal" in SU2 isan Area-Normal and is most likely not constant ---*/ + /*--- on the symmetry-plane. ---*/ + /*---------------------------------------------------------------------------------------------*/ + + /*--- Normal vector for a random vertex (zero) on this marker (negate for outward convention). ---*/ + geometry->vertex[val_marker][0]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) + Normal[iDim] = -Normal[iDim]; + + /*--- Compute unit normal, to be used for unit tangential, projected velocity and velocity component gradients. ---*/ + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Area += Normal[iDim]*Normal[iDim]; + Area = sqrt (Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; + + /*--- Preprocessing: Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ + if (config->GetViscous()) { + switch( nDim ) { + case 2: { + Tangential[0] = -UnitNormal[1]; + Tangential[1] = UnitNormal[0]; + break; + } + case 3: { + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ + unsigned short Largest, Arbitrary, Zero; + if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && + abs(UnitNormal[0]) >= abs(UnitNormal[2])) {Largest=0;Arbitrary=1;Zero=2;} + else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && + abs(UnitNormal[1]) >= abs(UnitNormal[2])) {Largest=1;Arbitrary=0;Zero=2;} + else {Largest=2;Arbitrary=1;Zero=0;} + + Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Zero] = 0.0; + break; + } + } } /*--- Loop over all the vertices on this boundary marker. ---*/ @@ -12516,7 +12565,6 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->node[iPoint]->GetDomain()) { /*-------------------------------------------------------------------------------*/ @@ -12539,27 +12587,18 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); - /*--- Compute unit normal, to be used for projected velocity and velocity component gradients. ---*/ - su2double Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt (Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Get current solution at this boundary node ---*/ V_domain = node[iPoint]->GetPrimitive(); /*--- Set the reflected state based on the boundary node. Scalars are copied and the velocity is mirrored along the symmetry boundary, i.e. the velocity in normal direction is substracted twice. ---*/ - ProjVelocity_i = 0.0; for(iVar = 0; iVar < nPrimVar; iVar++) V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from - velocity, in normal direction: v_r = v - 2 (v*n)n ---*/ + velocity in normal direction: v_r = v - 2 (v*n)n ---*/ + ProjVelocity_i = 0.0; for (iDim = 0; iDim < nDim; iDim++) ProjVelocity_i += node[iPoint]->GetVelocity(iDim)*UnitNormal[iDim]; @@ -12581,20 +12620,21 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - /*-------------------------------------------------------------------------------*/ - /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ - /*--- Primitive variables and their gradients. The viscous numerics ---*/ - /*--- container is filled just as the convective numerics container, ---*/ - /*--- but the primitive gradients of the reflected state have to be ---*/ - /*--- determined additionally such that symmetry at the boundary is ---*/ - /*--- enforced. Based on the Viscous_Residual routine. ---*/ - /*-------------------------------------------------------------------------------*/ if (config->GetViscous()) { + + /*-------------------------------------------------------------------------------*/ + /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ + /*--- Primitive variables and their gradients. The viscous numerics ---*/ + /*--- container is filled just as the convective numerics container, ---*/ + /*--- but the primitive gradients of the reflected state have to be ---*/ + /*--- determined additionally such that symmetry at the boundary is ---*/ + /*--- enforced. Based on the Viscous_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ /*--- Set the normal vector and the coordinates. ---*/ - visc_numerics->SetNormal(Normal); visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); - + visc_numerics->SetNormal(Normal); + /*--- Set the primitive and Secondary variables. ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); @@ -12608,54 +12648,33 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Get gradients of primitives of boundary cell ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) - Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); - + Grad_Reflected[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); + /*--- Reflect the gradients for all scalars including the velocity components. - The gradients of the velocity components are overriden later with the + The gradients of the velocity components are set later with the correct values: grad(V)_r = grad(V) - 2 [grad(V)*n]n, V beeing any primitive ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + if(iVar == 0 || iVar > nDim) { // Exclude velocity component gradients - /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; - - for (iDim = 0; iDim < nDim; iDim++) - Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; - } - - /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ - switch( nDim ) { - case 2: { - Tangential[0] = -UnitNormal[1]; - Tangential[1] = UnitNormal[0]; - break; - } - case 3: { - /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ - unsigned short Largest, Arbitrary, Zero; - if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} - else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} - else {Largest=2;Arbitrary=1;Zero=0;} - - Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); - Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); - Tangential[Zero] = 0.0; - break; + /*--- Compute projected part of the gradient in a dot product ---*/ + ProjGradient = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjGradient += Grad_Reflected[iVar][iDim]*UnitNormal[iDim]; + + for (iDim = 0; iDim < nDim; iDim++) + Grad_Reflected[iVar][iDim] = Grad_Reflected[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; } } /*--- Compute gradients of normal and tangential velocity: grad(v*n) = grad(v_x) n_x + grad(v_y) n_y (+ grad(v_z) n_z) grad(v*t) = grad(v_x) t_x + grad(v_y) t_y (+ grad(v_z) t_z) ---*/ - su2double GradNormVel[nDim]; - su2double GradTangVel[nDim]; for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components GradNormVel[iVar] = 0.0; GradTangVel[iVar] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normaal - GradNormVel[iVar] += Grad_Prim[iDim+1][iVar] * UnitNormal[iDim]; - GradTangVel[iVar] += Grad_Prim[iDim+1][iVar] * Tangential[iDim]; + for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normal/tangential + GradNormVel[iVar] += Grad_Reflected[iDim+1][iVar] * UnitNormal[iDim]; + GradTangVel[iVar] += Grad_Reflected[iDim+1][iVar] * Tangential[iDim]; } } @@ -12663,19 +12682,16 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, component twice, just as done with velocity above. grad(v*n)_r = grad(v*n) - 2 {grad([v*n])*t}t grad(v*t)_r = grad(v*t) - 2 {grad([v*t])*n}n ---*/ - su2double ReflGradNormVel[nDim]; - su2double ReflGradTangVel[nDim]; - su2double ProjNormVelGrad = 0.0; - su2double ProjTangVelGrad = 0.0; - + ProjNormVelGrad = 0.0; + ProjTangVelGrad = 0.0; for (iDim = 0; iDim < nDim; iDim++) { ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; //grad([v*n])*t ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; //grad([v*t])*n } for (iDim = 0; iDim < nDim; iDim++) { - ReflGradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; - ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; + GradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; + GradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; } /*--- Transfer reflected gradients back into the Cartesian Coordinate system: @@ -12684,7 +12700,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, ( grad(v_z)_r = grad(v*n)_r n_z + grad(v*t)_r t_z ) ---*/ for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above - Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; + Grad_Reflected[iVar+1][iDim] = GradNormVel[iDim]*UnitNormal[iVar] + GradTangVel[iDim]*Tangential[iVar]; /*--- Set the primitive gradients of the boundary and reflected state. ---*/ visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); @@ -12694,7 +12710,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0)); - /*--- Compute and update residual. Note that the shear stress tensor is computed in the + /*--- Compute and update residual. Note that the viscous shear stress tensor is computed in the following routine based upon the velocity-component gradients. ---*/ visc_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); @@ -12704,19 +12720,18 @@ void CEulerSolver::BC_Sym_Plane(CGeometry *geometry, if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); } - } + } } - + /*--- Free locally allocated memory ---*/ delete [] Normal; delete [] UnitNormal; delete [] Tangential; - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - delete [] Grad_Prim[iVar]; + delete [] GradNormVel; + delete [] GradTangVel; + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) delete [] Grad_Reflected[iVar]; - } - delete [] Grad_Prim; delete [] Grad_Reflected; } diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index e948f44e1c13..3a77dc4e0cac 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6065,22 +6065,71 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - su2double ProjVelocity_i, ProjGradient; + /*--- Allocation of variables necessary for convective fluxes. ---*/ + su2double Area, ProjVelocity_i; su2double *V_reflected, *V_domain; - su2double *Normal = new su2double[nDim]; su2double *UnitNormal = new su2double[nDim]; - su2double *Tangential = new su2double[nDim]; - /*--- Allocation of primitive gradient arrays. ---*/ + /*--- Allocation of variables necessary for viscous fluxes. ---*/ + su2double ProjGradient, ProjNormVelGrad, ProjTangVelGrad; + su2double *Tangential = new su2double[nDim]; + su2double *GradNormVel = new su2double[nDim]; + su2double *GradTangVel = new su2double[nDim]; + + /*--- Allocation of primitive gradient arrays for viscous fluxes. ---*/ su2double **Grad_Reflected = new su2double*[nPrimVarGrad]; - su2double **Grad_Prim = new su2double*[nPrimVarGrad]; - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + for (iVar = 0; iVar < nPrimVarGrad; iVar++) Grad_Reflected[iVar] = new su2double[nDim]; - Grad_Prim[iVar] = new su2double[nDim]; + + /*---------------------------------------------------------------------------------------------*/ + /*--- Preprocessing: On a symmetry-plane, the Unit-Normal is constant. Therefore a constant ---*/ + /*--- Unit-Tangential to that Unit-Normal can be prescribed. The computation ---*/ + /*--- of these vectors is done outside the loop (over all Marker-vertices). ---*/ + /*--- The "Normal" in SU2 isan Area-Normal and is most likely not constant ---*/ + /*--- on the symmetry-plane. ---*/ + /*---------------------------------------------------------------------------------------------*/ + + /*--- Normal vector for a random vertex (zero) on this marker (negate for outward convention). ---*/ + geometry->vertex[val_marker][0]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) + Normal[iDim] = -Normal[iDim]; + + /*--- Compute unit normal, to be used for unit tangential, projected velocity and velocity component gradients. ---*/ + Area = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Area += Normal[iDim]*Normal[iDim]; + Area = sqrt (Area); + + for (iDim = 0; iDim < nDim; iDim++) + UnitNormal[iDim] = -Normal[iDim]/Area; + + /*--- Preprocessing: Compute unit tangential, the direction is arbitrary as long as t*n=0 ---*/ + if (config->GetViscous()) { + switch( nDim ) { + case 2: { + Tangential[0] = -UnitNormal[1]; + Tangential[1] = UnitNormal[0]; + break; + } + case 3: { + /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ + unsigned short Largest, Arbitrary, Zero; + if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && + abs(UnitNormal[0]) >= abs(UnitNormal[2])) {Largest=0;Arbitrary=1;Zero=2;} + else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && + abs(UnitNormal[1]) >= abs(UnitNormal[2])) {Largest=1;Arbitrary=0;Zero=2;} + else {Largest=2;Arbitrary=1;Zero=0;} + + Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); + Tangential[Zero] = 0.0; + break; + } + } } - /*--- Loop over all the vertices on this boundary marker. ---*/ + /*--- Loop over all the vertices on this boundary marker. ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); @@ -6101,23 +6150,14 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Grid movement ---*/ if (config->GetGrid_Movement()) conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); - - /*--- Normal vector for this vertex (negate for outward convention). ---*/ + + /*--- Normal vector for this vertex (negate for outward convention). ---*/ geometry->vertex[val_marker][iVertex]->GetNormal(Normal); for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; conv_numerics->SetNormal(Normal); - - /*--- Compute unit normal, to be used for projected velocity and velocity component gradients. ---*/ - su2double Area = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Area += Normal[iDim]*Normal[iDim]; - Area = sqrt(Area); - - for (iDim = 0; iDim < nDim; iDim++) - UnitNormal[iDim] = -Normal[iDim]/Area; - /*--- Get current solution at this boundary node ---*/ + /*--- Get current solution at this boundary node ---*/ V_domain = node[iPoint]->GetPrimitive(); /*--- Set the reflected state based on the boundary node. Scalars are copied and @@ -6139,7 +6179,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, conv_numerics->SetPrimitive(V_domain, V_reflected); conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); - /*--- Compute the residual using an upwind scheme. ---*/ + /*--- Compute the residual using an upwind scheme. ---*/ conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); /*--- Update residual value ---*/ @@ -6150,21 +6190,22 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } - /*-------------------------------------------------------------------------------*/ - /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ - /*--- Primitive variables and their gradients. The viscous numerics ---*/ - /*--- container is filled just as the convective numerics container, ---*/ - /*--- but the primitive gradients of the reflected state have to be ---*/ - /*--- determined additionally such that symmetry at the boundary is ---*/ - /*--- enforced. Based on the Viscous_Residual routine. ---*/ - /*-------------------------------------------------------------------------------*/ if (config->GetViscous()) { + + /*-------------------------------------------------------------------------------*/ + /*--- Step 2: The viscous fluxes of the Navier-Stokes equations depend on the ---*/ + /*--- Primitive variables and their gradients. The viscous numerics ---*/ + /*--- container is filled just as the convective numerics container, ---*/ + /*--- but the primitive gradients of the reflected state have to be ---*/ + /*--- determined additionally such that symmetry at the boundary is ---*/ + /*--- enforced. Based on the Viscous_Residual routine. ---*/ + /*-------------------------------------------------------------------------------*/ /*--- Set the normal vector and the coordinates. ---*/ visc_numerics->SetCoord(geometry->node[iPoint]->GetCoord(), geometry->node[iPoint]->GetCoord()); visc_numerics->SetNormal(Normal); - /*--- Set the primitive and Secondary variables. ---*/ + /*--- Set the primitive and Secondary variables. ---*/ visc_numerics->SetPrimitive(V_domain, V_reflected); visc_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); @@ -6174,57 +6215,36 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, normal velocity in tangential direction is mirrored and the gradient of the tangential velocity in normal direction is mirrored. ---*/ - /*--- Get gradients of primitives of boundary cell ---*/ + /*--- Get gradients of primitives of boundary cell ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) for (iDim = 0; iDim < nDim; iDim++) - Grad_Prim[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); + Grad_Reflected[iVar][iDim] = node[iPoint]->GetGradient_Primitive(iVar, iDim); /*--- Reflect the gradients for all scalars including the velocity components. - The gradients of the velocity components are overriden later with the + The gradients of the velocity components are set later with the correct values: grad(V)_r = grad(V) - 2 [grad(V)*n]n, V beeing any primitive ---*/ for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + if(iVar == 0 || iVar > nDim) { // Exclude velocity component gradients - /*--- Compute projected part of the gradient in a dot product ---*/ - ProjGradient = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjGradient += Grad_Prim[iVar][iDim]*UnitNormal[iDim]; - - for (iDim = 0; iDim < nDim; iDim++) - Grad_Reflected[iVar][iDim] = Grad_Prim[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; - } - - /*--- Compute unit tangential, the direction is arbitrary as long as t*n=0. ---*/ - switch( nDim ) { - case 2: { - Tangential[0] = -UnitNormal[1]; - Tangential[1] = UnitNormal[0]; - break; - } - case 3: { - /*--- Find the largest entry index of the UnitNormal, and create Tangential vector based on that. ---*/ - unsigned short Largest, Arbitrary, Zero; - if (abs(UnitNormal[0]) >= abs(UnitNormal[1]) && abs(UnitNormal[0]) >= abs(UnitNormal[2])){Largest=0;Arbitrary=1;Zero=2;} - else if(abs(UnitNormal[1]) >= abs(UnitNormal[0]) && abs(UnitNormal[1]) >= abs(UnitNormal[2])){Largest=1;Arbitrary=0;Zero=2;} - else {Largest=2;Arbitrary=1;Zero=0;} - - Tangential[Largest] = -UnitNormal[Arbitrary]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); - Tangential[Arbitrary] = UnitNormal[Largest]/sqrt(pow(UnitNormal[Largest],2) + pow(UnitNormal[Arbitrary],2)); - Tangential[Zero] = 0.0; - break; + /*--- Compute projected part of the gradient in a dot product ---*/ + ProjGradient = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjGradient += Grad_Reflected[iVar][iDim]*UnitNormal[iDim]; + + for (iDim = 0; iDim < nDim; iDim++) + Grad_Reflected[iVar][iDim] = Grad_Reflected[iVar][iDim] - 2.0 * ProjGradient*UnitNormal[iDim]; } } /*--- Compute gradients of normal and tangential velocity: grad(v*n) = grad(v_x) n_x + grad(v_y) n_y (+ grad(v_z) n_z) grad(v*t) = grad(v_x) t_x + grad(v_y) t_y (+ grad(v_z) t_z) ---*/ - su2double GradNormVel[nDim]; - su2double GradTangVel[nDim]; for (iVar = 0; iVar < nDim; iVar++) { // counts gradient components GradNormVel[iVar] = 0.0; GradTangVel[iVar] = 0.0; for (iDim = 0; iDim < nDim; iDim++) { // counts sum with unit normal/tangential - GradNormVel[iVar] += Grad_Prim[iDim+1][iVar] * UnitNormal[iDim]; - GradTangVel[iVar] += Grad_Prim[iDim+1][iVar] * Tangential[iDim]; + GradNormVel[iVar] += Grad_Reflected[iDim+1][iVar] * UnitNormal[iDim]; + GradTangVel[iVar] += Grad_Reflected[iDim+1][iVar] * Tangential[iDim]; } } @@ -6232,19 +6252,16 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, component twice, just as done with velocity above. grad(v*n)_r = grad(v*n) - 2 {grad([v*n])*t}t grad(v*t)_r = grad(v*t) - 2 {grad([v*t])*n}n ---*/ - su2double ReflGradNormVel[nDim]; - su2double ReflGradTangVel[nDim]; - su2double ProjNormVelGrad = 0.0; - su2double ProjTangVelGrad = 0.0; - + ProjNormVelGrad = 0.0; + ProjTangVelGrad = 0.0; for (iDim = 0; iDim < nDim; iDim++) { ProjNormVelGrad += GradNormVel[iDim]*Tangential[iDim]; //grad([v*n])*t ProjTangVelGrad += GradTangVel[iDim]*UnitNormal[iDim]; //grad([v*t])*n } for (iDim = 0; iDim < nDim; iDim++) { - ReflGradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; - ReflGradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; + GradNormVel[iDim] = GradNormVel[iDim] - 2.0 * ProjNormVelGrad * Tangential[iDim]; + GradTangVel[iDim] = GradTangVel[iDim] - 2.0 * ProjTangVelGrad * UnitNormal[iDim]; } /*--- Transfer reflected gradients back into the Cartesian Coordinate system: @@ -6253,8 +6270,8 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, ( grad(v_z)_r = grad(v*n)_r n_z + grad(v*t)_r t_z ) ---*/ for (iVar = 0; iVar < nDim; iVar++) // loops over the velocity component gradients for (iDim = 0; iDim < nDim; iDim++) // loops over the entries of the above - Grad_Reflected[iVar+1][iDim] = ReflGradNormVel[iDim]*UnitNormal[iVar] + ReflGradTangVel[iDim]*Tangential[iVar]; - + Grad_Reflected[iVar+1][iDim] = GradNormVel[iDim]*UnitNormal[iVar] + GradTangVel[iDim]*Tangential[iVar]; + /*--- Set the primitive gradients of the boundary and reflected state. ---*/ visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected); @@ -6272,7 +6289,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, /*--- Jacobian contribution for implicit integration. ---*/ if (implicit) Jacobian.SubtractBlock(iPoint, iPoint, Jacobian_i); - } + } } } @@ -6280,15 +6297,13 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry *geometry, delete [] Normal; delete [] UnitNormal; delete [] Tangential; - - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - delete [] Grad_Prim[iVar]; + delete [] GradNormVel; + delete [] GradTangVel; + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) delete [] Grad_Reflected[iVar]; - } - delete [] Grad_Prim; delete [] Grad_Reflected; } - void CIncEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config) { From fa7ef0a6f0f6a32cd0512f738a1329bdc41ebe00 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Fri, 1 Mar 2019 16:53:31 +0100 Subject: [PATCH 6/8] Updated regression test values based on travis output for Sym-BC bugfix. --- .travis.yml | 2 +- SU2_CFD/src/solver_direct_turbulent.cpp | 2 +- TestCases/parallel_regression.py | 28 ++++++++++++------------- TestCases/parallel_regression_AD.py | 4 ++-- TestCases/serial_regression.py | 20 +++++++++--------- TestCases/serial_regression_AD.py | 4 ++-- 6 files changed, 30 insertions(+), 30 deletions(-) diff --git a/.travis.yml b/.travis.yml index 25d208b369af..8cc19774d0fe 100644 --- a/.travis.yml +++ b/.travis.yml @@ -82,7 +82,7 @@ install: before_script: # Get the test cases - - git clone -b develop https://github.com/su2code/TestCases.git ./TestData + - git clone -b bugfix_BC-sym-plane https://github.com/su2code/TestCases.git ./TestData - cp -R ./TestData/* ./TestCases/ # Get the tutorial cases diff --git a/SU2_CFD/src/solver_direct_turbulent.cpp b/SU2_CFD/src/solver_direct_turbulent.cpp index e933ca1aa578..b6c7990d54ce 100644 --- a/SU2_CFD/src/solver_direct_turbulent.cpp +++ b/SU2_CFD/src/solver_direct_turbulent.cpp @@ -620,7 +620,7 @@ void CTurbSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain void CTurbSolver::BC_Sym_Plane(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - /*--- Convective fluxes across symmetry plane are equal to zero. ---*/ + /*--- Convective and viscous fluxes across symmetry plane are equal to zero. ---*/ } diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 5c69142f630d..a441976dd2b7 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -90,7 +90,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-13.393796, -12.922372, 0.282557, 0.012706] #last 4 columns + oneram6.test_vals = [-10.392429, -9.840519, 0.282580, 0.012694] #last 4 columns oneram6.su2_exec = "parallel_computation.py -f" oneram6.timeout = 3200 oneram6.tol = 0.00001 @@ -139,7 +139,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 20 - flatplate.test_vals = [-4.648345, 0.813157, -0.130644, 0.024357] #last 4 columns + flatplate.test_vals = [-4.648252, 0.813253, -0.130643, 0.024357] #last 4 columns flatplate.su2_exec = "parallel_computation.py -f" flatplate.timeout = 1600 flatplate.tol = 0.00001 @@ -220,7 +220,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.146812, -6.734016, -0.176480, 0.057451] #last 4 columns + turb_flatplate.test_vals = [-4.145487, -6.734014, -0.176490, 0.057451] #last 4 columns turb_flatplate.su2_exec = "parallel_computation.py -f" turb_flatplate.timeout = 1600 turb_flatplate.tol = 0.00001 @@ -231,7 +231,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.327522, -6.564350, 0.230471, 0.155843] #last 4 columns + turb_oneram6.test_vals = [-2.327430, -6.564331, 0.230257, 0.155839] #last 4 columns turb_oneram6.su2_exec = "parallel_computation.py -f" turb_oneram6.timeout = 3200 turb_oneram6.tol = 0.00001 @@ -469,7 +469,7 @@ def main(): schubauer_klebanoff_transition.cfg_dir = "transition/Schubauer_Klebanoff" schubauer_klebanoff_transition.cfg_file = "transitional_BC_model_ConfigFile.cfg" schubauer_klebanoff_transition.test_iter = 10 - schubauer_klebanoff_transition.test_vals = [-8.219132, -14.278207, 0.000041, 0.007987] #last 4 columns + schubauer_klebanoff_transition.test_vals = [-7.994738, -14.278082, 0.000046, 0.007987] #last 4 columns schubauer_klebanoff_transition.su2_exec = "parallel_computation.py -f" schubauer_klebanoff_transition.timeout = 1600 schubauer_klebanoff_transition.tol = 0.00001 @@ -753,7 +753,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714721, -5.883008, -0.214968, 0.023783] #last 4 columns + ddes_flatplate.test_vals = [-2.714758, -5.883004, -0.215005, 0.023783] #last 4 columns ddes_flatplate.su2_exec = "parallel_computation.py -f" ddes_flatplate.timeout = 1600 ddes_flatplate.tol = 0.00001 @@ -854,7 +854,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 50 - uniform_flow.test_vals = [-0.368836, 5.156090, 0.000000, 0.000000] #last 4 columns + uniform_flow.test_vals = [-0.368877, 5.156053, 0.000000, 0.000000] #last 4 columns uniform_flow.su2_exec = "parallel_computation.py -f" uniform_flow.timeout = 1600 uniform_flow.tol = 0.000001 @@ -902,7 +902,7 @@ def main(): rotating_cylinders.cfg_dir = "sliding_interface/rotating_cylinders" rotating_cylinders.cfg_file = "rot_cylinders_WA.cfg" rotating_cylinders.test_iter = 3 - rotating_cylinders.test_vals = [-1.253430, 4.531359, 0.000000, 0.000000] #last 4 columns + rotating_cylinders.test_vals = [0.362623, 7.461154, 0.000000, 0.000000] #last 4 columns rotating_cylinders.su2_exec = "parallel_computation.py -f" rotating_cylinders.timeout = 1600 rotating_cylinders.tol = 0.00001 @@ -914,7 +914,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [-1.126767, 4.600769, 0.000000, 0.000000] #last 4 columns + supersonic_vortex_shedding.test_vals = [-1.124318, 4.605281, 0.000000, 0.000000] #last 4 columns supersonic_vortex_shedding.su2_exec = "parallel_computation.py -f" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -1001,7 +1001,7 @@ def main(): cht_incompressible.cfg_dir = "coupled_cht/incompressible" cht_incompressible.cfg_file = "config.cfg" cht_incompressible.test_iter = 10 - cht_incompressible.test_vals = [0.000000, 0.000000, -7.792549, -10248.498468] #last 4 columns + cht_incompressible.test_vals = [0.000000, 0.000000, -7.813888, -2543.238968] #last 4 columns cht_incompressible.su2_exec = "parallel_computation.py -f" cht_incompressible.timeout = 1600 cht_incompressible.tol = 0.00001 @@ -1125,7 +1125,7 @@ def main(): tutorial_inv_onera.cfg_dir = "../Tutorials/Inviscid_ONERAM6" tutorial_inv_onera.cfg_file = "inv_ONERAM6.cfg" tutorial_inv_onera.test_iter = 0 - tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.166172, 0.053116] #last 4 columns + tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.165766, 0.053239] #last 4 columns tutorial_inv_onera.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_inv_onera.timeout = 1600 tutorial_inv_onera.tol = 0.00001 @@ -1149,7 +1149,7 @@ def main(): tutorial_lam_flatplate.cfg_dir = "../Tutorials/Laminar_Flat_Plate" tutorial_lam_flatplate.cfg_file = "lam_flatplate.cfg" tutorial_lam_flatplate.test_iter = 0 - tutorial_lam_flatplate.test_vals = [-2.821818, 2.657591, -0.683901, 0.028634] #last 4 columns + tutorial_lam_flatplate.test_vals = [-2.821818, 2.657591, -0.683968, 0.028634] #last 4 columns tutorial_lam_flatplate.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_lam_flatplate.timeout = 1600 tutorial_lam_flatplate.tol = 0.00001 @@ -1161,7 +1161,7 @@ def main(): tutorial_turb_flatplate.cfg_dir = "../Tutorials/Turbulent_Flat_Plate" tutorial_turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" tutorial_turb_flatplate.test_iter = 0 - tutorial_turb_flatplate.test_vals = [-2.258584, -4.899476, -0.792617, 0.200320] #last 4 columns + tutorial_turb_flatplate.test_vals = [-2.258584, -4.899474, -0.753783, 0.200410] #last 4 columns tutorial_turb_flatplate.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_turb_flatplate.timeout = 1600 tutorial_turb_flatplate.tol = 0.00001 @@ -1185,7 +1185,7 @@ def main(): tutorial_turb_oneram6.cfg_dir = "../Tutorials/Turbulent_ONERAM6" tutorial_turb_oneram6.cfg_file = "turb_ONERAM6.cfg" tutorial_turb_oneram6.test_iter = 0 - tutorial_turb_oneram6.test_vals = [-4.499497, -11.518486, 0.391887, 0.343811] #last 4 columns + tutorial_turb_oneram6.test_vals = [-4.499497, -11.518421, 0.391293, 0.343702] #last 4 columns tutorial_turb_oneram6.su2_exec = "mpirun -np 2 SU2_CFD" tutorial_turb_oneram6.timeout = 1600 tutorial_turb_oneram6.tol = 0.00001 diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 37f71b1a47d8..c6648364d12e 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -79,7 +79,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [2.229071, 1.716910, 4.7258e+04, 0.0000e+00] #last 4 columns + discadj_arina2k.test_vals = [2.228982, 1.717042, 4.7258e+04, 0.0000e+00] #last 4 columns discadj_arina2k.su2_exec = "parallel_computation.py -f" discadj_arina2k.timeout = 8400 discadj_arina2k.tol = 0.00001 @@ -238,7 +238,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [3.162960, 0.923834, -223.148728, -3562.233908] #last 4 columns + discadj_heat.test_vals = [3.183713, 0.923840, -223.197830, -2059.808372] #last 4 columns discadj_heat.su2_exec = "parallel_computation.py -f" discadj_heat.timeout = 1600 discadj_heat.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 67189e248a1c..136002925bfa 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -88,7 +88,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-13.395738, -12.930653, 0.282557, 0.012706] #last 4 columns + oneram6.test_vals = [-10.384532, -9.835738, 0.282580, 0.012694] #last 4 columns oneram6.su2_exec = "SU2_CFD" oneram6.timeout = 9600 oneram6.tol = 0.00001 @@ -137,7 +137,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 20 - flatplate.test_vals = [-4.680896, 0.781111, -0.135957, 0.022978] #last 4 columns + flatplate.test_vals = [-4.680777, 0.781234, -0.135957, 0.022977] #last 4 columns flatplate.su2_exec = "SU2_CFD" flatplate.timeout = 1600 flatplate.tol = 0.00001 @@ -218,7 +218,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.158303, -6.737135, -0.176244, 0.057446] #last 4 columns + turb_flatplate.test_vals = [-4.157169, -6.737133, -0.176253, 0.057446] #last 4 columns turb_flatplate.su2_exec = "SU2_CFD" turb_flatplate.timeout = 1600 turb_flatplate.tol = 0.00001 @@ -229,7 +229,7 @@ def main(): turb_oneram6.cfg_dir = "rans/oneram6" turb_oneram6.cfg_file = "turb_ONERAM6.cfg" turb_oneram6.test_iter = 10 - turb_oneram6.test_vals = [-2.327523, -6.564349, 0.230471, 0.155843]#last 4 columns + turb_oneram6.test_vals = [-2.327431, -6.564331, 0.230257, 0.155839]#last 4 columns turb_oneram6.su2_exec = "SU2_CFD" turb_oneram6.timeout = 3200 turb_oneram6.tol = 0.00001 @@ -467,7 +467,7 @@ def main(): schubauer_klebanoff_transition.cfg_dir = "transition/Schubauer_Klebanoff" schubauer_klebanoff_transition.cfg_file = "transitional_BC_model_ConfigFile.cfg" schubauer_klebanoff_transition.test_iter = 10 - schubauer_klebanoff_transition.test_vals = [-8.287490, -14.278189, 0.000050, 0.007986] #last 4 columns + schubauer_klebanoff_transition.test_vals = [-8.029756, -14.278066, 0.000053, 0.007986] #last 4 columns schubauer_klebanoff_transition.su2_exec = "SU2_CFD" schubauer_klebanoff_transition.timeout = 1600 schubauer_klebanoff_transition.tol = 0.00001 @@ -751,7 +751,7 @@ def main(): ddes_flatplate.cfg_dir = "ddes/flatplate" ddes_flatplate.cfg_file = "ddes_flatplate.cfg" ddes_flatplate.test_iter = 10 - ddes_flatplate.test_vals = [-2.714721, -5.883008, -0.214968, 0.023783] #last 4 columns + ddes_flatplate.test_vals = [-2.714758, -5.883004, -0.215005, 0.023783] #last 4 columns ddes_flatplate.su2_exec = "SU2_CFD" ddes_flatplate.timeout = 1600 ddes_flatplate.tol = 0.00001 @@ -865,7 +865,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 50 - uniform_flow.test_vals = [-0.368836, 5.156090, 0.000000, 0.000000] #last 4 columns + uniform_flow.test_vals = [-0.368877, 5.156053, 0.000000, 0.000000] #last 4 columns uniform_flow.su2_exec = "SU2_CFD" uniform_flow.timeout = 1600 uniform_flow.tol = 0.000001 @@ -913,7 +913,7 @@ def main(): rotating_cylinders.cfg_dir = "sliding_interface/rotating_cylinders" rotating_cylinders.cfg_file = "rot_cylinders_WA.cfg" rotating_cylinders.test_iter = 3 - rotating_cylinders.test_vals = [-1.253498, 4.531302, 0.000000, 0.000000] #last 4 columns + rotating_cylinders.test_vals = [-1.254672, 4.530738, 0.000000, 0.000000] #last 4 columns rotating_cylinders.su2_exec = "SU2_CFD" rotating_cylinders.timeout = 1600 rotating_cylinders.tol = 0.00001 @@ -925,7 +925,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [-1.128085, 4.600597, 0.000000, 0.000000] #last 4 columns + supersonic_vortex_shedding.test_vals = [-1.130591, 4.595041, 0.000000, 0.000000] #last 4 columns supersonic_vortex_shedding.su2_exec = "SU2_CFD" supersonic_vortex_shedding.timeout = 1600 supersonic_vortex_shedding.tol = 0.00001 @@ -1023,7 +1023,7 @@ def main(): cht_incompressible.cfg_dir = "coupled_cht/incompressible" cht_incompressible.cfg_file = "config.cfg" cht_incompressible.test_iter = 10 - cht_incompressible.test_vals = [0.000000, 0.000000, -7.685301, -12947.783696] #last 4 columns + cht_incompressible.test_vals = [0.000000, 0.000000, -8.530925, -3091.634678] #last 4 columns cht_incompressible.su2_exec = "SU2_CFD" cht_incompressible.timeout = 1600 cht_incompressible.tol = 0.0001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index ea476886aaa0..1229c4a402dc 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -79,7 +79,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [-0.774805, -0.801209, 3.1979e+02, 0.0000e+00] #last 4 columns + discadj_arina2k.test_vals = [2.206814, 1.674391, 4.7250e+04, 0.0000e+00] #last 4 columns discadj_arina2k.su2_exec = "SU2_CFD_AD" discadj_arina2k.timeout = 8400 discadj_arina2k.tol = 0.00001 @@ -223,7 +223,7 @@ def main(): discadj_heat.cfg_dir = "disc_adj_heat" discadj_heat.cfg_file = "disc_adj_heat.cfg" discadj_heat.test_iter = 10 - discadj_heat.test_vals = [3.176483, 1.144873, -1040.512028, -3277.663739] #last 4 columns + discadj_heat.test_vals = [3.139355, 1.144919, -1040.637744, -2464.935518] #last 4 columns discadj_heat.su2_exec = "SU2_CFD_AD" discadj_heat.timeout = 1600 discadj_heat.tol = 0.00001 From ee8b26fd786727e50d78f8f059a918b6bd005d70 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Sun, 3 Mar 2019 22:36:12 +0100 Subject: [PATCH 7/8] Update regression test again for sym-plane bugfix. --- .travis.yml | 4 ++-- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression_AD.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8cc19774d0fe..2a85523ac883 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,11 +12,11 @@ compiler: notifications: email: recipients: - - tobias.kattmann@de.bosch.com + - su2code-dev@lists.stanford.edu branches: only: - - bugfix_BC-sym-plane + - develop python: - 2.7 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index a441976dd2b7..fe820aab08d5 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -902,7 +902,7 @@ def main(): rotating_cylinders.cfg_dir = "sliding_interface/rotating_cylinders" rotating_cylinders.cfg_file = "rot_cylinders_WA.cfg" rotating_cylinders.test_iter = 3 - rotating_cylinders.test_vals = [0.362623, 7.461154, 0.000000, 0.000000] #last 4 columns + rotating_cylinders.test_vals = [1.219987, 7.729743, 0.000000, 0.000000] #last 4 columns rotating_cylinders.su2_exec = "parallel_computation.py -f" rotating_cylinders.timeout = 1600 rotating_cylinders.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 1229c4a402dc..394518c89357 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -79,7 +79,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [2.206814, 1.674391, 4.7250e+04, 0.0000e+00] #last 4 columns + discadj_arina2k.test_vals = [-0.776022, -0.795092, 319.800000, 0.000000] #last 4 columns discadj_arina2k.su2_exec = "SU2_CFD_AD" discadj_arina2k.timeout = 8400 discadj_arina2k.tol = 0.00001 From c2e09d9f01cb892daec613a57b84a7b44bda839a Mon Sep 17 00:00:00 2001 From: "Thomas D. Economon" Date: Wed, 13 Mar 2019 15:15:08 -0700 Subject: [PATCH 8/8] Reverting travis. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 2a85523ac883..b974c79b1279 100644 --- a/.travis.yml +++ b/.travis.yml @@ -82,7 +82,7 @@ install: before_script: # Get the test cases - - git clone -b bugfix_BC-sym-plane https://github.com/su2code/TestCases.git ./TestData + - git clone -b develop https://github.com/su2code/TestCases.git ./TestData - cp -R ./TestData/* ./TestCases/ # Get the tutorial cases