fix #1175

Author: pcarruscag

SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp

@@ -92,6 +92,7 @@ class CFEM_DG_EulerSolver : public CSolver {
   AllBound_CEff_Inv;  /*!< \brief Total efficiency (Cl/Cd) (inviscid contribution) for all the boundaries. */
 
   su2double
+  AeroCoeffForceRef,   /*!< \brief Reference force for coefficients */
   Total_CL,      /*!< \brief Total lift coefficient for all the boundaries. */
   Total_CD,        /*!< \brief Total drag coefficient for all the boundaries. */
   Total_CSF,       /*!< \brief Total sideforce coefficient for all the boundaries. */
@@ -1127,6 +1128,11 @@ class CFEM_DG_EulerSolver : public CSolver {
    */
   inline void SetTotal_CL(su2double val_Total_CL) final { Total_CL = val_Total_CL; }
 
+  /*!
+   * \brief Get the reference force used to compute CL, CD, etc.
+   */
+  inline su2double GetAeroCoeffsReferenceForce() const final { return AeroCoeffForceRef; }
+
   /*!
    * \brief Provide the total (inviscid + viscous) non dimensional lift coefficient.
    * \return Value of the lift coefficient (inviscid + viscous contribution).

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -125,6 +125,8 @@ class CFVMFlowSolverBase : public CSolver {
   AeroCoeffsArray SurfaceCoeff; /*!< \brief Totals for each monitoring surface. */
   AeroCoeffs TotalCoeff;        /*!< \brief Totals for all boundaries. */
 
+  su2double AeroCoeffForceRef = 1.0;    /*!< \brief Reference force for aerodynamic coefficients. */
+
   su2double InverseDesign = 0.0;        /*!< \brief Inverse design functional for each boundary. */
   su2double Total_ComboObj = 0.0;       /*!< \brief Total 'combo' objective for all monitored boundaries */
   su2double Total_Custom_ObjFunc = 0.0; /*!< \brief Total custom objective function for all the boundaries. */
@@ -1559,6 +1561,11 @@ class CFVMFlowSolverBase : public CSolver {
    */
   inline su2double GetTotal_CEff() const final { return TotalCoeff.CEff; }
 
+  /*!
+   * \brief Get the reference force used to compute CL, CD, etc.
+   */
+  inline su2double GetAeroCoeffsReferenceForce() const final { return AeroCoeffForceRef; }
+
   /*!
    * \brief Provide the total (inviscid + viscous) non dimensional lift coefficient.
    * \return Value of the lift coefficient (inviscid + viscous contribution).

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1632,7 +1632,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Pressure_Forces(const CGeometry* geometr
     }
   }
 
-  factor = 1.0 / (0.5 * RefDensity * RefArea * RefVel2);
+  AeroCoeffForceRef = 0.5 * RefDensity * RefArea * RefVel2;
+  factor = 1.0 / AeroCoeffForceRef;
 
   /*--- Reference pressure is always the far-field value. ---*/
 
@@ -1938,57 +1939,18 @@ template <class V, ENUM_REGIME FlowRegime>
 void CFVMFlowSolverBase<V, FlowRegime>::Momentum_Forces(const CGeometry* geometry, const CConfig* config) {
   unsigned long iVertex, iPoint;
   unsigned short iDim, iMarker, Boundary, Monitoring, iMarker_Monitoring;
-  su2double factor, RefVel2 = 0.0, RefTemp, RefDensity = 0.0, Mach2Vel, Mach_Motion, MassFlow, Density;
+  su2double MassFlow, Density;
   const su2double *Normal = nullptr, *Coord = nullptr;
   string Marker_Tag, Monitoring_Tag;
   su2double AxiFactor;
 
   su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
   su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
-  su2double RefArea = config->GetRefArea();
   su2double RefLength = config->GetRefLength();
-  su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin = config->GetRefOriginMoment(0);
   bool axisymmetric = config->GetAxisymmetric();
 
-  /// TODO: Move these ifs to specialized functions.
-
-  if (FlowRegime == COMPRESSIBLE) {
-    /*--- Evaluate reference values for non-dimensionalization.
-     For dynamic meshes, use the motion Mach number as a reference value
-     for computing the force coefficients. Otherwise, use the freestream values,
-     which is the standard convention. ---*/
-
-    RefTemp = Temperature_Inf;
-    RefDensity = Density_Inf;
-    if (dynamic_grid) {
-      Mach2Vel = sqrt(Gamma * Gas_Constant * RefTemp);
-      Mach_Motion = config->GetMach_Motion();
-      RefVel2 = (Mach_Motion * Mach2Vel) * (Mach_Motion * Mach2Vel);
-    } else {
-      RefVel2 = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) RefVel2 += Velocity_Inf[iDim] * Velocity_Inf[iDim];
-    }
-  }
-
-  if (FlowRegime == INCOMPRESSIBLE) {
-    /*--- Evaluate reference values for non-dimensionalization.
-     For dimensional or non-dim based on initial values, use
-     the far-field state (inf). For a custom non-dim based
-     on user-provided reference values, use the ref values
-     to compute the forces. ---*/
-
-    if ((config->GetRef_Inc_NonDim() == DIMENSIONAL) || (config->GetRef_Inc_NonDim() == INITIAL_VALUES)) {
-      RefDensity = Density_Inf;
-      RefVel2 = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) RefVel2 += Velocity_Inf[iDim] * Velocity_Inf[iDim];
-    } else if (config->GetRef_Inc_NonDim() == REFERENCE_VALUES) {
-      RefDensity = config->GetInc_Density_Ref();
-      RefVel2 = config->GetInc_Velocity_Ref() * config->GetInc_Velocity_Ref();
-    }
-  }
-
-  factor = 1.0 / (0.5 * RefDensity * RefArea * RefVel2);
+  const su2double factor = 1.0 / AeroCoeffForceRef;
 
   /*-- Variables initialization ---*/
 
@@ -2263,8 +2225,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0, dTn, dTven,
-            RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
+  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefVel2 = 0.0, dTn, dTven,
+            RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
@@ -2273,7 +2235,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
   su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
-  su2double RefArea = config->GetRefArea();
   su2double RefLength = config->GetRefLength();
   su2double RefHeatFlux = config->GetHeat_Flux_Ref();
   su2double Gas_Constant = config->GetGas_ConstantND();
@@ -2294,44 +2255,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     VEL_INDEX     = nSpecies+2;
   }
 
-  /// TODO: Move these ifs to specialized functions.
-
-  if (FlowRegime == COMPRESSIBLE) {
-    /*--- Evaluate reference values for non-dimensionalization.
-     For dynamic meshes, use the motion Mach number as a reference value
-     for computing the force coefficients. Otherwise, use the freestream values,
-     which is the standard convention. ---*/
-
-    RefTemp = Temperature_Inf;
-    RefDensity = Density_Inf;
-    if (dynamic_grid) {
-      Mach2Vel = sqrt(Gamma * Gas_Constant * RefTemp);
-      Mach_Motion = config->GetMach_Motion();
-      RefVel2 = (Mach_Motion * Mach2Vel) * (Mach_Motion * Mach2Vel);
-    } else {
-      RefVel2 = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) RefVel2 += Velocity_Inf[iDim] * Velocity_Inf[iDim];
-    }
-  }
-
-  if (FlowRegime == INCOMPRESSIBLE) {
-    /*--- Evaluate reference values for non-dimensionalization.
-     For dimensional or non-dim based on initial values, use
-     the far-field state (inf). For a custom non-dim based
-     on user-provided reference values, use the ref values
-     to compute the forces. ---*/
-
-    if ((config->GetRef_Inc_NonDim() == DIMENSIONAL) || (config->GetRef_Inc_NonDim() == INITIAL_VALUES)) {
-      RefDensity = Density_Inf;
-      RefVel2 = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) RefVel2 += Velocity_Inf[iDim] * Velocity_Inf[iDim];
-    } else if (config->GetRef_Inc_NonDim() == REFERENCE_VALUES) {
-      RefDensity = config->GetInc_Density_Ref();
-      RefVel2 = config->GetInc_Velocity_Ref() * config->GetInc_Velocity_Ref();
-    }
-  }
-
-  const su2double factor = 1.0 / (0.5 * RefDensity * RefArea * RefVel2);
+  const su2double factor = 1.0 / AeroCoeffForceRef;
 
   /*--- Variables initialization ---*/
 
