Refactor of the NEMO solvers: shared helpers and removal of unreachable code

AUTHORS.md

@@ -50,6 +50,7 @@ Alessandro Gastaldi
 Aman uz zaman Baig
 Amit Sachdeva
 Ana Lourenco
+Andrea Rausa
 Andrew Burkett
 Andrew Wendorff
 Angelo Passariello

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -411,6 +411,22 @@ class CFVMFlowSolverBase : public CSolver {
   void LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter, bool update_geo,
                         su2double* RestartSolution = nullptr, unsigned short nVar_Restart = 0);
 
+  /*!
+   * \brief Compute the reference values and the non-dimensional free-stream values that are common
+   *        to the compressible solvers, and store them in the config.
+   * \note The reference pressure, density and temperature follow from the option REF_DIMENSIONALIZATION.
+   *       All other reference values are derived from these three, with a reference length of one.
+   * \param[in,out] config - Definition of the particular problem.
+   * \param[in] Pressure_FreeStream - Dimensional free-stream pressure.
+   * \param[in] Density_FreeStream - Dimensional free-stream density.
+   * \param[in] Temperature_FreeStream - Dimensional free-stream temperature.
+   * \param[in] Viscosity_FreeStream - Dimensional free-stream laminar viscosity.
+   * \param[in] ModVel_FreeStream - Dimensional magnitude of the free-stream velocity.
+   */
+  void SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream, su2double Density_FreeStream,
+                                   su2double Temperature_FreeStream, su2double Viscosity_FreeStream,
+                                   su2double ModVel_FreeStream) const;
+
   /*!
    * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -850,6 +850,112 @@ void CFVMFlowSolverBase<V, R>::UpdateCustomBoundaryConditions(
   }
 }
 
+template <class V, ENUM_REGIME R>
+void CFVMFlowSolverBase<V, R>::SetCompressibleNonDimValues(CConfig *config, su2double Pressure_FreeStream,
+                                                           su2double Density_FreeStream,
+                                                           su2double Temperature_FreeStream,
+                                                           su2double Viscosity_FreeStream,
+                                                           su2double ModVel_FreeStream) const {
+
+  const su2double Mach = config->GetMach();
+
+  /*--- Reference pressure, density and temperature. ---*/
+
+  su2double Pressure_Ref = 0.0, Density_Ref = 0.0, Temperature_Ref = 0.0;
+
+  if (config->GetRef_NonDim() == DIMENSIONAL) {
+    Pressure_Ref      = 1.0;
+    Density_Ref       = 1.0;
+    Temperature_Ref   = 1.0;
+  }
+  else if (config->GetRef_NonDim() == FREESTREAM_PRESS_EQ_ONE) {
+    Pressure_Ref      = Pressure_FreeStream;     // Pressure_FreeStream = 1.0
+    Density_Ref       = Density_FreeStream;      // Density_FreeStream = 1.0
+    Temperature_Ref   = Temperature_FreeStream;  // Temperature_FreeStream = 1.0
+  }
+  else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_MACH) {
+    Pressure_Ref      = Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/Gamma
+    Density_Ref       = Density_FreeStream;        // Density_FreeStream = 1.0
+    Temperature_Ref   = Temperature_FreeStream;    // Temp_FreeStream = 1.0
+  }
+  else if (config->GetRef_NonDim() == FREESTREAM_VEL_EQ_ONE) {
+    Pressure_Ref      = Mach*Mach*Gamma*Pressure_FreeStream; // Pressure_FreeStream = 1.0/(Gamma*(M_inf)^2)
+    Density_Ref       = Density_FreeStream;        // Density_FreeStream = 1.0
+    Temperature_Ref   = Temperature_FreeStream;    // Temp_FreeStream = 1.0
+  }
+  config->SetPressure_Ref(Pressure_Ref);
+  config->SetDensity_Ref(Density_Ref);
+  config->SetTemperature_Ref(Temperature_Ref);
+
+  /*--- Reference values that are derived from the three above. ---*/
+
+  const su2double Length_Ref       = 1.0;
+  const su2double Velocity_Ref     = sqrt(config->GetPressure_Ref()/config->GetDensity_Ref());
+  const su2double Time_Ref         = Length_Ref/Velocity_Ref;
+  const su2double Omega_Ref        = Velocity_Ref/Length_Ref;
+  const su2double Force_Ref        = config->GetDensity_Ref()*Velocity_Ref*Velocity_Ref*Length_Ref*Length_Ref;
+  const su2double Heat_Flux_Ref    = Density_Ref*Velocity_Ref*Velocity_Ref*Velocity_Ref;
+  const su2double Gas_Constant_Ref = Velocity_Ref*Velocity_Ref/config->GetTemperature_Ref();
+  const su2double Viscosity_Ref    = config->GetDensity_Ref()*Velocity_Ref*Length_Ref;
+  const su2double Conductivity_Ref = Viscosity_Ref*Gas_Constant_Ref;
+  const su2double Froude           = ModVel_FreeStream/sqrt(STANDARD_GRAVITY*Length_Ref);
+
+  config->SetLength_Ref(Length_Ref);
+  config->SetVelocity_Ref(Velocity_Ref);
+  config->SetTime_Ref(Time_Ref);
+  config->SetOmega_Ref(Omega_Ref);
+  config->SetForce_Ref(Force_Ref);
+  config->SetHeat_Flux_Ref(Heat_Flux_Ref);
+  config->SetGas_Constant_Ref(Gas_Constant_Ref);
+  config->SetViscosity_Ref(Viscosity_Ref);
+  config->SetConductivity_Ref(Conductivity_Ref);
+  config->SetFroude(Froude);
+
+  /*--- Divide by reference values, to compute the non-dimensional free-stream values. ---*/
+
+  const su2double Pressure_FreeStreamND = Pressure_FreeStream/config->GetPressure_Ref();
+  config->SetPressure_FreeStreamND(Pressure_FreeStreamND);
+
+  const su2double Density_FreeStreamND = Density_FreeStream/config->GetDensity_Ref();
+  config->SetDensity_FreeStreamND(Density_FreeStreamND);
+
+  su2double Velocity_FreeStreamND[MAXNDIM] = {0.0};
+  for (auto iDim = 0u; iDim < nDim; iDim++) {
+    Velocity_FreeStreamND[iDim] = config->GetVelocity_FreeStream()[iDim]/Velocity_Ref;
+    config->SetVelocity_FreeStreamND(Velocity_FreeStreamND[iDim], iDim);
+  }
+
+  const su2double Temperature_FreeStreamND = Temperature_FreeStream/config->GetTemperature_Ref();
+  config->SetTemperature_FreeStreamND(Temperature_FreeStreamND);
+
+  const su2double Gas_ConstantND = config->GetGas_Constant()/Gas_Constant_Ref;
+  config->SetGas_ConstantND(Gas_ConstantND);
+
+  const su2double ModVel_FreeStreamND = GeometryToolbox::Norm(nDim, Velocity_FreeStreamND);
+  config->SetModVel_FreeStreamND(ModVel_FreeStreamND);
+
+  const su2double Viscosity_FreeStreamND = Viscosity_FreeStream / Viscosity_Ref;
+  config->SetViscosity_FreeStreamND(Viscosity_FreeStreamND);
+
+  /*--- Free-stream turbulence quantities, from the turbulence intensity and the viscosity ratio. ---*/
+
+  const su2double TurbIntensity = config->GetTurbulenceIntensity_FreeStream();
+
+  const su2double Tke_FreeStream = 3.0/2.0*(ModVel_FreeStream*ModVel_FreeStream*TurbIntensity*TurbIntensity);
+  config->SetTke_FreeStream(Tke_FreeStream);
+
+  const su2double Tke_FreeStreamND = 3.0/2.0*(ModVel_FreeStreamND*ModVel_FreeStreamND*TurbIntensity*TurbIntensity);
+  config->SetTke_FreeStreamND(Tke_FreeStreamND);
+
+  const su2double Omega_FreeStream =
+      Density_FreeStream*Tke_FreeStream/max(Viscosity_FreeStream*config->GetTurb2LamViscRatio_FreeStream(), EPS);
+  config->SetOmega_FreeStream(Omega_FreeStream);
+
+  const su2double Omega_FreeStreamND =
+      Density_FreeStreamND*Tke_FreeStreamND/max(Viscosity_FreeStreamND*config->GetTurb2LamViscRatio_FreeStream(), EPS);
+  config->SetOmega_FreeStreamND(Omega_FreeStreamND);
+}
+
 template <class V, ENUM_REGIME R>
 void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver ***solver, CConfig *config, int iter,
                                                 bool update_geo, su2double* SolutionRestart,

