Merge branch 'develop' into feature_small_restart

Author: bigfooted

Common/include/toolboxes/geometry_toolbox.hpp

@@ -79,6 +79,14 @@ inline T Norm(Int nDim, const T* a) {
   return sqrt(SquaredNorm(nDim, a));
 }
 
+/*! \brief dn = max(abs(n.(a-b)), 0.05 * ||a-b|| */
+template <class T, typename Int>
+inline T NormalDistance(Int nDim, const T* n, const T* a, const T* b) {
+  T d[3] = {0};
+  Distance(nDim, a, b, d);
+  return fmax(fabs(DotProduct(nDim, n, d)), 0.05 * Norm(nDim, d));
+}
+
 /*! \brief c = a x b */
 template <class T>
 inline void CrossProduct(const T* a, const T* b, T* c) {

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -345,18 +345,19 @@ struct Bsl {
 
     /*--- Limiting of \hat{S} based on "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model"
      * Note 1 option c in https://turbmodels.larc.nasa.gov/spalart.html ---*/
+    const su2double d_Sbar = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
     if (Sbar >= - c2 * var.Omega) {
       var.Shat = var.Omega + Sbar;
+      var.d_Shat = d_Sbar;
     } else {
       const su2double Num = var.Omega * (c2 * c2 * var.Omega + c3 * Sbar);
       const su2double Den = (c3 - 2 * c2) * var.Omega - Sbar;
       var.Shat = var.Omega + Num / Den;
+      var.d_Shat = d_Sbar * (c3 * var.Omega + Num / Den) / Den;
     }
     if (var.Shat <= 1e-10) {
       var.Shat = 1e-10;
       var.d_Shat = 0.0;
-    } else {
-      var.d_Shat = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
     }
   }
 };

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2357,12 +2357,11 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
   const su2double minYPlus = config->GetwallModel_MinYPlus();
 
-  string Marker_Tag, Monitoring_Tag;
-
   const su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0;
   const su2double Beta = config->GetAoS() * PI_NUMBER / 180.0;
   const su2double RefLength = config->GetRefLength();
   const su2double RefHeatFlux = config->GetHeat_Flux_Ref();
+  const su2double RefTemperature = config->GetTemperature_Ref();
   const su2double Gas_Constant = config->GetGas_ConstantND();
   auto Origin = config->GetRefOriginMoment(0);
 
@@ -2393,15 +2392,17 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
   for (iMarker = 0; iMarker < nMarker; iMarker++) {
 
-    Marker_Tag = config->GetMarker_All_TagBound(iMarker);
     if (!config->GetViscous_Wall(iMarker)) continue;
+    const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker);
+
+    const bool py_custom = config->GetMarker_All_PyCustom(iMarker);
 
     /*--- Obtain the origin for the moment computation for a particular marker ---*/
 
     const auto Monitoring = config->GetMarker_All_Monitoring(iMarker);
     if (Monitoring == YES) {
       for (iMarker_Monitoring = 0; iMarker_Monitoring < config->GetnMarker_Monitoring(); iMarker_Monitoring++) {
-        Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
+        const auto Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
         if (Marker_Tag == Monitoring_Tag) Origin = config->GetRefOriginMoment(iMarker_Monitoring);
       }
     }
@@ -2506,26 +2507,47 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
       /*--- Compute total and maximum heat flux on the wall ---*/
 
-      su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
-
-      if (!nemo){
-
+      if (!nemo) {
         if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
-
           Cp = (Gamma / Gamma_Minus_One) * Gas_Constant;
           thermal_conductivity = Cp * Viscosity / Prandtl_Lam;
         }
         if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE) {
-          if (!energy) dTdn = 0.0;
           thermal_conductivity = nodes->GetThermalConductivity(iPoint);
         }
-        HeatFlux[iMarker][iVertex] = -thermal_conductivity * dTdn * RefHeatFlux;
 