@@ -2452,8 +2376,6 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 
-
-
       /// TODO: Move these ifs to specialized functions.
       if (!nemo){
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -2516,6 +2516,11 @@ class CSolver {
    */
   inline virtual void SetTotal_CNearFieldOF(su2double val_cnearfieldpress) { }
 
+  /*!
+   * \brief Get the reference force used to compute CL, CD, etc.
+   */
+  inline virtual su2double GetAeroCoeffsReferenceForce() const { return 0; }
+
   /*!
    * \brief A virtual member.
    * \return Value of the lift coefficient (inviscid + viscous contribution).

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1022,7 +1022,7 @@ void CFlowOutput::WriteMetaData(CConfig *config){
 
 void CFlowOutput::WriteForcesBreakdown(CConfig *config, CGeometry *geometry, CSolver **solver_container){
 
-  unsigned short iDim, iMarker_Monitoring;
+  unsigned short iMarker_Monitoring;
 
   const bool compressible    = (config->GetKind_Regime() == COMPRESSIBLE);
   const bool incompressible  = (config->GetKind_Regime() == INCOMPRESSIBLE);
@@ -2144,26 +2144,8 @@ void CFlowOutput::WriteForcesBreakdown(CConfig *config, CGeometry *geometry, CSo
 
     /*--- Reference area and force factors. ---*/
 
-    su2double RefDensity, RefArea, RefVel, Factor, Ref;
-    RefArea     = config->GetRefArea();
-    if (compressible) {
-      RefDensity  = solver_container[FLOW_SOL]->GetDensity_Inf();
-      RefVel = solver_container[FLOW_SOL]->GetModVelocity_Inf();
-    } else {
-      if ((config->GetRef_Inc_NonDim() == DIMENSIONAL) ||
-          (config->GetRef_Inc_NonDim() == INITIAL_VALUES)) {
-        RefDensity  = solver_container[FLOW_SOL]->GetDensity_Inf();
-        RefVel = 0.0;
-        for (iDim = 0; iDim < nDim; iDim++)
-          RefVel  += solver_container[FLOW_SOL]->GetVelocity_Inf(iDim)*solver_container[FLOW_SOL]->GetVelocity_Inf(iDim);
-        RefVel = sqrt(RefVel);
-      } else {
-        RefDensity = config->GetInc_Density_Ref();
-        RefVel    = config->GetInc_Velocity_Ref();
-      }
-    }
-    Factor = (0.5*RefDensity*RefArea*RefVel*RefVel);
-    Ref = config->GetDensity_Ref() * config->GetVelocity_Ref() * config->GetVelocity_Ref() * 1.0 * 1.0;
+    const su2double Factor = solver_container[FLOW_SOL]->GetAeroCoeffsReferenceForce();
+    const su2double Ref = config->GetDensity_Ref() * pow(config->GetVelocity_Ref(),2);
 
     Breakdown_file << "NOTE: Multiply forces by the non-dimensional factor: " << Factor << ", and the reference factor: " << Ref  << "\n";
     Breakdown_file << "to obtain the dimensional force."  << "\n" << "\n";

SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp

@@ -6752,7 +6752,8 @@ void CFEM_DG_EulerSolver::Pressure_Forces(const CGeometry* geometry, const CConf
       RefVel2 += Velocity_Inf[iDim]*Velocity_Inf[iDim];
   }
 
-  const su2double factor = 1.0/(0.5*RefDensity*RefArea*RefVel2);
+  AeroCoeffForceRef = 0.5 * RefDensity * RefArea * RefVel2;
+  const su2double factor = 1.0 / AeroCoeffForceRef;
 
   /*-- Variables initialization ---*/
   Total_CD = 0.0; Total_CL = 0.0; Total_CSF = 0.0; Total_CEff = 0.0;

SU2_CFD/src/solvers/CFEM_DG_NSSolver.cpp

@@ -209,33 +209,13 @@ void CFEM_DG_NSSolver::Friction_Forces(const CGeometry* geometry, const CConfig*
   /*--- Get the information of the angle of attack, reference area, etc. ---*/
   const su2double Alpha        = config->GetAoA()*PI_NUMBER/180.0;
   const su2double Beta         = config->GetAoS()*PI_NUMBER/180.0;
-  const su2double RefArea      = config->GetRefArea();
   const su2double RefLength    = config->GetRefLength();
-  const su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin      = config->GetRefOriginMoment(0);
-  const bool grid_movement     = config->GetGrid_Movement();
-
-/*--- Evaluate reference values for non-dimensionalization.
-        For dynamic meshes, use the motion Mach number as a reference value
-        for computing the force coefficients. Otherwise, use the freestream
-        values, which is the standard convention. ---*/
-  const su2double RefTemp     = Temperature_Inf;
-  const su2double RefDensity  = Density_Inf;
-  const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
 
-  su2double RefVel2;
-  if (grid_movement) {
-    const su2double Mach2Vel = sqrt(Gamma*Gas_Constant*RefTemp);
-    const su2double Mach_Motion = config->GetMach_Motion();
-    RefVel2 = (Mach_Motion*Mach2Vel)*(Mach_Motion*Mach2Vel);
-  }
-  else {
-    RefVel2 = 0.0;
-    for(unsigned short iDim=0; iDim<nDim; ++iDim)
-      RefVel2 += Velocity_Inf[iDim]*Velocity_Inf[iDim];
-  }
+  /*--- Evaluate reference values for non-dimensionalization. ---*/
+  const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
 
-  const su2double factor = 1.0/(0.5*RefDensity*RefArea*RefVel2);
+  const su2double factor = 1.0 / AeroCoeffForceRef;
 
   /*--- Variables initialization ---*/
   AllBound_CD_Visc = 0.0;  AllBound_CL_Visc  = 0.0; AllBound_CSF_Visc = 0.0;