Removal of code related to the Split Velocity Method

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.

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.

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.

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);
       }

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, Gust[3] = {0.0}, NewGridVel[3] = {0.0};
+  const su2double* GridVel = nullptr;
 
   su2double Physical_dt = config->GetDelta_UnstTime();
   unsigned long TimeIter = config->GetTimeIter();
@@ -329,10 +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]();
-  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) {
 	  cout << endl << "Setting up a wind gust type " << Gust_Type << " with amplitude of " << gust_amp << " in direction " << GustDir << endl;
@@ -370,10 +361,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,23 +385,13 @@ 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;
 
           case ONE_M_COSINE:
             // 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 +426,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 ---*/
@@ -475,10 +436,6 @@ void CFluidIteration::SetWind_GustField(CConfig* config, CGeometry** geometry, C
       }
     }
   }
-
-  delete[] Gust;
-  delete[] GustDer;
-  delete[] NewGridVel;
 }
 
 void CFluidIteration::InitializeVortexDistribution(unsigned long& nVortex, vector<su2double>& x0, vector<su2double>& y0,

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) :

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];

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()) {