+        if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::HEAT_FLUX) {
+          if (py_custom) {
+            HeatFlux[iMarker][iVertex] = -geometry->GetCustomBoundaryHeatFlux(iMarker, iVertex);
+          } else {
+            HeatFlux[iMarker][iVertex] = -config->GetWall_HeatFlux(Marker_Tag);
+            if (config->GetIntegrated_HeatFlux()) {
+              HeatFlux[iMarker][iVertex] /= geometry->GetSurfaceArea(config, iMarker);
+            }
+          }
+        } else if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::ISOTHERMAL) {
+          su2double Twall = 0.0;
+          if (py_custom) {
+            Twall = geometry->GetCustomBoundaryTemperature(iMarker, iVertex) / RefTemperature;
+          } else {
+            Twall = config->GetIsothermal_Temperature(Marker_Tag) / RefTemperature;
+          }
+          iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
+          Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
+          const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord, Coord_Normal);
+          const su2double There = nodes->GetTemperature(iPointNormal);
+          HeatFlux[iMarker][iVertex] = thermal_conductivity * (There - Twall) / dist_ij * RefHeatFlux;
+        } else {
+          su2double dTdn = GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
+          if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE && !energy) dTdn = 0.0;
+          HeatFlux[iMarker][iVertex] = thermal_conductivity * dTdn * RefHeatFlux;
+        }
       } else {
 
         const auto& thermal_conductivity_tr = nodes->GetThermalConductivity(iPoint);
         const auto& thermal_conductivity_ve = nodes->GetThermalConductivity_ve(iPoint);
 
+        const su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
         const su2double dTvedn = -GeometryToolbox::DotProduct(nDim, Grad_Temp_ve, UnitNormal);
 
         /*--- Surface energy balance: trans-rot heat flux, vib-el heat flux ---*/
@@ -2653,8 +2675,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       /*--- Compute the coefficients per surface ---*/
 
       for (iMarker_Monitoring = 0; iMarker_Monitoring < config->GetnMarker_Monitoring(); iMarker_Monitoring++) {
-        Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
-        Marker_Tag = config->GetMarker_All_TagBound(iMarker);
+        const auto Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
         if (Marker_Tag == Monitoring_Tag) {
           SurfaceViscCoeff.CL[iMarker_Monitoring] += ViscCoeff.CL[iMarker];
           SurfaceViscCoeff.CD[iMarker_Monitoring] += ViscCoeff.CD[iMarker];

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -341,7 +341,6 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
   /*--- Variables for streamwise periodicity ---*/
   const bool streamwise_periodic = (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE);
   const bool streamwise_periodic_temperature = config->GetStreamwise_Periodic_Temperature();
-  su2double Cp, thermal_conductivity, dot_product, scalar_factor;
 
   /*--- Identify the boundary by string name ---*/
 
@@ -434,15 +433,17 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
       /*--- With streamwise periodic flow and heatflux walls an additional term is introduced in the boundary formulation ---*/
       if (streamwise_periodic && streamwise_periodic_temperature) {
 
-        Cp = nodes->GetSpecificHeatCp(iPoint);
-        thermal_conductivity = nodes->GetThermalConductivity(iPoint);
+        const su2double Cp = nodes->GetSpecificHeatCp(iPoint);
+        const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
 
         /*--- Scalar factor of the residual contribution ---*/
         const su2double norm2_translation = GeometryToolbox::SquaredNorm(nDim, config->GetPeriodic_Translation(0));
-        scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow*thermal_conductivity / (SPvals.Streamwise_Periodic_MassFlow * Cp * norm2_translation);
+        const su2double scalar_factor =
+            SPvals.Streamwise_Periodic_IntegratedHeatFlow*thermal_conductivity /
+            (SPvals.Streamwise_Periodic_MassFlow * Cp * norm2_translation);
 
         /*--- Dot product ---*/
-        dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);
+        const su2double dot_product = GeometryToolbox::DotProduct(nDim, config->GetPeriodic_Translation(0), Normal);
 
         LinSysRes(iPoint, nDim+1) += scalar_factor*dot_product;
       } // if streamwise_periodic
@@ -472,18 +473,17 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
 
       const auto Coord_i = geometry->nodes->GetCoord(iPoint);
       const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
-      su2double Edge_Vector[MAXNDIM];
-      GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector);
-      su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector);
-      su2double dist_ij = sqrt(dist_ij_2);
+      su2double UnitNormal[MAXNDIM] = {0.0};
+      for (auto iDim = 0u; iDim < nDim; ++iDim) UnitNormal[iDim] = Normal[iDim] / Area;
+      const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord_i, Coord_j);
 
       /*--- Compute the normal gradient in temperature using Twall ---*/
 
-      su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;
+      const su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij;
 
       /*--- Get thermal conductivity ---*/
 
-      su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
+      const su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint);
 
       /*--- Apply a weak boundary condition for the energy equation.
       Compute the residual due to the prescribed heat flux. ---*/