SU2_CFD/include/solvers/CNEMOEulerSolver.hpp

@@ -60,6 +60,48 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
 
   CNEMOEulerVariable* node_infty = nullptr;
 
+  /*!
+   * \brief Pass the thermodynamic derivatives and vib.-el. quantities of two points to a numerics object.
+   * \note The ratio of specific heats (Gamma) is not passed.
+   * \param[in,out] numerics - Numerics object that receives the data.
+   * \param[in] nodes_i - Variable container of point i.
+   * \param[in] iPoint - Index of point i in its container.
+   * \param[in] nodes_j - Variable container of point j.
+   * \param[in] jPoint - Index of point j in its container.
+   */
+  static void SetNumericsThermoState(CNumerics* numerics, CNEMOEulerVariable* nodes_i, unsigned long iPoint,
+                                     CNEMOEulerVariable* nodes_j, unsigned long jPoint) {
+    numerics->SetdPdU  (nodes_i->GetdPdU(iPoint),   nodes_j->GetdPdU(jPoint));
+    numerics->SetdTdU  (nodes_i->GetdTdU(iPoint),   nodes_j->GetdTdU(jPoint));
+    numerics->SetdTvedU(nodes_i->GetdTvedU(iPoint), nodes_j->GetdTvedU(jPoint));
+    numerics->SetEve   (nodes_i->GetEve(iPoint),    nodes_j->GetEve(jPoint));
+    numerics->SetCvve  (nodes_i->GetCvve(iPoint),   nodes_j->GetCvve(jPoint));
+  }
+
+  /*!
+   * \brief Pass the transport properties of two points of the solver to a viscous numerics object.
+   * \param[in,out] numerics - Numerics object that receives the data.
+   * \param[in] iPoint - Index of point i.
+   * \param[in] jPoint - Index of point j.
+   */
+  void SetNumericsTransportState(CNumerics *numerics, unsigned long iPoint, unsigned long jPoint);
+
+  /*!
+   * \brief Compute the viscous flux between a boundary point and a prescribed boundary state,
+   *        and add it to the residual and to the Jacobian.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] visc_numerics - Description of the numerical method for viscous terms.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] iPoint - Index of the boundary point.
+   * \param[in] Point_Normal - Index of the closest interior neighbor of the boundary point.
+   * \param[in] Normal - Normal vector of the boundary face, the same that is passed to the convective numerics.
+   * \param[in] nodes_bound - Variable container that holds the boundary state.
+   * \param[in] iPoint_bound - Index of the boundary state in its container.
+   */
+  void AddViscousBoundaryResidual(CGeometry *geometry, CNumerics *visc_numerics, const CConfig *config,
+                                  unsigned long iPoint, unsigned long Point_Normal, const su2double *Normal,
+                                  CNEMOEulerVariable *nodes_bound, unsigned long iPoint_bound);
+
   /*!
    * \brief Set the maximum value of the eigenvalue.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CNEMONSSolver.hpp

@@ -72,6 +72,29 @@ class CNEMONSSolver final : public CNEMOEulerSolver {
    */
   void SetTau_Wall_WF(CGeometry *geometry, CSolver** solver_container, const CConfig* config);
 
