From 0a0dc9e5dc4bcf99aa1356010991fcfcc7fed476 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Arne=20Vo=C3=9F?= Date: Fri, 9 Jun 2023 15:36:30 +0200 Subject: [PATCH 1/2] First step in removing the unused code of the Split Velocity Method. --- SU2_CFD/include/numerics/CNumerics.hpp | 24 +----- SU2_CFD/include/variables/CEulerVariable.hpp | 31 -------- SU2_CFD/include/variables/CVariable.hpp | 28 ------- SU2_CFD/src/drivers/CDriver.cpp | 3 - SU2_CFD/src/iteration/CFluidIteration.cpp | 41 +--------- SU2_CFD/src/numerics/flow/flow_sources.cpp | 84 -------------------- SU2_CFD/src/solvers/CEulerSolver.cpp | 32 -------- SU2_CFD/src/variables/CEulerVariable.cpp | 1 - 8 files changed, 3 insertions(+), 241 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index e8682cc8da58..c2068ef311b5 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -171,9 +171,7 @@ class CNumerics { su2double vel2_inf; /*!< \brief value of the square of freestream speed. */ const su2double *WindGust_i, /*!< \brief Wind gust at point i. */ - *WindGust_j, /*!< \brief Wind gust at point j. */ - *WindGustDer_i, /*!< \brief Wind gust derivatives at point i. */ - *WindGustDer_j; /*!< \brief Wind gust derivatives at point j. */ + *WindGust_j; /*!< \brief Wind gust at point j. */ const su2double *Vorticity_i, *Vorticity_j; /*!< \brief Vorticity. */ su2double StrainMag_i, StrainMag_j; /*!< \brief Strain rate magnitude. */ su2double Dissipation_i, Dissipation_j; /*!< \brief Dissipation. */ @@ -860,26 +858,6 @@ class CNumerics { GridVel_j = val_gridvel_j; } - /*! - * \brief Set the wind gust value. - * \param[in] val_windgust_i - Wind gust of the point i. - * \param[in] val_windgust_j - Wind gust of the point j. - */ - inline void SetWindGust(const su2double *val_windgust_i, const su2double *val_windgust_j) { - WindGust_i = val_windgust_i; - WindGust_j = val_windgust_j; - } - - /*! - * \brief Set the wind gust derivatives values. - * \param[in] val_windgust_i - Wind gust derivatives of the point i. - * \param[in] val_windgust_j - Wind gust derivatives of the point j. - */ - inline void SetWindGustDer(const su2double *val_windgustder_i, const su2double *val_windgustder_j) { - WindGustDer_i = val_windgustder_i; - WindGustDer_j = val_windgustder_j; - } - /*! * \brief Set the value of the pressure. * \param[in] val_pressure_i - Value of the pressure at point i. diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 44f07ffa0a2e..43641a27199f 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -72,7 +72,6 @@ class CEulerVariable : public CFlowVariable { MatrixType Secondary; MatrixType WindGust; /*! < \brief Wind gust value */ - MatrixType WindGustDer; /*! < \brief Wind gust derivatives value */ public: /*! @@ -285,36 +284,6 @@ class CEulerVariable : public CFlowVariable { for (unsigned long iDim = 0; iDim < nDim; iDim++) Res_TruncError(iPoint,iDim+1) = 0.0; } - /*! - * \brief Get the value of the wind gust - * \return Value of the wind gust - */ - inline su2double* GetWindGust(unsigned long iPoint) final { return WindGust[iPoint]; } - - /*! - * \brief Set the value of the wind gust - * \param[in] Value of the wind gust - */ - inline void SetWindGust(unsigned long iPoint, const su2double* val_WindGust) final { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - WindGust(iPoint,iDim) = val_WindGust[iDim]; - } - - /*! - * \brief Get the value of the derivatives of the wind gust - * \return Value of the derivatives of the wind gust - */ - inline su2double* GetWindGustDer(unsigned long iPoint) final { return WindGustDer[iPoint]; } - - /*! - * \brief Set the value of the derivatives of the wind gust - * \param[in] Value of the derivatives of the wind gust - */ - inline void SetWindGustDer(unsigned long iPoint, const su2double* val_WindGustDer) final { - for (unsigned long iDim = 0; iDim < nDim+1; iDim++) - WindGustDer(iPoint,iDim) = val_WindGustDer[iDim]; - } - /*! * \brief Specify a vector to set the velocity components of the solution. Multiplied by density for compressible cases. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 5167ef953f69..a1c89d5a4145 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -794,34 +794,6 @@ class CVariable { inline MatrixType& GetSolution_Min() { return Solution_Min; } inline const MatrixType& GetSolution_Min() const { return Solution_Min; } - /*! - * \brief Get the value of the wind gust - * \param[in] iPoint - Point index. - * \return Value of the wind gust - */ - inline virtual su2double* GetWindGust(unsigned long iPoint) { return nullptr; } - - /*! - * \brief Set the value of the wind gust - * \param[in] iPoint - Point index. - * \param[in] val_WindGust - Value of the wind gust - */ - inline virtual void SetWindGust(unsigned long iPoint, const su2double* val_WindGust) {} - - /*! - * \brief Get the value of the derivatives of the wind gust - * \param[in] iPoint - Point index. - * \return Value of the derivatives of the wind gust - */ - inline virtual su2double* GetWindGustDer(unsigned long iPoint) { return nullptr;} - - /*! - * \brief Set the value of the derivatives of the wind gust - * \param[in] iPoint - Point index. - * \param[in] val_WindGust - Value of the derivatives of the wind gust - */ - inline virtual void SetWindGustDer(unsigned long iPoint, const su2double* val_WindGust) {} - /*! * \brief Set the value of the time step. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 1eee1af5e852..5b85ead935eb 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1873,9 +1873,6 @@ void CDriver::InitializeNumerics(CConfig *config, CGeometry **geometry, CSolver else if (config->GetGravityForce() == YES) { numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGravity(nDim, nVar_Flow, config); } - else if (config->GetWind_Gust() == YES) { - numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceWindGust(nDim, nVar_Flow, config); - } else { numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceNothing(nDim, nVar_Flow, config); } diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index da63e7812a51..a5948656b550 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -297,11 +297,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C // described in the NASA TM–2012-217771 - Development, Verification and Use of Gust Modeling in the NASA Computational // Fluid Dynamics Code FUN3D the desired gust is prescribed as the negative of the grid velocity. - // If a source term is included to account for the gust field, the method is described by Jones et al. as the Split - // Velocity Method in Simulation of Airfoil Gust Responses Using Prescribed Velocities. In this routine the gust - // derivatives needed for the source term are calculated when applicable. If the gust derivatives are zero the source - // term is also zero. The source term itself is implemented in the class CSourceWindGust - unsigned short iDim, nDim = geometry[MESH_0]->GetnDim(); /*--- Gust Parameters from config ---*/ @@ -318,8 +313,8 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C unsigned long iPoint; unsigned short iMGlevel, nMGlevel = config->GetnMGLevels(); - su2double x, y, x_gust, dgust_dx, dgust_dy, dgust_dz, dgust_dt; - su2double *Gust, *GridVel, *NewGridVel, *GustDer; + su2double x, y, x_gust; + su2double *Gust, *GridVel, *NewGridVel; su2double Physical_dt = config->GetDelta_UnstTime(); unsigned long TimeIter = config->GetTimeIter(); @@ -331,7 +326,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C Gust = new su2double[nDim](); NewGridVel = new su2double[nDim](); - GustDer = new su2double[nDim+1](); // Print some information to check that we are doing the right thing. Not sure how to convert the index back to a string... if (rank == MASTER_NODE) { @@ -370,10 +364,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C for (iDim = 0; iDim < nDim; iDim++) { Gust[iDim] = 0.0; } - dgust_dx = 0.0; - dgust_dy = 0.0; - dgust_dz = 0.0; - dgust_dt = 0.0; /*--- Begin applying the gust ---*/ @@ -398,11 +388,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C // Check if we are in the region where the gust is active if (x_gust > 0 && x_gust < n) { Gust[GustDir] = gust_amp * (sin(2 * PI_NUMBER * x_gust)); - - // Gust derivatives - // dgust_dx = gust_amp*2*PI_NUMBER*(cos(2*PI_NUMBER*x_gust))/L; - // dgust_dy = 0; - // dgust_dt = gust_amp*2*PI_NUMBER*(cos(2*PI_NUMBER*x_gust))*(-Uinf)/L; } break; @@ -410,11 +395,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C // Check if we are in the region where the gust is active if (x_gust > 0 && x_gust < n) { Gust[GustDir] = gust_amp * 0.5 * (1 - cos(2 * PI_NUMBER * x_gust)); - - // Gust derivatives - // dgust_dx = gust_amp*2*PI_NUMBER*(sin(2*PI_NUMBER*x_gust))/L; - // dgust_dy = 0; - // dgust_dt = gust_amp*2*PI_NUMBER*(sin(2*PI_NUMBER*x_gust))*(-Uinf)/L; } break; @@ -449,22 +429,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C } } - /*--- Set the Wind Gust, Wind Gust Derivatives and the Grid Velocities ---*/ - if (nDim == 2) { - GustDer[0] = dgust_dx; - GustDer[1] = dgust_dy; - GustDer[2] = dgust_dt; - } - else { - GustDer[0] = dgust_dx; - GustDer[1] = dgust_dy; - GustDer[2] = dgust_dz; - GustDer[3] = dgust_dt; - } - // I think we don't need to set any source terms because they depend on the derivatives, which are zero in all cases from above. - solver[iMGlevel][FLOW_SOL]->GetNodes()->SetWindGust(iPoint, Gust); - solver[iMGlevel][FLOW_SOL]->GetNodes()->SetWindGustDer(iPoint, GustDer); - GridVel = geometry[iMGlevel]->nodes->GetGridVel(iPoint); /*--- Store new grid velocity ---*/ @@ -477,7 +441,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C } delete[] Gust; - delete[] GustDer; delete[] NewGridVel; } diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index cd8e8d1d4ad4..64f8145cf150 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -715,90 +715,6 @@ CNumerics::ResidualType<> CSourceVorticityConfinement::ComputeResidual(const CCo return ResidualType<>(residual, jacobian, nullptr); } -CSourceWindGust::CSourceWindGust(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config) : - CSourceBase_Flow(val_nDim, val_nVar, config) { } - -CNumerics::ResidualType<> CSourceWindGust::ComputeResidual(const CConfig* config) { - - su2double u_gust, v_gust, du_gust_dx, du_gust_dy, du_gust_dt, dv_gust_dx, dv_gust_dy, dv_gust_dt; - su2double smx, smy, se, rho, u, v, p; - unsigned short GustDir = config->GetGust_Dir(); //Gust direction - - u_gust = WindGust_i[0]; - v_gust = WindGust_i[1]; -// w_gust = WindGust_i[2]; - - if (GustDir == X_DIR) { - du_gust_dx = WindGustDer_i[0]; - du_gust_dy = WindGustDer_i[1]; - //du_gust_dz = WindGustDer_i[2]; - du_gust_dt = WindGustDer_i[2]; - - dv_gust_dx = 0.0; - dv_gust_dy = 0.0; - //dv_gust_dz = 0.0; - dv_gust_dt = 0.0; - - //dw_gust_dx = 0.0; - //dw_gust_dy = 0.0; - //dw_gust_dz = 0.0; - //dw_gust_dt = 0.0; - } else { - du_gust_dx = 0.0; - du_gust_dy = 0.0; - //du_gust_dz = 0.0; - du_gust_dt = 0.0; - dv_gust_dx = WindGustDer_i[0]; - dv_gust_dy = WindGustDer_i[1]; - //dv_gust_dz = WindGustDer_i[2] - dv_gust_dt = WindGustDer_i[2]; - - //dw_gust_dx = 0.0; - //dw_gust_dy = 0.0; - //dw_gust_dz = 0.0; - //dw_gust_dt = 0.0; - // - - } - - /*--- Primitive variables at point i ---*/ - u = V_i[1]; - v = V_i[2]; - // w = V_i[3] - - p = V_i[nDim+1]; - rho = V_i[nDim+2]; - - /*--- Source terms ---*/ - smx = rho*(du_gust_dt + (u+u_gust)*du_gust_dx + (v+v_gust)*du_gust_dy); - smy = rho*(dv_gust_dt + (u+u_gust)*dv_gust_dx + (v+v_gust)*dv_gust_dy); - //smz = rho*(dw_gust_dt + (u+u_gust)*dw_gust_dx + (v+v_gust)*dw_gust_dy) + (w+w_gust)*dw_gust_dz; - - se = u*smx + v*smy + p*(du_gust_dx + dv_gust_dy); - //se = u*smx + v*smy + w*smz + p*(du_gust_dx + dv_gust_dy + dw_gust_dz); - - if (nDim == 2) { - residual[0] = 0.0; - residual[1] = smx*Volume; - residual[2] = smy*Volume; - //residual[3] = smz*Volume; - residual[3] = se*Volume; - } - else { - residual[0] = 0.0; - residual[1] = 0.0; - residual[2] = 0.0; - residual[3] = 0.0; - residual[4] = 0.0; - } - - /*--- For now the source term Jacobian is just set to zero ---*/ - //bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - - return ResidualType<>(residual, jacobian, nullptr); -} - - CSourceIncStreamwise_Periodic::CSourceIncStreamwise_Periodic(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 5db213329f9b..a8dedcba6d7f 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2010,7 +2010,6 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain const bool axisymmetric = config->GetAxisymmetric(); const bool gravity = (config->GetGravityForce() == YES); const bool harmonic_balance = (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE); - const bool windgust = config->GetWind_Gust(); const bool body_force = config->GetBody_Force(); const bool vorticity_confinement = config->GetVorticityConfinement(); const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || @@ -2180,37 +2179,6 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain END_SU2_OMP_FOR } - if (windgust) { - - /*--- Loop over all points ---*/ - SU2_OMP_FOR_DYN(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Load the wind gust ---*/ - numerics->SetWindGust(nodes->GetWindGust(iPoint), nodes->GetWindGust(iPoint)); - - /*--- Load the wind gust derivatives ---*/ - numerics->SetWindGustDer(nodes->GetWindGustDer(iPoint), nodes->GetWindGustDer(iPoint)); - - /*--- Load the primitive variables ---*/ - numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint)); - - /*--- Load the volume of the dual mesh cell ---*/ - numerics->SetVolume(geometry->nodes->GetVolume(iPoint)); - - /*--- Compute the rotating frame source residual ---*/ - auto residual = numerics->ComputeResidual(config); - - /*--- Add the source residual to the total ---*/ - LinSysRes.AddBlock(iPoint, residual); - - /*--- Add the implicit Jacobian contribution ---*/ - if (implicit) Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - - } - END_SU2_OMP_FOR - } - if (vorticity_confinement) { CNumerics* second_numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num()*MAX_TERMS]; diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 4a666042c961..cf3432ebd071 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -70,7 +70,6 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 if (config->GetWind_Gust()) { WindGust.resize(nPoint,nDim); - WindGustDer.resize(nPoint,nDim+1); } if (config->GetVorticityConfinement()) { From 33ba4183870633562beda47a950514a88540a475 Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Sun, 11 Jun 2023 22:12:02 +0100 Subject: [PATCH 2/2] Apply suggestions from code review --- SU2_CFD/src/iteration/CFluidIteration.cpp | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index a5948656b550..cde93a085bbb 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -313,8 +313,8 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C unsigned long iPoint; unsigned short iMGlevel, nMGlevel = config->GetnMGLevels(); - su2double x, y, x_gust; - su2double *Gust, *GridVel, *NewGridVel; + su2double x, y, x_gust, Gust[3] = {0.0}, NewGridVel[3] = {0.0}; + const su2double* GridVel = nullptr; su2double Physical_dt = config->GetDelta_UnstTime(); unsigned long TimeIter = config->GetTimeIter(); @@ -324,9 +324,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C su2double Uinf = solver[MESH_0][FLOW_SOL]->GetVelocity_Inf(0); // Assumption gust moves at infinity velocity - Gust = new su2double[nDim](); - NewGridVel = new su2double[nDim](); - // Print some information to check that we are doing the right thing. Not sure how to convert the index back to a string... if (rank == MASTER_NODE) { cout << endl << "Setting up a wind gust type " << Gust_Type << " with amplitude of " << gust_amp << " in direction " << GustDir << endl; @@ -439,9 +436,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C } } } - - delete[] Gust; - delete[] NewGridVel; } void CFluidIteration::InitializeVortexDistribution(unsigned long& nVortex, vector& x0, vector& y0,