@@ -493,10 +493,7 @@ void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *con
       /*--- Jacobian contribution for temperature equation. ---*/
 
       if (implicit) {
-        su2double proj_vector_ij = 0.0;
-        if (dist_ij_2 > 0.0)
-          proj_vector_ij = GeometryToolbox::DotProduct(nDim, Edge_Vector, Normal) / dist_ij_2;
-        Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity*proj_vector_ij);
+        Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity * Area / dist_ij);
       }
       break;
     } // switch

SU2_CFD/src/solvers/CNSSolver.cpp

@@ -672,22 +672,19 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver
 
     const auto Coord_i = geometry->nodes->GetCoord(iPoint);
     const auto Coord_j = geometry->nodes->GetCoord(Point_Normal);
-
-    su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j);
+    const su2double dist_ij = GeometryToolbox::NormalDistance(nDim, UnitNormal, Coord_i, Coord_j);
 
     /*--- Store the corrected velocity at the wall which will
      be zero (v = 0), unless there is grid motion (v = u_wall)---*/
 
     if (dynamic_grid) {
       nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint));
-    }
-    else {
+    } else {
       su2double zero[MAXNDIM] = {0.0};
       nodes->SetVelocity_Old(iPoint, zero);
     }
 
-    for (auto iDim = 0u; iDim < nDim; iDim++)
-      LinSysRes(iPoint, iDim+1) = 0.0;
+    for (auto iDim = 0u; iDim < nDim; iDim++) LinSysRes(iPoint, iDim+1) = 0.0;
     nodes->SetVel_ResTruncError_Zero(iPoint);
 
     /*--- Get transport coefficients ---*/
@@ -708,8 +705,7 @@ void CNSSolver::BC_Isothermal_Wall_Generic(CGeometry *geometry, CSolver **solver
     if (cht_mode) {
       Twall = GetCHTWallTemperature(config, val_marker, iVertex, dist_ij,
                                     thermal_conductivity, There, Temperature_Ref);
-    }
-    else if (config->GetMarker_All_PyCustom(val_marker)) {
+    } else if (config->GetMarker_All_PyCustom(val_marker)) {
       Twall = geometry->GetCustomBoundaryTemperature(val_marker, iVertex) / Temperature_Ref;
     }
 

SU2_CFD/src/solvers/CSolver.cpp

@@ -3163,7 +3163,8 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c
   const unsigned long nFields = Restart_Vars[1];
   const unsigned long nPointFile = Restart_Vars[2];
   const auto t0 = SU2_MPI::Wtime();
-  auto nRecurse = 0;
+  int nRecurse = 0;
+  const int maxNRecurse = 128;
 
   if (rank == MASTER_NODE) {
     cout << "\nThe number of points in the restart file (" << nPointFile << ") does not match "
@@ -3262,7 +3263,7 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c
   bool done = false;
 
   SU2_OMP_PARALLEL
-  while (!done) {
+  while (!done && nRecurse < maxNRecurse) {
     SU2_OMP_FOR_DYN(roundUpDiv(nPointDomain,2*omp_get_num_threads()))
     for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) {
       /*--- Do not change points that are already interpolated. ---*/
@@ -3273,7 +3274,8 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c
 
       for (const auto jPoint : geometry->nodes->GetPoints(iPoint)) {
         if (!isMapped[jPoint]) continue;
-        if (boundary_i != geometry->nodes->GetSolidBoundary(jPoint)) continue;
+        /*--- Take data from anywhere if we are looping too many times. ---*/
+        if (boundary_i != geometry->nodes->GetSolidBoundary(jPoint) && nRecurse < 8) continue;
 
         nDonor[iPoint]++;
 
@@ -3315,6 +3317,10 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c
 
   } // everything goes out of scope except "localVars"
 
+  if (nRecurse == maxNRecurse) {
+    SU2_MPI::Error("Limit number of recursions reached, the meshes may be too different.", CURRENT_FUNCTION);
+  }
+
   /*--- Move to Restart_Data in ascending order of global index, which is how a matching restart would have been read. ---*/
 
   Restart_Data.resize(nPointDomain*nFields);
@@ -3329,9 +3335,11 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c
       counter++;
     }
   }
+  int nRecurseMax = 0;
+  SU2_MPI::Reduce(&nRecurse, &nRecurseMax, 1, MPI_INT, MPI_MAX, MASTER_NODE, SU2_MPI::GetComm());
 
   if (rank == MASTER_NODE) {
-    cout << "Number of recursions: " << nRecurse << ".\n"
+    cout << "Number of recursions: " << nRecurseMax << ".\n"
             "Elapsed time: " << SU2_MPI::Wtime()-t0 << "s.\n" << endl;
   }
 }