+  /*!
+   * \brief Impose the velocity of a wall point as a strong (Dirichlet) boundary condition.
+   *        The velocity is fixed and the momentum residual and truncation error are set to zero.
+   * \param[in] iPoint - Index of the wall point.
+   * \param[in] velocity - Velocity to impose (zero for no-slip, slip velocity otherwise).
+   */
+  void SetWallVelocity_Strong(unsigned long iPoint, const su2double *velocity);
+
+  /*!
+   * \brief Compute the residual of the two energy equations that weakly imposes the wall temperatures.
+   * \param[in] iPoint - Index of the wall point.
+   * \param[in] jPoint - Index of the closest normal neighbor of the wall point.
+   * \param[in] Twall - Target translational-rotational temperature at the wall.
+   * \param[in] Tvewall - Target vibrational-electronic temperature at the wall.
+   * \param[in] pcontrol - "Proportional control" constant of the temperature difference term.
+   * \param[in] area - Area of the dual-grid face on the wall.
+   * \param[in] dist - Distance between the wall point and its normal neighbor.
+   * \param[out] residual - Residual, only the two energy entries are set.
+   */
+  void SetWallEnergyFlux_Weak(unsigned long iPoint, unsigned long jPoint, su2double Twall, su2double Tvewall,
+                              su2double pcontrol, su2double area, su2double dist,
+                              su2double *residual) const;
+
 public:
 
   /*!