From bef2a1e63e171cce88f780ba972fda88fe05c438 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 8 Dec 2024 14:08:44 -0800
Subject: [PATCH 01/21] measure the heat flux used in BCs instead of computing
 from the normal gradient

---
 .../include/solvers/CFVMFlowSolverBase.inl    | 46 ++++++++++++++++---
 SU2_CFD/src/solvers/CIncNSSolver.cpp          | 18 ++++----
 SU2_CFD/src/solvers/CNSSolver.cpp             | 16 ++++---
 3 files changed, 57 insertions(+), 23 deletions(-)

diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 4797293529f..10166b18978 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -2363,6 +2363,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   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);
 
@@ -2396,6 +2397,8 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     Marker_Tag = config->GetMarker_All_TagBound(iMarker);
     if (!config->GetViscous_Wall(iMarker)) continue;
 
+    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);
@@ -2506,26 +2509,55 @@ 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);
+          su2double Vec_ij[MAXNDIM] = {0.0};
+          GeometryToolbox::Distance(nDim, Coord, Coord_Normal, Vec_ij);
+
+          /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+          const su2double dist_ij = fmax(
+            fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
+            fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
+
+          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 ---*/
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 41cedcdb5f5..67b8ed863c4 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -472,18 +472,21 @@ 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];
+      su2double Edge_Vector[MAXNDIM] = {0.0};
       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);
+
+      /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+      const su2double dist_ij = fmax(
+        fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Edge_Vector, Normal)) / Area,
+        fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Edge_Vector), EPS));
 
       /*--- 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 +496,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
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
index 127ff3f5b08..3d6c6609e3b 100644
--- a/SU2_CFD/src/solvers/CNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -672,22 +672,25 @@ 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 Vec_ij[MAXNDIM] = {0.0};
+    GeometryToolbox::Distance(nDim, Coord_i, Coord_j, Vec_ij);
 
-    su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_i, Coord_j);
+    /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
+    const su2double dist_ij = fmax(
+      fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
+      fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
 
     /*--- 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 +711,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;
     }
 

From 43cb5f23d68b14f6f773d2130ef7a52d907ada40 Mon Sep 17 00:00:00 2001
From: emaberman <emaberman@gmail.com>
Date: Mon, 9 Dec 2024 17:30:56 +0200
Subject: [PATCH 02/21] fix bug with dshat for SA jacobian

---
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index c9464b1dfc3..0614e451714 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -347,16 +347,18 @@ struct Bsl {
      * Note 1 option c in https://turbmodels.larc.nasa.gov/spalart.html ---*/
     if (Sbar >= - c2 * var.Omega) {
       var.Shat = var.Omega + Sbar;
+      var.d_Shat = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
     } 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;
+      
+      const su2double d_SB =var.d_Shat = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
+      var.d_Shat = d_SB * (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;
     }
   }
 };

From fefd2d343eb98a6abec7da6464d332915256fe29 Mon Sep 17 00:00:00 2001
From: emaberman <emaberman@gmail.com>
Date: Mon, 9 Dec 2024 17:38:21 +0200
Subject: [PATCH 03/21] fix bug with dshat for SA jacobian

---
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 0614e451714..5a173f6aeea 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -353,7 +353,7 @@ struct Bsl {
       const su2double Den = (c3 - 2 * c2) * var.Omega - Sbar;
       var.Shat = var.Omega + Num / Den;
       
-      const su2double d_SB =var.d_Shat = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
+      const su2double d_SB = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
       var.d_Shat = d_SB * (c3 * var.Omega + Num / Den) / Den;
     }
     if (var.Shat <= 1e-10) {

From 147da5d4eda752f84512a7b260b2f2594d7042a7 Mon Sep 17 00:00:00 2001
From: Eitan Aberman <139676851+emaberman@users.noreply.github.com>
Date: Tue, 10 Dec 2024 08:02:44 +0200
Subject: [PATCH 04/21] Update
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
---
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 5a173f6aeea..8698e2822c7 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -345,6 +345,7 @@ 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 = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;

From b6ebcc5c725b6c72f13141626d115937e239b297 Mon Sep 17 00:00:00 2001
From: Eitan Aberman <139676851+emaberman@users.noreply.github.com>
Date: Tue, 10 Dec 2024 08:02:52 +0200
Subject: [PATCH 05/21] Update
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
---
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 8698e2822c7..f2e3d1c3809 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -348,7 +348,7 @@ struct Bsl {
     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 = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
+      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;

From 8214318f56150b3e9c7648e13aef1689fef4b628 Mon Sep 17 00:00:00 2001
From: Eitan Aberman <139676851+emaberman@users.noreply.github.com>
Date: Tue, 10 Dec 2024 08:02:58 +0200
Subject: [PATCH 06/21] Update
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp

Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com>
---
 SU2_CFD/include/numerics/turbulent/turb_sources.hpp | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index f2e3d1c3809..4e5f1cd51f3 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -353,9 +353,7 @@ struct Bsl {
       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;
-      
-      const su2double d_SB = (var.fv2 + nue * var.d_fv2) * var.inv_k2_d2;
-      var.d_Shat = d_SB * (c3 * var.Omega + Num / Den) / Den;
+      var.d_Shat = d_Sbar * (c3 * var.Omega + Num / Den) / Den;
     }
     if (var.Shat <= 1e-10) {
       var.Shat = 1e-10;

From 397c22819d41b7669d7fb234db5115917acd05ef Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 15 Dec 2024 11:48:21 -0800
Subject: [PATCH 07/21] recover old sign convention and use helper function

---
 Common/include/toolboxes/geometry_toolbox.hpp |  8 ++++++++
 .../include/solvers/CFVMFlowSolverBase.inl    | 20 ++++++-------------
 SU2_CFD/src/solvers/CIncNSSolver.cpp          | 10 +++-------
 SU2_CFD/src/solvers/CNSSolver.cpp             |  8 +-------
 4 files changed, 18 insertions(+), 28 deletions(-)

diff --git a/Common/include/toolboxes/geometry_toolbox.hpp b/Common/include/toolboxes/geometry_toolbox.hpp
index 773f1d553dc..ddc2fd4434a 100644
--- a/Common/include/toolboxes/geometry_toolbox.hpp
+++ b/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) {
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 10166b18978..40206768e4c 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -2520,9 +2520,9 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
 
         if (config->GetMarker_All_KindBC(iMarker) == BC_TYPE::HEAT_FLUX) {
           if (py_custom) {
-            HeatFlux[iMarker][iVertex] = geometry->GetCustomBoundaryHeatFlux(iMarker, iVertex);
+            HeatFlux[iMarker][iVertex] = -geometry->GetCustomBoundaryHeatFlux(iMarker, iVertex);
           } else {
-            HeatFlux[iMarker][iVertex] = config->GetWall_HeatFlux(Marker_Tag);
+            HeatFlux[iMarker][iVertex] = -config->GetWall_HeatFlux(Marker_Tag);
             if (config->GetIntegrated_HeatFlux()) {
               HeatFlux[iMarker][iVertex] /= geometry->GetSurfaceArea(config, iMarker);
             }
@@ -2536,21 +2536,13 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
           }
           iPointNormal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor();
           Coord_Normal = geometry->nodes->GetCoord(iPointNormal);
-          su2double Vec_ij[MAXNDIM] = {0.0};
-          GeometryToolbox::Distance(nDim, Coord, Coord_Normal, Vec_ij);
-
-          /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
-          const su2double dist_ij = fmax(
-            fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
-            fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
-
+          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;
+          HeatFlux[iMarker][iVertex] = thermal_conductivity * (There - Twall) / dist_ij * RefHeatFlux;
         } else {
-          su2double dTdn = -GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
+          su2double dTdn = GeometryToolbox::DotProduct(nDim, Grad_Temp, UnitNormal);
           if (FlowRegime == ENUM_REGIME::INCOMPRESSIBLE && !energy) dTdn = 0.0;
-          HeatFlux[iMarker][iVertex] = -thermal_conductivity * dTdn * RefHeatFlux;
+          HeatFlux[iMarker][iVertex] = thermal_conductivity * dTdn * RefHeatFlux;
         }
       } else {
 
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 67b8ed863c4..8db972191c3 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -472,13 +472,9 @@ 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] = {0.0};
-      GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector);
-
-      /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
-      const su2double dist_ij = fmax(
-        fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Edge_Vector, Normal)) / Area,
-        fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Edge_Vector), EPS));
+      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 ---*/
 
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
index 3d6c6609e3b..50f243c7f0c 100644
--- a/SU2_CFD/src/solvers/CNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -672,13 +672,7 @@ 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 Vec_ij[MAXNDIM] = {0.0};
-    GeometryToolbox::Distance(nDim, Coord_i, Coord_j, Vec_ij);
-
-    /*--- Prevent divisions by 0 by limiting the normal projection. ---*/
-    const su2double dist_ij = fmax(
-      fabs(GeometryToolbox::DotProduct(int(MAXNDIM), Vec_ij, UnitNormal)),
-      fmax(0.05 * GeometryToolbox::Norm(int(MAXNDIM), Vec_ij), EPS));
+    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)---*/

From 82ac79102671c041c339668316ff2177f1419ecd Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 15 Dec 2024 13:40:19 -0800
Subject: [PATCH 08/21] update tests

---
 TestCases/hybrid_regression.py      | 34 +++++++++++-----------
 TestCases/hybrid_regression_AD.py   |  4 +--
 TestCases/parallel_regression.py    | 44 ++++++++++++++---------------
 TestCases/parallel_regression_AD.py |  4 +--
 TestCases/serial_regression.py      | 42 +++++++++++++--------------
 TestCases/serial_regression_AD.py   |  4 +--
 6 files changed, 66 insertions(+), 66 deletions(-)

diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index fe2634fb2a5..a5fbd5961b6 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -103,7 +103,7 @@ def main():
     flatplate.cfg_dir   = "navierstokes/flatplate"
     flatplate.cfg_file  = "lam_flatplate.cfg"
     flatplate.test_iter = 100
-    flatplate.test_vals = [-7.679131, -2.206953, 0.001084, 0.036233, 2.361500, -2.325300, -1.984700, -1.984700]
+    flatplate.test_vals = [-7.679131, -2.206953, 0.001084, 0.036233, 2.361500, -2.325300, 0, 0]
     test_list.append(flatplate)
 
     # Laminar cylinder (steady)
@@ -111,7 +111,7 @@ def main():
     cylinder.cfg_dir   = "navierstokes/cylinder"
     cylinder.cfg_file  = "lam_cylinder.cfg"
     cylinder.test_iter = 25
-    cylinder.test_vals = [-8.266513, -2.783904, -0.019899, 1.615668, -0.010207]
+    cylinder.test_vals = [-8.266513, -2.783904, -0.019899, 1.615668, 0]
     test_list.append(cylinder)
 
     # Laminar cylinder (low Mach correction)
@@ -119,8 +119,8 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals         = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
-    cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0.002457]
+    cylinder_lowmach.test_vals         = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
+    cylinder_lowmach.test_vals_aarch64 = [-6.830996, -1.368850, -0.143956, 73.963354, 0]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)
@@ -128,7 +128,7 @@ def main():
     poiseuille.cfg_dir   = "navierstokes/poiseuille"
     poiseuille.cfg_file  = "lam_poiseuille.cfg"
     poiseuille.test_iter = 10
-    poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735818, -2.142500]
+    poiseuille.test_vals = [-5.046131, 0.652984, 0.008355, 13.735818, 0]
     test_list.append(poiseuille)
 
     # 2D Poiseuille flow (inlet profile file)
@@ -157,7 +157,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, -80602.000000]
+    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, 0]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -165,7 +165,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-0.510363, 4.872736, 0.815617, 0.060920, -73391.000000]
+    rae2822_sst.test_vals = [-0.510363, 4.872736, 0.815617, 0.060920, 0]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -189,7 +189,7 @@ def main():
     turb_oneram6.cfg_dir   = "rans/oneram6"
     turb_oneram6.cfg_file  = "turb_ONERAM6.cfg"
     turb_oneram6.test_iter = 10
-    turb_oneram6.test_vals = [-2.408523, -6.662833, 0.238333, 0.158910, -52718]
+    turb_oneram6.test_vals = [-2.408523, -6.662833, 0.238333, 0.158910, 0]
     test_list.append(turb_oneram6)
 
     # NACA0012 (SA, FUN3D finest grid results: CL=1.0983, CD=0.01242)
@@ -197,8 +197,8 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals         = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, -44.871000]
-    turb_naca0012_sa.test_vals_aarch64 = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, -44.871000]
+    turb_naca0012_sa.test_vals         = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0]
+    turb_naca0012_sa.test_vals_aarch64 = [-12.098325, -14.149988, 1.057665, 0.022971, 20.000000, -2.292707, 0.000000, -12.068169, 0]
     test_list.append(turb_naca0012_sa)
 
     # NACA0012 (SST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
@@ -206,7 +206,7 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals = [-12.105781, -15.277738, -6.210248, 1.049757, 0.019249, -2.807857, -38.976000]
+    turb_naca0012_sst.test_vals = [-12.105781, -15.277738, -6.210248, 1.049757, 0.019249, -2.807857, 0]
     test_list.append(turb_naca0012_sst)
 
     # NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
@@ -250,7 +250,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-12.070954, -7.407644, -8.698118, -4.008751, -3572.100000]
+    axi_rans_air_nozzle_restart.test_vals = [-12.070954, -7.407644, -8.698118, -4.008751, 0]
     test_list.append(axi_rans_air_nozzle_restart)
 
     #################################
@@ -381,8 +381,8 @@ def main():
     inc_poly_cylinder.cfg_dir   = "incomp_navierstokes/cylinder"
     inc_poly_cylinder.cfg_file  = "poly_cylinder.cfg"
     inc_poly_cylinder.test_iter = 20
-    inc_poly_cylinder.test_vals         = [-7.851512, -2.093420, 0.029974, 1.921595, -175.300000]
-    inc_poly_cylinder.test_vals_aarch64 = [-7.851510, -2.093419, 0.029974, 1.921595, -175.300000]
+    inc_poly_cylinder.test_vals         = [-7.827942, -2.061513, 0.029525, 1.953498, -174.780000]
+    inc_poly_cylinder.test_vals_aarch64 = [-7.827942, -2.061513, 0.029525, 1.953498, -174.780000]
     test_list.append(inc_poly_cylinder)
 
     # X-coarse laminar bend as a mixed element CGNS test
@@ -452,8 +452,8 @@ def main():
     square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     square_cylinder.cfg_file  = "turb_square.cfg"
     square_cylinder.test_iter = 3
-    square_cylinder.test_vals = [-2.560839, -1.173497, 0.061188, 1.399403, 2.220575, 1.399351, 2.218781, -0.584690]
-    square_cylinder.test_vals_aarch64 = [-2.557902, -1.173574, 0.058050, 1.399794, 2.220402, 1.399748, 2.218604, -0.453270]
+    square_cylinder.test_vals = [-2.560839, -1.173497, 0.061188, 1.399403, 2.220575, 1.399351, 2.218781, 0]
+    square_cylinder.test_vals_aarch64 = [-2.557902, -1.173574, 0.058050, 1.399794, 2.220402, 1.399748, 2.218604, 0]
     square_cylinder.unsteady  = True
     test_list.append(square_cylinder)
 
@@ -503,7 +503,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882681, -0.215041, 0.023758, -617.440000]
+    ddes_flatplate.test_vals = [-2.714785, -5.882681, -0.215041, 0.023758, 0]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py
index a0e29d54c84..3cc8adf1a68 100644
--- a/TestCases/hybrid_regression_AD.py
+++ b/TestCases/hybrid_regression_AD.py
@@ -110,8 +110,8 @@ def main():
     discadj_incomp_cylinder.cfg_dir   = "disc_adj_incomp_navierstokes/cylinder"
     discadj_incomp_cylinder.cfg_file  = "heated_cylinder.cfg"
     discadj_incomp_cylinder.test_iter = 20
-    discadj_incomp_cylinder.test_vals         = [20.000000, -2.705921, -2.837904, 0.000000]
-    discadj_incomp_cylinder.test_vals_aarch64 = [20.000000, -2.705918, -2.837766, 0.000000]
+    discadj_incomp_cylinder.test_vals         = [20.000000, -2.746353, -2.934792, 0.000000]
+    discadj_incomp_cylinder.test_vals_aarch64 = [20.000000, -2.746353, -2.934792, 0.000000]
     test_list.append(discadj_incomp_cylinder)
 
     ######################################
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 5462560119e..e5f3f7ae029 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -303,7 +303,7 @@ def main():
     flatplate.cfg_dir   = "navierstokes/flatplate"
     flatplate.cfg_file  = "lam_flatplate.cfg"
     flatplate.test_iter = 100
-    flatplate.test_vals = [-7.613122, -2.140941, 0.001084, 0.036230, 2.361500, -2.325300, -1.975600, -1.975600]
+    flatplate.test_vals = [-7.613122, -2.140941, 0.001084, 0.036230, 2.361500, -2.325300, 0.000000, 0.000000]
     test_list.append(flatplate)
 
     # Custom objective function
@@ -319,7 +319,7 @@ def main():
     cylinder.cfg_dir   = "navierstokes/cylinder"
     cylinder.cfg_file  = "lam_cylinder.cfg"
     cylinder.test_iter = 25
-    cylinder.test_vals = [-8.422091, -2.930561, -0.003396, 1.608418, -0.009920]
+    cylinder.test_vals = [-8.422091, -2.930561, -0.003396, 1.608418, 0.000000]
     test_list.append(cylinder)
 
     # Laminar cylinder (low Mach correction)
@@ -327,7 +327,7 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals = [-6.841604, -1.379532, -1.266739, 76.118218, 0.000274]
+    cylinder_lowmach.test_vals = [-6.841604, -1.379532, -1.266739, 76.118218, 0.000000]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)
@@ -335,7 +335,7 @@ def main():
     poiseuille.cfg_dir   = "navierstokes/poiseuille"
     poiseuille.cfg_file  = "lam_poiseuille.cfg"
     poiseuille.test_iter = 10
-    poiseuille.test_vals = [-5.050889, 0.648196, 0.000199, 13.639173, -2.114500]
+    poiseuille.test_vals = [-5.050889, 0.648196, 0.000199, 13.639173, 0.000000]
     poiseuille.tol       = 0.001
     test_list.append(poiseuille)
 
@@ -357,7 +357,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.004689, -5.265782, 0.809463, 0.062016, -80576.000000]
+    rae2822_sa.test_vals = [-2.004689, -5.265782, 0.809463, 0.062016, 0.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -365,7 +365,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-0.510369, 4.870643, 0.816647, 0.061833, -71779.000000]
+    rae2822_sst.test_vals = [-0.510369, 4.870643, 0.816647, 0.061833, 0.000000]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -413,7 +413,7 @@ def main():
     turb_oneram6.cfg_dir   = "rans/oneram6"
     turb_oneram6.cfg_file  = "turb_ONERAM6.cfg"
     turb_oneram6.test_iter = 10
-    turb_oneram6.test_vals = [-2.408532, -6.662836, 0.238334, 0.158910, -52718]
+    turb_oneram6.test_vals = [-2.408532, -6.662836, 0.238334, 0.158910, 0]
     turb_oneram6.timeout   = 3200
     test_list.append(turb_oneram6)
 
@@ -422,7 +422,7 @@ def main():
     turb_oneram6_vc.cfg_dir = "rans/oneram6"
     turb_oneram6_vc.cfg_file = "turb_ONERAM6_vc.cfg"
     turb_oneram6_vc.test_iter = 15
-    turb_oneram6_vc.test_vals = [-2.281896, -6.614616, 0.233785, 0.142861, -47942]
+    turb_oneram6_vc.test_vals = [-2.281896, -6.614616, 0.233785, 0.142861, 0]
     turb_oneram6_vc.timeout = 3200
     test_list.append(turb_oneram6_vc)
 
@@ -441,7 +441,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.094721, -14.685268, 1.057665, 0.022971, 20.000000, -1.692925, 20.000000, -4.037672, -44.871000]
+    turb_naca0012_sa.test_vals = [-12.094721, -14.685268, 1.057665, 0.022971, 20.000000, -1.692925, 20.000000, -4.037672, 0]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
 
@@ -450,8 +450,8 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, -38.976000]
-    turb_naca0012_sst.test_vals_aarch64 = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, -38.976000]
+    turb_naca0012_sst.test_vals = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, 0]
+    turb_naca0012_sst.test_vals_aarch64 = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, 0]
     turb_naca0012_sst.timeout   = 3200
     test_list.append(turb_naca0012_sst)
 
@@ -530,7 +530,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-12.071395, -7.467871, -8.649076, -3.995810, -3572.100000]
+    axi_rans_air_nozzle_restart.test_vals = [-12.071395, -7.467871, -8.649076, -3.995810, 0]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)
 
@@ -602,7 +602,7 @@ def main():
     inc_poly_cylinder.cfg_dir   = "incomp_navierstokes/cylinder"
     inc_poly_cylinder.cfg_file  = "poly_cylinder.cfg"
     inc_poly_cylinder.test_iter = 20
-    inc_poly_cylinder.test_vals = [-7.791831, -2.062292, 0.013040, 1.913997, -171.120000]
+    inc_poly_cylinder.test_vals = [-7.786564, -2.036735, 0.012980, 1.944887, -170.930000]
     test_list.append(inc_poly_cylinder)
 
     # X-coarse laminar bend as a mixed element CGNS test
@@ -732,8 +732,8 @@ def main():
     turbmod_sa_neg_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_neg_rae2822.cfg_file  = "turb_SA_NEG_RAE2822.cfg"
     turbmod_sa_neg_rae2822.test_iter = 10
-    turbmod_sa_neg_rae2822.test_vals         = [-1.466238, 3.169232, 2.756518, 4.722767, 1.120253, 0.378675, -83444.000000]
-    turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.359612, 1.493629, 1.218367, -1.441703, 1.248499, 0.457987, -86467.000000]
+    turbmod_sa_neg_rae2822.test_vals         = [-1.466238, 3.169232, 2.756518, 4.722767, 1.120253, 0.378675, 0]
+    turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.359612, 1.493629, 1.218367, -1.441703, 1.248499, 0.457987, 0]
     test_list.append(turbmod_sa_neg_rae2822)
 
     # SA Compressibility Correction
@@ -974,7 +974,7 @@ def main():
     square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     square_cylinder.cfg_file  = "turb_square.cfg"
     square_cylinder.test_iter = 3
-    square_cylinder.test_vals = [-1.173495, 0.061186, 1.399404, 2.220578, 1.399352, 2.218783, -0.584750]
+    square_cylinder.test_vals = [-1.173495, 0.061186, 1.399404, 2.220578, 1.399352, 2.218783, 0]
     square_cylinder.unsteady  = True
     test_list.append(square_cylinder)
 
@@ -1001,7 +1001,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882679, -0.215041, 0.023758, -617.450000]
+    ddes_flatplate.test_vals = [-2.714785, -5.882679, -0.215041, 0.023758, 0]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
@@ -1282,7 +1282,7 @@ def main():
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
     p1rad.test_iter = 100
-    p1rad.test_vals = [-7.743666, -7.921411, -2.111848, 0.098302, -45.023000]
+    p1rad.test_vals = [-7.743666, -7.921411, -2.111848, 0.098302, -47.897000]
     test_list.append(p1rad)
 
 
@@ -1359,8 +1359,8 @@ def main():
     pywrapper_turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     pywrapper_turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     pywrapper_turb_naca0012_sst.test_iter = 10
-    pywrapper_turb_naca0012_sst.test_vals = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, -38.976000]
-    pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, -38.976000]
+    pywrapper_turb_naca0012_sst.test_vals = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, 0]
+    pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.107692, -15.277743, -6.210238, 1.049757, 0.019249, -2.357984, 0]
     pywrapper_turb_naca0012_sst.command   = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f")
     pywrapper_turb_naca0012_sst.timeout   = 3200
     test_list.append(pywrapper_turb_naca0012_sst)
@@ -1370,7 +1370,7 @@ def main():
     pywrapper_square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     pywrapper_square_cylinder.cfg_file  = "turb_square.cfg"
     pywrapper_square_cylinder.test_iter = 10
-    pywrapper_square_cylinder.test_vals = [-1.176405, -0.354027, 1.407859, 2.360784, 1.404715, 2.302615, -0.394750]
+    pywrapper_square_cylinder.test_vals = [-1.176405, -0.354027, 1.407859, 2.360784, 1.404715, 2.302615, 0]
     pywrapper_square_cylinder.command   = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f")
     pywrapper_square_cylinder.unsteady  = True
     test_list.append(pywrapper_square_cylinder)
@@ -1422,7 +1422,7 @@ def main():
     pywrapper_unsteadyCHT.cfg_dir       = "py_wrapper/flatPlate_unsteady_CHT"
     pywrapper_unsteadyCHT.cfg_file      = "unsteady_CHT_FlatPlate_Conf.cfg"
     pywrapper_unsteadyCHT.test_iter     = 5
-    pywrapper_unsteadyCHT.test_vals     = [-1.614167, 2.240663, -0.001316, 0.177491]
+    pywrapper_unsteadyCHT.test_vals     = [-1.614167, 2.240663, -0.001392, 0.177499]
     pywrapper_unsteadyCHT.command       = TestCase.Command("mpirun -np 2", "python", "launch_unsteady_CHT_FlatPlate.py --parallel -f")
     pywrapper_unsteadyCHT.unsteady      = True
     test_list.append(pywrapper_unsteadyCHT)
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 7d0f0d10049..dbc1ee9bd9b 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -116,7 +116,7 @@ def main():
     discadj_incomp_cylinder.cfg_dir   = "disc_adj_incomp_navierstokes/cylinder"
     discadj_incomp_cylinder.cfg_file  = "heated_cylinder.cfg"
     discadj_incomp_cylinder.test_iter = 20
-    discadj_incomp_cylinder.test_vals = [20.000000, -2.195581, -2.162081, 0.000000]
+    discadj_incomp_cylinder.test_vals = [20.000000, -2.082673, -2.013587, 0.000000]
     test_list.append(discadj_incomp_cylinder)
 
     ######################################
@@ -256,7 +256,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals         = [-2.117791, 0.784475, 0.000000, -0.574700]
+    discadj_heat.test_vals         = [-1.840132, 0.751835, 0.000000, 0.006760]
     discadj_heat.test_vals_aarch64 = [-2.226539, 0.605868, 0.000000, -6.256400]
     test_list.append(discadj_heat)
 
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 10312c92cc3..0ea5038cfc5 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -166,7 +166,7 @@ def main():
     flatplate.cfg_dir   = "navierstokes/flatplate"
     flatplate.cfg_file  = "lam_flatplate.cfg"
     flatplate.test_iter = 20
-    flatplate.test_vals = [-5.097707, 0.381809, 0.001324, 0.027932, 2.361600, -2.333600, -2.845200, -2.845200]
+    flatplate.test_vals = [-5.097707, 0.381809, 0.001324, 0.027932, 2.361600, -2.333600, 0, 0]
     test_list.append(flatplate)
 
     # Laminar cylinder (steady)
@@ -174,7 +174,7 @@ def main():
     cylinder.cfg_dir   = "navierstokes/cylinder"
     cylinder.cfg_file  = "lam_cylinder.cfg"
     cylinder.test_iter = 25
-    cylinder.test_vals = [-8.363995, -2.882536, -0.017780, 1.607979, -0.010080]
+    cylinder.test_vals = [-8.363995, -2.882536, -0.017780, 1.607979, 0]
     test_list.append(cylinder)
 
     # Laminar cylinder (low Mach correction)
@@ -182,7 +182,7 @@ def main():
     cylinder_lowmach.cfg_dir   = "navierstokes/cylinder"
     cylinder_lowmach.cfg_file  = "cylinder_lowmach.cfg"
     cylinder_lowmach.test_iter = 25
-    cylinder_lowmach.test_vals = [-6.830989, -1.368842, -0.143838, 73.962440, 0.002454]
+    cylinder_lowmach.test_vals = [-6.830989, -1.368842, -0.143838, 73.962440, 0]
     test_list.append(cylinder_lowmach)
 
     # 2D Poiseuille flow (body force driven with periodic inlet / outlet)
@@ -190,7 +190,7 @@ def main():
     poiseuille.cfg_dir   = "navierstokes/poiseuille"
     poiseuille.cfg_file  = "lam_poiseuille.cfg"
     poiseuille.test_iter = 10
-    poiseuille.test_vals = [-5.050753, 0.648333, 0.012273, 13.643141, -2.114700]
+    poiseuille.test_vals = [-5.050753, 0.648333, 0.012273, 13.643141, 0]
     test_list.append(poiseuille)
 
     # 2D Poiseuille flow (inlet profile file)
@@ -218,7 +218,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, -80602.000000]
+    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, 0]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -226,7 +226,7 @@ def main():
     rae2822_sst.cfg_dir   = "rans/rae2822"
     rae2822_sst.cfg_file  = "turb_SST_RAE2822.cfg"
     rae2822_sst.test_iter = 20
-    rae2822_sst.test_vals = [-0.510376, 4.873617, 0.816197, 0.060504, -71886.000000]
+    rae2822_sst.test_vals = [-0.510376, 4.873617, 0.816197, 0.060504, 0]
     test_list.append(rae2822_sst)
 
     # RAE2822 SST_SUST
@@ -266,7 +266,7 @@ def main():
     turb_oneram6.cfg_dir   = "rans/oneram6"
     turb_oneram6.cfg_file  = "turb_ONERAM6.cfg"
     turb_oneram6.test_iter = 10
-    turb_oneram6.test_vals = [-2.408533, -6.662837, 0.238334, 0.158910, -52718]
+    turb_oneram6.test_vals = [-2.408533, -6.662837, 0.238334, 0.158910, 0]
     turb_oneram6.timeout   = 3200
     test_list.append(turb_oneram6)
 
@@ -275,7 +275,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.091696, -14.685322, 1.057665, 0.022971, 20.000000, -2.686306, 20.000000, -4.459916, -44.871000]
+    turb_naca0012_sa.test_vals = [-12.091696, -14.685322, 1.057665, 0.022971, 20.000000, -2.686306, 20.000000, -4.459916, 0]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
 
@@ -284,8 +284,8 @@ def main():
     turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     turb_naca0012_sst.test_iter = 10
-    turb_naca0012_sst.test_vals         = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, -38.97600]
-    turb_naca0012_sst.test_vals_aarch64 = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, -38.97600]
+    turb_naca0012_sst.test_vals         = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, 0]
+    turb_naca0012_sst.test_vals_aarch64 = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, 0]
     turb_naca0012_sst.timeout   = 3200
     test_list.append(turb_naca0012_sst)
 
@@ -294,7 +294,7 @@ def main():
     turb_naca0012_sst_2003m.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_2003m.cfg_file  = "turb_NACA0012_sst_2003m.cfg"
     turb_naca0012_sst_2003m.test_iter = 10
-    turb_naca0012_sst_2003m.test_vals = [-7.688139, -10.046053, -3.410061, 1.048970, 0.019798, -2.208236, -45.678000]
+    turb_naca0012_sst_2003m.test_vals = [-7.688139, -10.046053, -3.410061, 1.048970, 0.019798, -2.208236, 0]
     turb_naca0012_sst_2003m.timeout   = 3200
     test_list.append(turb_naca0012_sst_2003m)
 
@@ -335,8 +335,8 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-12.071702, -7.474599, -8.646498, -3.988633, -3572.100000]
-    axi_rans_air_nozzle_restart.test_vals_aarch64 = [-12.071702, -7.474599, -8.646498, -3.988633, -3572.100000]
+    axi_rans_air_nozzle_restart.test_vals = [-12.071702, -7.474599, -8.646498, -3.988633, 0]
+    axi_rans_air_nozzle_restart.test_vals_aarch64 = [-12.071702, -7.474599, -8.646498, -3.988633, 0]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)
 
@@ -414,7 +414,7 @@ def main():
     inc_poly_cylinder.cfg_dir   = "incomp_navierstokes/cylinder"
     inc_poly_cylinder.cfg_file  = "poly_cylinder.cfg"
     inc_poly_cylinder.test_iter = 20
-    inc_poly_cylinder.test_vals = [-8.106741, -2.160042, 0.019225, 1.902421, -172.750000]
+    inc_poly_cylinder.test_vals = [-8.083556, -2.134369, 0.018999, 1.932938, -173.730000]
     test_list.append(inc_poly_cylinder)
 
     # X-coarse laminar bend as a mixed element CGNS test
@@ -457,7 +457,7 @@ def main():
     inc_turb_wallfunction_flatplate_sst.cfg_dir   = "wallfunctions/flatplate/incompressible_SST"
     inc_turb_wallfunction_flatplate_sst.cfg_file  = "turb_SST_flatplate.cfg"
     inc_turb_wallfunction_flatplate_sst.test_iter = 10
-    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.507362, -5.693894, -6.434063, -4.223774, -7.008049, -1.954102, 10.000000, -3.047554, 0.001081, 0.003644, 0.618140]
+    inc_turb_wallfunction_flatplate_sst.test_vals = [-6.507362, -5.693894, -6.434063, -4.223774, -7.008049, -1.954102, 10.000000, -3.047554, 0.001081, 0.003644, 0]
     test_list.append(inc_turb_wallfunction_flatplate_sst)
 
     # FLAT PLATE, WALL FUNCTIONS, INCOMPRESSIBLE SA
@@ -769,7 +769,7 @@ def main():
     square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     square_cylinder.cfg_file  = "turb_square.cfg"
     square_cylinder.test_iter = 3
-    square_cylinder.test_vals = [-2.560840, -1.173495, 0.061186, 1.399403, 2.220585, 1.399351, 2.218790, -0.584750]
+    square_cylinder.test_vals = [-2.560840, -1.173495, 0.061186, 1.399403, 2.220585, 1.399351, 2.218790, 0]
     square_cylinder.unsteady  = True
     test_list.append(square_cylinder)
 
@@ -796,7 +796,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882680, -0.215041, 0.023758, -617.450000]
+    ddes_flatplate.test_vals = [-2.714785, -5.882680, -0.215041, 0.023758, 0]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
@@ -1093,7 +1093,7 @@ def main():
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
     p1rad.test_iter = 100
-    p1rad.test_vals = [-7.751309, -7.923059, -2.119084, 0.091733, -44.537000]
+    p1rad.test_vals = [-7.751309, -7.923059, -2.119084, 0.091733, -47.387000]
     test_list.append(p1rad)
 
     # ###############################
@@ -1520,8 +1520,8 @@ def main():
     pywrapper_turb_naca0012_sst.cfg_dir   = "rans/naca0012"
     pywrapper_turb_naca0012_sst.cfg_file  = "turb_NACA0012_sst.cfg"
     pywrapper_turb_naca0012_sst.test_iter = 10
-    pywrapper_turb_naca0012_sst.test_vals         = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, -38.976000]
-    pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, -38.976000]
+    pywrapper_turb_naca0012_sst.test_vals         = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, 0]
+    pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.107132, -15.277740, -6.210248, 1.049757, 0.019249, -3.173936, 0]
     pywrapper_turb_naca0012_sst.command   =  TestCase.Command(exec = "SU2_CFD.py", param = "-f")
     pywrapper_turb_naca0012_sst.timeout   = 3200
     pywrapper_turb_naca0012_sst.tol       = 0.00001
@@ -1534,7 +1534,7 @@ def main():
     pywrapper_square_cylinder.cfg_dir   = "unsteady/square_cylinder"
     pywrapper_square_cylinder.cfg_file  = "turb_square.cfg"
     pywrapper_square_cylinder.test_iter = 3
-    pywrapper_square_cylinder.test_vals = [-2.560840, -1.173495, 0.061186, 1.399403, 2.220585, 1.399351, 2.218790, -0.584750]
+    pywrapper_square_cylinder.test_vals = [-2.560840, -1.173495, 0.061186, 1.399403, 2.220585, 1.399351, 2.218790, 0]
     pywrapper_square_cylinder.command   =  TestCase.Command(exec = "SU2_CFD.py", param = "-f")
     pywrapper_square_cylinder.timeout   = 1600
     pywrapper_square_cylinder.tol       = 0.00001
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 9817708b300..9bca48cd310 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -110,7 +110,7 @@ def main():
     discadj_incomp_cylinder.cfg_dir   = "disc_adj_incomp_navierstokes/cylinder"
     discadj_incomp_cylinder.cfg_file  = "heated_cylinder.cfg"
     discadj_incomp_cylinder.test_iter = 20
-    discadj_incomp_cylinder.test_vals = [20.000000, -2.373367, -2.368305, 0.000000]
+    discadj_incomp_cylinder.test_vals = [20.000000, -2.386480, -2.408986, 0.000000]
     test_list.append(discadj_incomp_cylinder)
 
     #######################################################
@@ -183,7 +183,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals = [-2.122406, 0.693852, 0.000000, -0.869010]
+    discadj_heat.test_vals = [-1.956349, 0.722272, 0.000000, 0.007024]
     test_list.append(discadj_heat)
 
     ###################################

From d476c8ba2b634dadfe39ab56667702abd9c322d5 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 15 Dec 2024 20:36:26 -0800
Subject: [PATCH 09/21] update file regression

---
 .../include/solvers/CFVMFlowSolverBase.inl    |   9 +-
 SU2_CFD/src/solvers/CIncNSSolver.cpp          |  11 +-
 .../radiation/p1adjoint/of_grad_cd.csv.ref    | 100 +++++++++---------
 3 files changed, 59 insertions(+), 61 deletions(-)

diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 40206768e4c..d7a634e8227 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -2357,8 +2357,6 @@ 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();
@@ -2394,8 +2392,8 @@ 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);
 
@@ -2404,7 +2402,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
     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);
       }
     }
@@ -2677,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];
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
index 8db972191c3..e0a4fda96ef 100644
--- a/SU2_CFD/src/solvers/CIncNSSolver.cpp
+++ b/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
diff --git a/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref b/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref
index 8c422aae719..0d3e4c92dbc 100644
--- a/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref
+++ b/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref
@@ -1,51 +1,51 @@
 "VARIABLE"      , "GRADIENT"      , "FINDIFF_STEP"  
-0               , -0.000954218    , 0.001           
-1               , -0.00432384     , 0.001           
-2               , -0.0111485      , 0.001           
-3               , -0.0214656      , 0.001           
-4               , -0.0344212      , 0.001           
-5               , -0.0490221      , 0.001           
-6               , -0.0643278      , 0.001           
-7               , -0.0792041      , 0.001           
-8               , -0.0924318      , 0.001           
-9               , -0.103025       , 0.001           
-10              , -0.110239       , 0.001           
-11              , -0.113383       , 0.001           
-12              , -0.111898       , 0.001           
-13              , -0.105745       , 0.001           
-14              , -0.0956773      , 0.001           
-15              , -0.0830288      , 0.001           
-16              , -0.069201       , 0.001           
-17              , -0.0553022      , 0.001           
-18              , -0.0421316      , 0.001           
-19              , -0.0303252      , 0.001           
-20              , -0.0203976      , 0.001           
-21              , -0.0126387      , 0.001           
-22              , -0.00701749     , 0.001           
-23              , -0.0032361      , 0.001           
-24              , -0.000955756    , 0.001           
-25              , 0.00128953      , 0.001           
-26              , 0.00518424      , 0.001           
-27              , 0.0125712       , 0.001           
-28              , 0.0234942       , 0.001           
-29              , 0.0370801       , 0.001           
-30              , 0.052231        , 0.001           
-31              , 0.06787         , 0.001           
-32              , 0.0827503       , 0.001           
-33              , 0.0955927       , 0.001           
-34              , 0.105415        , 0.001           
-35              , 0.11154         , 0.001           
-36              , 0.113399        , 0.001           
-37              , 0.110583        , 0.001           
-38              , 0.103215        , 0.001           
-39              , 0.0921957       , 0.001           
-40              , 0.0789703       , 0.001           
-41              , 0.0650043       , 0.001           
-42              , 0.0514107       , 0.001           
-43              , 0.0389322       , 0.001           
-44              , 0.0280921       , 0.001           
-45              , 0.0192441       , 0.001           
-46              , 0.0124768       , 0.001           
-47              , 0.00752972      , 0.001           
-48              , 0.00391057      , 0.001           
-49              , 0.00131276      , 0.001           
+0               , -0.000958205    , 0.001           
+1               , -0.00432733     , 0.001           
+2               , -0.0111443      , 0.001           
+3               , -0.0214515      , 0.001           
+4               , -0.0343997      , 0.001           
+5               , -0.0489968      , 0.001           
+6               , -0.0643023      , 0.001           
+7               , -0.0791818      , 0.001           
+8               , -0.092416       , 0.001           
+9               , -0.103017       , 0.001           
+10              , -0.11024        , 0.001           
+11              , -0.113392       , 0.001           
+12              , -0.111913       , 0.001           
+13              , -0.105766       , 0.001           
+14              , -0.0957015      , 0.001           
+15              , -0.0830517      , 0.001           
+16              , -0.0692182      , 0.001           
+17              , -0.0553101      , 0.001           
+18              , -0.042129       , 0.001           
+19              , -0.0303137      , 0.001           
+20              , -0.0203813      , 0.001           
+21              , -0.0126226      , 0.001           
+22              , -0.00700441     , 0.001           
+23              , -0.00322726     , 0.001           
+24              , -0.000951904    , 0.001           
+25              , 0.00128803      , 0.001           
+26              , 0.0051774       , 0.001           
+27              , 0.0125547       , 0.001           
+28              , 0.0234674       , 0.001           
+29              , 0.0370465       , 0.001           
+30              , 0.0521956       , 0.001           
+31              , 0.0678379       , 0.001           
+32              , 0.0827257       , 0.001           
+33              , 0.0955787       , 0.001           
+34              , 0.105412        , 0.001           
+35              , 0.111548        , 0.001           
+36              , 0.113415        , 0.001           
+37              , 0.110605        , 0.001           
+38              , 0.103242        , 0.001           
+39              , 0.0922244       , 0.001           
+40              , 0.0789963       , 0.001           
+41              , 0.0650234       , 0.001           
+42              , 0.0514199       , 0.001           
+43              , 0.0389312       , 0.001           
+44              , 0.0280831       , 0.001           
+45              , 0.0192316       , 0.001           
+46              , 0.0124655       , 0.001           
+47              , 0.00752232      , 0.001           
+48              , 0.0039072       , 0.001           
+49              , 0.00131195      , 0.001           

From bc59826ee2a4ac70a8fbaac45a966f61d61dd949 Mon Sep 17 00:00:00 2001
From: emaberman <emaberman@gmail.com>
Date: Wed, 18 Dec 2024 20:38:55 +0200
Subject: [PATCH 10/21] Fix regression testing

---
 TestCases/hybrid_regression.py      | 10 +++++-----
 TestCases/parallel_regression.py    | 28 ++++++++++++++--------------
 TestCases/parallel_regression_AD.py |  4 ++--
 TestCases/serial_regression.py      | 16 ++++++++--------
 TestCases/serial_regression_AD.py   |  2 +-
 TestCases/tutorials.py              |  6 +++---
 TestCases/vandv.py                  |  4 ++--
 7 files changed, 35 insertions(+), 35 deletions(-)

diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index fe2634fb2a5..e5b0a847e61 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -157,7 +157,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, -80602.000000]
+    rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, -80603.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -181,7 +181,7 @@ def main():
     turb_flatplate.cfg_dir   = "rans/flatplate"
     turb_flatplate.cfg_file  = "turb_SA_flatplate.cfg"
     turb_flatplate.test_iter = 20
-    turb_flatplate.test_vals = [-4.312826, -6.736053, -0.187467, 0.057454]
+    turb_flatplate.test_vals = [-4.312725, -6.737976, -0.187467, 0.057468]
     test_list.append(turb_flatplate)
 
     # ONERA M6 Wing
@@ -238,7 +238,7 @@ def main():
     propeller.cfg_dir   = "rans/propeller"
     propeller.cfg_file  = "propeller.cfg"
     propeller.test_iter = 10
-    propeller.test_vals = [-3.389724, -8.409223, 0.000048, 0.056344]
+    propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344]
     test_list.append(propeller)
 
     #######################################
@@ -403,7 +403,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788405, -11.039465, 0.000008, 0.309490]
+    inc_turb_naca0012.test_vals = [-4.788405, -11.040560, 0.000008, 0.309505]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST
@@ -503,7 +503,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882681, -0.215041, 0.023758, -617.440000]
+    ddes_flatplate.test_vals = [-2.714786, -5.882652, -0.215041, 0.023758, -617.470000]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index 5462560119e..f19b323f761 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -357,7 +357,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.004689, -5.265782, 0.809463, 0.062016, -80576.000000]
+    rae2822_sa.test_vals = [-2.004689, -5.265797, 0.809463, 0.062016, -80577.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -381,7 +381,7 @@ def main():
     turb_flatplate.cfg_dir   = "rans/flatplate"
     turb_flatplate.cfg_file  = "turb_SA_flatplate.cfg"
     turb_flatplate.test_iter = 20
-    turb_flatplate.test_vals = [-4.293562, -6.728390, -0.187643, 0.057686]
+    turb_flatplate.test_vals = [-4.293470, -6.730438, -0.187644, 0.057700]
     test_list.append(turb_flatplate)
 
     # Flat plate (compressible) with species inlet
@@ -389,7 +389,7 @@ def main():
     turb_flatplate_species.cfg_dir   = "rans/flatplate"
     turb_flatplate_species.cfg_file  = "turb_SA_flatplate_species.cfg"
     turb_flatplate_species.test_iter = 20
-    turb_flatplate_species.test_vals = [-4.243136, -0.634954, -1.706721, 1.231179, -3.266212, 9.000000, -6.632862, 5.000000, -6.979197, 10.000000, -6.007240, 0.996237, 0.996237]
+    turb_flatplate_species.test_vals = [-4.243064, -0.634797, -1.706652, 1.231264, -3.266203, 9.000000, -6.632972, 5.000000, -6.985977, 10.000000, -6.007208, 0.996237, 0.996237]
     test_list.append(turb_flatplate_species)
 
     # Flat plate SST compressibility correction Wilcox
@@ -431,7 +431,7 @@ def main():
     turb_oneram6_nk.cfg_dir   = "rans/oneram6"
     turb_oneram6_nk.cfg_file  = "turb_ONERAM6_nk.cfg"
     turb_oneram6_nk.test_iter = 20
-    turb_oneram6_nk.test_vals = [-4.843772, -4.444210, -11.473964, 0.216337, 0.049646, 5.000000, -0.613234, 23.568000]
+    turb_oneram6_nk.test_vals = [-4.851388, -4.457414, -11.468726, 0.217228, 0.049043, 5.000000, -0.533763, 23.567000]
     turb_oneram6_nk.timeout   = 600
     turb_oneram6_nk.tol       = 0.0001
     test_list.append(turb_oneram6_nk)
@@ -441,7 +441,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.094721, -14.685268, 1.057665, 0.022971, 20.000000, -1.692925, 20.000000, -4.037672, -44.871000]
+    turb_naca0012_sa.test_vals = [-12.094695, -14.685268, 1.057665, 0.022971, 20.000000, -1.692967, 20.000000, -4.037673, -44.871000]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
 
@@ -507,7 +507,7 @@ def main():
     propeller.cfg_dir   = "rans/propeller"
     propeller.cfg_file  = "propeller.cfg"
     propeller.test_iter = 10
-    propeller.test_vals = [-3.389724, -8.409223, 0.000048, 0.056344]
+    propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344]
     propeller.timeout   = 3200
     test_list.append(propeller)
 
@@ -638,7 +638,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788596, -11.039529, -0.000002, 0.309504]
+    inc_turb_naca0012.test_vals = [-4.788595, -11.040625, -0.000002, 0.309519]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST
@@ -724,7 +724,7 @@ def main():
     turbmod_sa_bsl_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_bsl_rae2822.cfg_file  = "turb_SA_BSL_RAE2822.cfg"
     turbmod_sa_bsl_rae2822.test_iter = 20
-    turbmod_sa_bsl_rae2822.test_vals = [-2.004689, 0.742306, 0.497308, -5.265782, 0.809463, 0.062016]
+    turbmod_sa_bsl_rae2822.test_vals = [-2.004689, 0.742306, 0.497308, -5.265797, 0.809463, 0.062016]
     test_list.append(turbmod_sa_bsl_rae2822)
 
     # SA Negative
@@ -732,7 +732,7 @@ def main():
     turbmod_sa_neg_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_neg_rae2822.cfg_file  = "turb_SA_NEG_RAE2822.cfg"
     turbmod_sa_neg_rae2822.test_iter = 10
-    turbmod_sa_neg_rae2822.test_vals         = [-1.466238, 3.169232, 2.756518, 4.722767, 1.120253, 0.378675, -83444.000000]
+    turbmod_sa_neg_rae2822.test_vals         = [-1.204800, 1.611685, 1.349330, 1.489602, 1.263603, 0.466487, -83309.000000]
     turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.359612, 1.493629, 1.218367, -1.441703, 1.248499, 0.457987, -86467.000000]
     test_list.append(turbmod_sa_neg_rae2822)
 
@@ -741,7 +741,7 @@ def main():
     turbmod_sa_comp_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_comp_rae2822.cfg_file  = "turb_SA_COMP_RAE2822.cfg"
     turbmod_sa_comp_rae2822.test_iter = 20
-    turbmod_sa_comp_rae2822.test_vals = [-2.004687, 0.742304, 0.497309, -5.266067, 0.809467, 0.062029]
+    turbmod_sa_comp_rae2822.test_vals = [-2.004687, 0.742304, 0.497309, -5.266084, 0.809467, 0.062029]
     test_list.append(turbmod_sa_comp_rae2822)
 
     # SA Edwards
@@ -765,7 +765,7 @@ def main():
     turbmod_sa_qcr_rae2822.cfg_dir   = "turbulence_models/sa/rae2822"
     turbmod_sa_qcr_rae2822.cfg_file  = "turb_SA_QCR_RAE2822.cfg"
     turbmod_sa_qcr_rae2822.test_iter = 20
-    turbmod_sa_qcr_rae2822.test_vals = [-2.004793, 0.742353, 0.497315, -5.265962, 0.807841, 0.062027]
+    turbmod_sa_qcr_rae2822.test_vals = [-2.004793, 0.742353, 0.497315, -5.265977, 0.807841, 0.062027]
     test_list.append(turbmod_sa_qcr_rae2822)
 
     ############################
@@ -1001,7 +1001,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882679, -0.215041, 0.023758, -617.450000]
+    ddes_flatplate.test_vals = [-2.714786, -5.882637, -0.215041, 0.023758, -617.470000]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
@@ -1010,7 +1010,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-3.004012, -6.876220, 1.487882, 0.421861]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.004011, -6.876250, 1.487888, 0.421869]
     unst_inc_turb_naca0015_sa.unsteady  = True
     test_list.append(unst_inc_turb_naca0015_sa)
 
@@ -1530,7 +1530,7 @@ def main():
     species2_primitiveVenturi_mixingmodel_viscosity.cfg_dir   = "species_transport/venturi_primitive_3species"
     species2_primitiveVenturi_mixingmodel_viscosity.cfg_file  = "species2_primitiveVenturi_mixingmodel_viscosity.cfg"
     species2_primitiveVenturi_mixingmodel_viscosity.test_iter = 50
-    species2_primitiveVenturi_mixingmodel_viscosity.test_vals = [-4.857397, -3.646605, -3.737462, -7.602922, -5.008846, 5.000000, -1.756226, 5.000000, -3.163353, 5.000000, -2.189723, 2.476808, 0.976999, 0.609280, 0.890529]
+    species2_primitiveVenturi_mixingmodel_viscosity.test_vals = [-4.857405, -3.646534, -3.737422, -7.602756, -5.008835, 5.000000, -1.756256, 5.000000, -3.163548, 5.000000, -2.189690, 2.476807, 0.977000, 0.609279, 0.890528]
     test_list.append(species2_primitiveVenturi_mixingmodel_viscosity)
 
     # 2 species (1 eq) primitive venturi mixing using mixing model including heat capacity and mass diffusivity
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 7d0f0d10049..c549f753c0f 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -83,7 +83,7 @@ def main():
     discadj_rans_naca0012_sa.cfg_dir   = "disc_adj_rans/naca0012"
     discadj_rans_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     discadj_rans_naca0012_sa.test_iter = 10
-    discadj_rans_naca0012_sa.test_vals = [-2.996963, -0.196020, 0.000004, -0.000000, 5.000000, -3.430615, 5.000000, -7.411396]
+    discadj_rans_naca0012_sa.test_vals = [-2.996963, -0.196020, 0.000004, -0.000000, 5.000000, -3.430615, 5.000000, -7.411381]
     test_list.append(discadj_rans_naca0012_sa)
 
     # Adjoint turbulent NACA0012 SST
@@ -256,7 +256,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals         = [-2.117791, 0.784475, 0.000000, -0.574700]
+    discadj_heat.test_vals         = [-2.117792, 0.784117, 0.000000, -0.574700]
     discadj_heat.test_vals_aarch64 = [-2.226539, 0.605868, 0.000000, -6.256400]
     test_list.append(discadj_heat)
 
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 10312c92cc3..4aa5f883481 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -218,7 +218,7 @@ def main():
     rae2822_sa.cfg_dir   = "rans/rae2822"
     rae2822_sa.cfg_file  = "turb_SA_RAE2822.cfg"
     rae2822_sa.test_iter = 20
-    rae2822_sa.test_vals = [-2.020123, -5.269324, 0.807147, 0.060499, -80602.000000]
+    rae2822_sa.test_vals = [-2.020123, -5.269339, 0.807147, 0.060499, -80603.000000]
     test_list.append(rae2822_sa)
 
     # RAE2822 SST
@@ -242,7 +242,7 @@ def main():
     turb_flatplate.cfg_dir   = "rans/flatplate"
     turb_flatplate.cfg_file  = "turb_SA_flatplate.cfg"
     turb_flatplate.test_iter = 20
-    turb_flatplate.test_vals = [-4.312826, -6.736053, -0.187467, 0.057454]
+    turb_flatplate.test_vals = [-4.312724, -6.737976, -0.187467, 0.057468]
     test_list.append(turb_flatplate)
 
     # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SST
@@ -258,7 +258,7 @@ def main():
     turb_wallfunction_flatplate_sa.cfg_dir   = "wallfunctions/flatplate/compressible_SA"
     turb_wallfunction_flatplate_sa.cfg_file  = "turb_SA_flatplate.cfg"
     turb_wallfunction_flatplate_sa.test_iter = 10
-    turb_wallfunction_flatplate_sa.test_vals = [-4.460537, -2.033605, -2.117564, 0.889626, -5.381903, 10.000000, -1.517487, 0.034212, 0.002636]
+    turb_wallfunction_flatplate_sa.test_vals = [-4.460657, -2.033641, -2.118149, 0.889562, -5.382249, 10.000000, -1.517453, 0.034213, 0.002636]
     test_list.append(turb_wallfunction_flatplate_sa)
 
     # ONERA M6 Wing
@@ -275,7 +275,7 @@ def main():
     turb_naca0012_sa.cfg_dir   = "rans/naca0012"
     turb_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     turb_naca0012_sa.test_iter = 5
-    turb_naca0012_sa.test_vals = [-12.091696, -14.685322, 1.057665, 0.022971, 20.000000, -2.686306, 20.000000, -4.459916, -44.871000]
+    turb_naca0012_sa.test_vals = [-12.091778, -14.685322, 1.057665, 0.022971, 20.000000, -2.686203, 20.000000, -4.459916, -44.871000]
     turb_naca0012_sa.timeout   = 3200
     test_list.append(turb_naca0012_sa)
 
@@ -322,7 +322,7 @@ def main():
     propeller.cfg_dir   = "rans/propeller"
     propeller.cfg_file  = "propeller.cfg"
     propeller.test_iter = 10
-    propeller.test_vals = [-3.389724, -8.409223, 0.000048, 0.056344]
+    propeller.test_vals = [-3.389724, -8.409502, 0.000048, 0.056344]
     propeller.timeout   = 3200
     test_list.append(propeller)
 
@@ -441,7 +441,7 @@ def main():
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
     inc_turb_naca0012.test_iter = 20
-    inc_turb_naca0012.test_vals = [-4.788496, -11.039482, 0.000023, 0.309488]
+    inc_turb_naca0012.test_vals = [-4.788495, -11.040578, 0.000023, 0.309503]
     test_list.append(inc_turb_naca0012)
 
     # NACA0012, SST_SUST
@@ -796,7 +796,7 @@ def main():
     ddes_flatplate.cfg_dir   = "ddes/flatplate"
     ddes_flatplate.cfg_file  = "ddes_flatplate.cfg"
     ddes_flatplate.test_iter = 10
-    ddes_flatplate.test_vals = [-2.714785, -5.882680, -0.215041, 0.023758, -617.450000]
+    ddes_flatplate.test_vals = [-2.714786, -5.882646, -0.215041, 0.023758, -617.470000]
     ddes_flatplate.unsteady  = True
     test_list.append(ddes_flatplate)
 
@@ -805,7 +805,7 @@ def main():
     unst_inc_turb_naca0015_sa.cfg_dir   = "unsteady/pitching_naca0015_rans_inc"
     unst_inc_turb_naca0015_sa.cfg_file  = "config_incomp_turb_sa.cfg"
     unst_inc_turb_naca0015_sa.test_iter = 1
-    unst_inc_turb_naca0015_sa.test_vals = [-3.007635, -6.879778, 1.445293, 0.419274]
+    unst_inc_turb_naca0015_sa.test_vals = [-3.007635, -6.879809, 1.445300, 0.419281]
     unst_inc_turb_naca0015_sa.unsteady  = True
     test_list.append(unst_inc_turb_naca0015_sa)
 
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 9817708b300..8dc7ca3704c 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -183,7 +183,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals = [-2.122406, 0.693852, 0.000000, -0.869010]
+    discadj_heat.test_vals = [-2.122406, 0.693159, 0.000000, -0.869010]
     test_list.append(discadj_heat)
 
     ###################################
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
index eb018f65e73..e819ff246c1 100644
--- a/TestCases/tutorials.py
+++ b/TestCases/tutorials.py
@@ -228,7 +228,7 @@ def main():
     tutorial_trans_e387_sa.cfg_dir    = "../Tutorials/compressible_flow/Transitional_Airfoil/Langtry_and_Menter/E387"
     tutorial_trans_e387_sa.cfg_file   = "transitional_SA_LM_model_ConfigFile.cfg"
     tutorial_trans_e387_sa.test_iter  = 20
-    tutorial_trans_e387_sa.test_vals  = [-6.527027, -5.081558, -0.795261, 1.022606, 0.150125, 2.000000, -9.580659]
+    tutorial_trans_e387_sa.test_vals  = [-6.527027, -5.081541, -0.795261, 1.022607, 0.150175, 2.000000, -9.580660]
     tutorial_trans_e387_sa.no_restart = True
     test_list.append(tutorial_trans_e387_sa)
 
@@ -263,7 +263,7 @@ def main():
     tutorial_unst_naca0012.cfg_dir       = "../Tutorials/compressible_flow/Unsteady_NACA0012"
     tutorial_unst_naca0012.cfg_file      = "unsteady_naca0012.cfg"
     tutorial_unst_naca0012.test_iter     = 520
-    tutorial_unst_naca0012.test_vals         = [520.000000, 0.000000, -5.291711, 0.000000, 0.305248, 0.810326, 0.001814, 0.006573]
+    tutorial_unst_naca0012.test_vals         = [520.000000, 0.000000, -5.290694, 0.000000, 0.317272, 0.820972, 0.002144, 0.012805]
     tutorial_unst_naca0012.test_vals_aarch64 = [520.000000, 0.000000, -5.298777, 0.000000, 0.288956, 0.736706, 0.002419, 0.007134]
     tutorial_unst_naca0012.unsteady      = True
     test_list.append(tutorial_unst_naca0012)
@@ -273,7 +273,7 @@ def main():
     propeller_var_load.cfg_dir   = "../Tutorials/compressible_flow/ActuatorDisk_VariableLoad"
     propeller_var_load.cfg_file  = "propeller_variable_load.cfg"
     propeller_var_load.test_iter = 20
-    propeller_var_load.test_vals = [-1.830276, -4.535127, -0.000323, 0.171623]
+    propeller_var_load.test_vals = [-1.830257, -4.535041, -0.000323, 0.171647]
     propeller_var_load.timeout   = 3200
     test_list.append(propeller_var_load)
 
diff --git a/TestCases/vandv.py b/TestCases/vandv.py
index 01359163885..f9da6969994 100644
--- a/TestCases/vandv.py
+++ b/TestCases/vandv.py
@@ -45,8 +45,8 @@ def main():
     p30n30.cfg_dir   = "vandv/rans/30p30n"
     p30n30.cfg_file  = "config.cfg"
     p30n30.test_iter = 20
-    p30n30.test_vals         = [-10.582171, -10.106603, -10.474926, -10.182536, -12.679336, 0.052181, 2.829820, 1.318613, -0.221214]
-    p30n30.test_vals_aarch64 = [-10.582171, -10.106603, -10.474926, -10.182536, -12.679336, 0.052181, 2.829820, 1.318613, -0.221214]
+    p30n30.test_vals         = [-10.582183, -10.106601, -10.474910, -10.182549, -12.679336, 0.052181, 2.829820, 1.318613, -0.221374]
+    p30n30.test_vals_aarch64 = [-10.582183, -10.106601, -10.474910, -10.182549, -12.679336, 0.052181, 2.829820, 1.318613, -0.221374]
     test_list.append(p30n30)
 
     # flat plate - sst-v1994m

From 1b8fad3391a687f748507b5fe017d42f529d9867 Mon Sep 17 00:00:00 2001
From: emaberman <emaberman@gmail.com>
Date: Wed, 18 Dec 2024 21:57:19 +0200
Subject: [PATCH 11/21] fixes

---
 TestCases/serial_regression_AD.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 8dc7ca3704c..71d4ccf90ea 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -78,7 +78,7 @@ def main():
     discadj_rans_naca0012_sa.cfg_dir   = "disc_adj_rans/naca0012"
     discadj_rans_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     discadj_rans_naca0012_sa.test_iter = 10
-    discadj_rans_naca0012_sa.test_vals = [-2.996976, -0.196055, 0.000004, -0.000000, 5.000000, -3.971736, 5.000000, -10.464337]
+    discadj_rans_naca0012_sa.test_vals = [-2.996976, -0.196055, 0.000004, -0.000000, 5.000000, -3.971736, 5.000000, -10.464319v]
     test_list.append(discadj_rans_naca0012_sa)
 
     # Adjoint turbulent NACA0012 SST

From cea32627c0e7476b2c9688ea12c07d66e36667fa Mon Sep 17 00:00:00 2001
From: emaberman <emaberman@gmail.com>
Date: Wed, 18 Dec 2024 22:03:53 +0200
Subject: [PATCH 12/21] fixes

---
 TestCases/serial_regression_AD.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 71d4ccf90ea..47cb87fc98e 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -78,7 +78,7 @@ def main():
     discadj_rans_naca0012_sa.cfg_dir   = "disc_adj_rans/naca0012"
     discadj_rans_naca0012_sa.cfg_file  = "turb_NACA0012_sa.cfg"
     discadj_rans_naca0012_sa.test_iter = 10
-    discadj_rans_naca0012_sa.test_vals = [-2.996976, -0.196055, 0.000004, -0.000000, 5.000000, -3.971736, 5.000000, -10.464319v]
+    discadj_rans_naca0012_sa.test_vals = [-2.996976, -0.196055, 0.000004, -0.000000, 5.000000, -3.971736, 5.000000, -10.464319]
     test_list.append(discadj_rans_naca0012_sa)
 
     # Adjoint turbulent NACA0012 SST

From 46405cdd5ef8dab7ac6676052c7a38d6b4263423 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sat, 21 Dec 2024 14:48:10 -0800
Subject: [PATCH 13/21] add thermal expansion to CConfig

---
 Common/include/CConfig.hpp                       | 11 +++++++++--
 Common/src/CConfig.cpp                           | 14 +++++++++++---
 .../numerics/elasticity/CFEAElasticity.hpp       |  2 ++
 .../src/numerics/elasticity/CFEAElasticity.cpp   | 13 ++++++-------
 .../numerics/elasticity/CFEALinearElasticity.cpp |  7 +++++++
 .../src/numerics/elasticity/nonlinear_models.cpp |  2 +-
 TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg  | 16 ++++++++--------
 7 files changed, 44 insertions(+), 21 deletions(-)

diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 0c32f044288..9eb6ef74e71 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1013,10 +1013,12 @@ class CConfig {
   su2double *Int_Coeffs;         /*!< \brief Time integration coefficients for structural method. */
   unsigned short nElasticityMod, /*!< \brief Number of different values for the elasticity modulus. */
   nPoissonRatio,                    /*!< \brief Number of different values for the Poisson ratio modulus. */
-  nMaterialDensity;                 /*!< \brief Number of different values for the Material density. */
+  nMaterialDensity,                 /*!< \brief Number of different values for the Material density. */
+  nMaterialThermalExpansion;        /*!< \brief Number of different values for thermal expansion coefficient. */
   su2double *ElasticityMod,         /*!< \brief Value of the elasticity moduli. */
   *PoissonRatio,                    /*!< \brief Value of the Poisson ratios. */
-  *MaterialDensity;                 /*!< \brief Value of the Material densities. */
+  *MaterialDensity,                 /*!< \brief Value of the Material densities. */
+  *MaterialThermalExpansion;        /*!< \brief Value of the thermal expansion coefficients. */
   unsigned short nElectric_Field,   /*!< \brief Number of different values for the electric field in the membrane. */
   nDim_Electric_Field;              /*!< \brief Dimensionality of the problem. */
   unsigned short nDim_RefNode;      /*!< \brief Dimensionality of the vector . */
@@ -2389,6 +2391,11 @@ class CConfig {
    */
   su2double GetMaterialDensity(unsigned short id_val) const { return MaterialDensity[id_val]; }
 
+  /*!
+   * \brief Get the thermal expansion coefficient.
+   */
+  su2double GetMaterialThermalExpansion(unsigned short id_val) const { return MaterialThermalExpansion[id_val]; }
+
   /*!
    * \brief Compressibility/incompressibility of the solids analysed using the structural solver.
    * \return Compressible or incompressible.
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index c5bd653b04a..90af5948d25 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -2440,6 +2440,8 @@ void CConfig::SetConfig_Options() {
   addDoubleListOption("POISSON_RATIO", nPoissonRatio, PoissonRatio);
   /* DESCRIPTION: Material density */
   addDoubleListOption("MATERIAL_DENSITY", nMaterialDensity, MaterialDensity);
+  /* DESCRIPTION: Material thermal expansion coefficient */
+  addDoubleListOption("MATERIAL_THERMAL_EXPANSION_COEFF", nMaterialThermalExpansion, MaterialThermalExpansion);
   /* DESCRIPTION: Knowles B constant */
   addDoubleOption("KNOWLES_B", Knowles_B, 1.0);
   /* DESCRIPTION: Knowles N constant */
@@ -4834,9 +4836,15 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854;
   }
 
-  if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity) {
-    SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, and MATERIAL_DENSITY need to have the same number "
-                   "of entries (the number of materials).", CURRENT_FUNCTION);
+  if (nMaterialThermalExpansion == 0) {
+    nMaterialThermalExpansion = 1;
+    MaterialThermalExpansion = new su2double[1]();
+  }
+
+  if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity ||
+      nElasticityMod != nMaterialThermalExpansion) {
+    SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, MATERIAL_DENSITY, and THERMAL_EXPANSION_COEFF need "
+                   "to have the same number of entries (the number of materials).", CURRENT_FUNCTION);
   }
 
   if (nElectric_Constant == 0) {
diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
index 11470fd0fba..43d5ebd5a63 100644
--- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
+++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
@@ -54,6 +54,7 @@ class CFEAElasticity : public CNumerics {
   su2double Nu        = 0.0;              /*!< \brief Aux. variable, Poisson's ratio. */
   su2double Rho_s     = 0.0;              /*!< \brief Aux. variable, Structural density. */
   su2double Rho_s_DL  = 0.0;              /*!< \brief Aux. variable, Structural density (for dead loads). */
+  su2double Alpha     = 0.0;              /*!< \brief Aux. variable, thermal expansion coefficient. */
 
   su2double Mu        = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
   su2double Lambda    = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
@@ -63,6 +64,7 @@ class CFEAElasticity : public CNumerics {
   su2double *Nu_i     = nullptr;          /*!< \brief Poisson's ratio. */
   su2double *Rho_s_i  = nullptr;          /*!< \brief Structural density. */
   su2double *Rho_s_DL_i = nullptr;        /*!< \brief Structural density (for dead loads). */
+  su2double *Alpha_i  = nullptr;          /*!< \brief Thermal expansion coefficient. */
 
   su2double **Ba_Mat = nullptr;           /*!< \brief Matrix B for node a - Auxiliary. */
   su2double **Bb_Mat = nullptr;           /*!< \brief Matrix B for node b - Auxiliary. */
diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
index 44014591204..b91be37223e 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
@@ -44,19 +44,16 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   const auto nProp = config->GetnElasticityMat();
 
   E_i = new su2double[nProp];
-  for (iVar = 0; iVar < nProp; iVar++)
-    E_i[iVar] = config->GetElasticyMod(iVar);
-
   Nu_i = new su2double[nProp];
-  for (iVar = 0; iVar < nProp; iVar++)
-    Nu_i[iVar] = config->GetPoissonRatio(iVar);
-
   Rho_s_i = new su2double[nProp];     // For inertial effects
   Rho_s_DL_i = new su2double[nProp];  // For dead loads
-
+  Alpha_i = new su2double[nProp];
   for (iVar = 0; iVar < nProp; iVar++) {
+    E_i[iVar] = config->GetElasticyMod(iVar);
+    Nu_i[iVar] = config->GetPoissonRatio(iVar);
     Rho_s_DL_i[iVar] = config->GetMaterialDensity(iVar);
     Rho_s_i[iVar] = pseudo_static ? 0.0 : config->GetMaterialDensity(iVar);
+    Alpha_i[iVar] = config->GetMaterialThermalExpansion(iVar);
   }
 
   // Initialization
@@ -64,6 +61,7 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   Nu  = Nu_i[0];
   Rho_s = Rho_s_i[0];
   Rho_s_DL = Rho_s_DL_i[0];
+  Alpha = Alpha_i[0];
 
   Compute_Lame_Parameters();
 
@@ -302,6 +300,7 @@ void CFEAElasticity::SetElement_Properties(const CElement *element, const CConfi
   Nu  = Nu_i[element->Get_iProp()];
   Rho_s = Rho_s_i[element->Get_iProp()];
   Rho_s_DL = Rho_s_DL_i[element->Get_iProp()];
+  Alpha = Alpha_i[element->Get_iProp()];
 
   switch (config->GetDV_FEA()) {
     case YOUNG_MODULUS:
diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
index a6763bed0b8..de7ad86219f 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
@@ -99,6 +99,13 @@ void CFEALinearElasticity::Compute_Tangent_Matrix(CElement *element, const CConf
 
     for (iNode = 0; iNode < nNode; iNode++) {
 
+      su2double KAux_t_a[3] = {0.0};
+      const su2double thermalStress = -E * 2e-4 / (1 - 2 * Nu);
+      for (iVar = 0; iVar < nDim; iVar++) {
+        KAux_t_a[iVar] += Weight * thermalStress * GradNi_Ref_Mat[iNode][iVar] * Jac_X;
+      }
+      element->Add_Kt_a(iNode, KAux_t_a);
+
       if (nDim == 2) {
         Ba_Mat[0][0] = GradNi_Ref_Mat[iNode][0];
         Ba_Mat[1][1] = GradNi_Ref_Mat[iNode][1];
diff --git a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
index d21fc228300..55d48a29010 100644
--- a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
+++ b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
@@ -112,7 +112,7 @@ void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfi
   for (iVar = 0; iVar < 3; iVar++) {
     for (jVar = 0; jVar < 3; jVar++) {
       su2double dij = deltaij(iVar,jVar);
-      Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + Lambda_J * log(J_F) * dij;
+      Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) - E * 2e-4 / (1 - 2 * Nu)) * dij;
     }
   }
 
diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
index a6c9cde3c84..128958ecb19 100644
--- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
+++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
@@ -15,13 +15,13 @@ MESH_FILENAME= meshBeam_3d.su2
 ELASTICITY_MODULUS=3E7
 POISSON_RATIO=0.3
 MATERIAL_DENSITY=7854
-MARKER_CLAMPED = ( left , right )
-MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0)
-MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0)
-LINEAR_SOLVER= FGMRES
-LINEAR_SOLVER_PREC= LU_SGS
-LINEAR_SOLVER_ERROR= 1E-8
-LINEAR_SOLVER_ITER= 500
+MARKER_CLAMPED = ( left )
+MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0)
+%MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0)
+LINEAR_SOLVER= CONJUGATE_GRADIENT
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-9
+LINEAR_SOLVER_ITER= 1000
 MESH_FORMAT= SU2
 TABULAR_FORMAT= CSV
 CONV_FILENAME= history_beam
@@ -29,4 +29,4 @@ VOLUME_FILENAME= beam
 RESTART_FILENAME= restart_beam.dat
 SOLUTION_FILENAME= restart_beam.dat
 OUTPUT_WRT_FREQ= 1
-INNER_ITER=1
+INNER_ITER=10

From aa88aba520bae2858f6980d9cd01430d8131fdc1 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sat, 21 Dec 2024 14:49:05 -0800
Subject: [PATCH 14/21] error if interpolating takes too long

---
 SU2_CFD/src/solvers/CSolver.cpp | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index e446c3f1133..2c6322392c7 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/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;
   }
 }

From 0fc99e944651233ba6ffb75fd2ca860c660f4bb6 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sat, 21 Dec 2024 18:17:59 -0800
Subject: [PATCH 15/21] re-update

---
 TestCases/parallel_regression_AD.py | 2 +-
 TestCases/serial_regression_AD.py   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
index 0cb0d452d57..f54849610f4 100644
--- a/TestCases/parallel_regression_AD.py
+++ b/TestCases/parallel_regression_AD.py
@@ -256,7 +256,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals         = [-2.117792, 0.784117, 0.000000, -0.574700]
+    discadj_heat.test_vals         = [-1.840134, 0.750337, 0.000000, 0.006760]
     discadj_heat.test_vals_aarch64 = [-2.226539, 0.605868, 0.000000, -6.256400]
     test_list.append(discadj_heat)
 
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
index 1a1267e7306..620a3b20e8f 100644
--- a/TestCases/serial_regression_AD.py
+++ b/TestCases/serial_regression_AD.py
@@ -183,7 +183,7 @@ def main():
     discadj_heat.cfg_dir   = "disc_adj_heat"
     discadj_heat.cfg_file  = "disc_adj_heat.cfg"
     discadj_heat.test_iter = 10
-    discadj_heat.test_vals = [-1.956349, 0.722272, 0.000000, 0.007024]
+    discadj_heat.test_vals = [-1.956346, 0.721746, 0.000000, 0.007024]
     test_list.append(discadj_heat)
 
     ###################################

From 458e73224ccdee72f99d79635745197b879c4973 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 22 Dec 2024 18:00:05 -0800
Subject: [PATCH 16/21] thermal stress term, not coupled to heat solver

---
 Common/include/CConfig.hpp                    |  8 +++-
 Common/include/geometry/elements/CElement.hpp | 41 ++++++++++++++---
 Common/src/CConfig.cpp                        |  2 +
 Common/src/geometry/elements/CElement.cpp     |  2 +
 Common/src/geometry/elements/CQUAD4.cpp       |  1 -
 .../numerics/elasticity/CFEAElasticity.hpp    | 11 +++--
 .../elasticity/CFEANonlinearElasticity.hpp    |  3 +-
 .../numerics/elasticity/nonlinear_models.hpp  |  8 ++--
 .../numerics/elasticity/CFEAElasticity.cpp    | 46 ++++++-------------
 .../elasticity/CFEALinearElasticity.cpp       | 10 ++--
 .../elasticity/CFEANonlinearElasticity.cpp    |  6 +--
 .../numerics/elasticity/nonlinear_models.cpp  | 26 ++++++-----
 SU2_CFD/src/solvers/CFEASolver.cpp            |  8 ++++
 .../fea_fsi/StatBeam_3d/configBeam_3d.cfg     |  8 ++--
 14 files changed, 110 insertions(+), 70 deletions(-)

diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 9eb6ef74e71..d91a6011611 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1018,7 +1018,8 @@ class CConfig {
   su2double *ElasticityMod,         /*!< \brief Value of the elasticity moduli. */
   *PoissonRatio,                    /*!< \brief Value of the Poisson ratios. */
   *MaterialDensity,                 /*!< \brief Value of the Material densities. */
-  *MaterialThermalExpansion;        /*!< \brief Value of the thermal expansion coefficients. */
+  *MaterialThermalExpansion,        /*!< \brief Value of the thermal expansion coefficients. */
+  MaterialReferenceTemperature;     /*!< \brief Value of the reference temperature for thermal expansion. */
   unsigned short nElectric_Field,   /*!< \brief Number of different values for the electric field in the membrane. */
   nDim_Electric_Field;              /*!< \brief Dimensionality of the problem. */
   unsigned short nDim_RefNode;      /*!< \brief Dimensionality of the vector . */
@@ -2396,6 +2397,11 @@ class CConfig {
    */
   su2double GetMaterialThermalExpansion(unsigned short id_val) const { return MaterialThermalExpansion[id_val]; }
 
+  /*!
+   * \brief Temperature at which there is no stress from thermal expansion.
+   */
+  su2double GetMaterialReferenceTemperature() const { return MaterialReferenceTemperature; }
+
   /*!
    * \brief Compressibility/incompressibility of the solids analysed using the structural solver.
    * \return Compressible or incompressible.
diff --git a/Common/include/geometry/elements/CElement.hpp b/Common/include/geometry/elements/CElement.hpp
index 535323ec0c0..6f7bd435560 100644
--- a/Common/include/geometry/elements/CElement.hpp
+++ b/Common/include/geometry/elements/CElement.hpp
@@ -67,6 +67,8 @@ class CElement {
   su2activematrix NodalExtrap;  /*!< \brief Coordinates of the nodal points for Gaussian extrapolation. */
   su2activematrix NodalStress;  /*!< \brief Stress at the nodes. */
 
+  su2activevector NodalTemperature; /*!< \brief Temperature at the nodes. */
+
   /*--- Stiffness and load matrices. ---*/
   std::vector<su2activematrix> Kab; /*!< \brief Structure for the constitutive component of the tangent matrix. */
   su2activematrix Mab;              /*!< \brief Structure for the nodal components of the mass matrix. */
@@ -151,7 +153,7 @@ class CElement {
    * \param[in] iDim - Dimension.
    * \param[in] val_CoordRef - Value of the coordinate.
    */
-  inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordRef) {
+  inline void SetRef_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordRef) {
     RefCoord(iNode, iDim) = val_CoordRef;
   }
 
@@ -161,10 +163,22 @@ class CElement {
    * \param[in] iDim - Dimension.
    * \param[in] val_CoordRef - Value of the coordinate.
    */
-  inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, su2double val_CoordCurr) {
+  inline void SetCurr_Coord(unsigned short iNode, unsigned short iDim, const su2double& val_CoordCurr) {
     CurrentCoord(iNode, iDim) = val_CoordCurr;
   }
 
+  /*!
+   * \brief Set the value of the temperature of a node.
+   */
+  inline void SetTemperature(unsigned short iNode, const su2double& val_Temperature) {
+    NodalTemperature[iNode] = val_Temperature;
+  }
+
+  /*!
+   * \brief Set the value of the temperature of all nodes.
+   */
+  inline void SetTemperature(const su2double& val_Temperature) { NodalTemperature = val_Temperature; }
+
   /*!
    * \brief Get the value of the coordinate of the nodes in the reference configuration.
    * \param[in] iNode - Index of node.
@@ -181,6 +195,17 @@ class CElement {
    */
   inline su2double GetCurr_Coord(unsigned short iNode, unsigned short iDim) const { return CurrentCoord(iNode, iDim); }
 
+  /*!
+   * \brief Get the value of the temperature at a Gaussian integration point.
+   */
+  inline su2double GetTemperature(unsigned short iGauss) const {
+    su2double Temperature = 0;
+    for (auto iNode = 0u; iNode < nNodes; ++iNode) {
+      Temperature += GetNi(iNode, iGauss) * NodalTemperature[iNode];
+    }
+    return Temperature;
+  }
+
   /*!
    * \brief Get the weight of the corresponding Gaussian Point.
    * \param[in] iGauss - index of the Gaussian point.
@@ -214,7 +239,9 @@ class CElement {
    * \param[in] nodeB - index of Node b.
    * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution.
    */
-  inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, su2double val_Mab) { Mab(nodeA, nodeB) += val_Mab; }
+  inline void Add_Mab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Mab) {
+    Mab(nodeA, nodeB) += val_Mab;
+  }
 
   /*!
    * \brief Add the value of a submatrix K relating nodes a and b, for the constitutive term.
@@ -243,7 +270,7 @@ class CElement {
    * \param[in] nodeB - index of Node b.
    * \param[in] val_Ks_ab - value of the term that will constitute the diagonal of the stress contribution.
    */
-  inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, su2double val_Ks_ab) {
+  inline void Add_Ks_ab(unsigned short nodeA, unsigned short nodeB, const su2double& val_Ks_ab) {
     Ks_ab(nodeA, nodeB) += val_Ks_ab;
   }
 
@@ -354,7 +381,7 @@ class CElement {
    * \param[in] iVar - Variable index.
    * \param[in] val_Stress - Value of the stress added.
    */
-  inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, su2double val_Stress) {
+  inline void Add_NodalStress(unsigned short iNode, unsigned short iVar, const su2double& val_Stress) {
     NodalStress(iNode, iVar) += val_Stress;
   }
 
@@ -789,7 +816,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> {
   /*!
    * \brief Shape functions (Ni) evaluated at point Xi,Eta.
    */
-  inline static void ShapeFunctions(su2double Xi, su2double Eta, su2double* Ni) {
+  inline static void ShapeFunctions(const su2double& Xi, const su2double& Eta, su2double* Ni) {
     Ni[0] = 0.25 * (1.0 - Xi) * (1.0 - Eta);
     Ni[1] = 0.25 * (1.0 + Xi) * (1.0 - Eta);
     Ni[2] = 0.25 * (1.0 + Xi) * (1.0 + Eta);
@@ -799,7 +826,7 @@ class CQUAD4 final : public CElementWithKnownSizes<4, 4, 2> {
   /*!
    * \brief Shape function Jacobian (dNi) evaluated at point Xi,Eta.
    */
-  inline static void ShapeFunctionJacobian(su2double Xi, su2double Eta, su2double dNi[][2]) {
+  inline static void ShapeFunctionJacobian(const su2double& Xi, const su2double& Eta, su2double dNi[][2]) {
     dNi[0][0] = -0.25 * (1.0 - Eta);
     dNi[0][1] = -0.25 * (1.0 - Xi);
     dNi[1][0] = 0.25 * (1.0 - Eta);
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 90af5948d25..5c220686271 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -2442,6 +2442,8 @@ void CConfig::SetConfig_Options() {
   addDoubleListOption("MATERIAL_DENSITY", nMaterialDensity, MaterialDensity);
   /* DESCRIPTION: Material thermal expansion coefficient */
   addDoubleListOption("MATERIAL_THERMAL_EXPANSION_COEFF", nMaterialThermalExpansion, MaterialThermalExpansion);
+  /* DESCRIPTION: Temperature at which there is no stress from thermal expansion */
+  addDoubleOption("MATERIAL_REFERENCE_TEMPERATURE", MaterialReferenceTemperature, 288.15);
   /* DESCRIPTION: Knowles B constant */
   addDoubleOption("KNOWLES_B", Knowles_B, 1.0);
   /* DESCRIPTION: Knowles N constant */
diff --git a/Common/src/geometry/elements/CElement.cpp b/Common/src/geometry/elements/CElement.cpp
index c42a3707562..5acf96a3ac9 100644
--- a/Common/src/geometry/elements/CElement.cpp
+++ b/Common/src/geometry/elements/CElement.cpp
@@ -47,6 +47,8 @@ CElement::CElement(unsigned short ngauss, unsigned short nnodes, unsigned short
 
   NodalStress.resize(nNodes, 6) = su2double(0.0);
 
+  NodalTemperature.resize(nNodes) = su2double(0.0);
+
   Mab.resize(nNodes, nNodes);
   Ks_ab.resize(nNodes, nNodes);
   Kab.resize(nNodes);
diff --git a/Common/src/geometry/elements/CQUAD4.cpp b/Common/src/geometry/elements/CQUAD4.cpp
index 07725a3af19..29d3a621281 100644
--- a/Common/src/geometry/elements/CQUAD4.cpp
+++ b/Common/src/geometry/elements/CQUAD4.cpp
@@ -67,7 +67,6 @@ CQUAD4::CQUAD4() : CElementWithKnownSizes<NGAUSS, NNODE, NDIM>() {
   /*--- Store the extrapolation functions (used to compute nodal stresses) ---*/
 
   su2double ExtrapCoord[4][2], sqrt3 = 1.732050807568877;
-  ;
 
   ExtrapCoord[0][0] = -sqrt3;
   ExtrapCoord[0][1] = -sqrt3;
diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
index 43d5ebd5a63..bbaa805e046 100644
--- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
+++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
@@ -59,6 +59,7 @@ class CFEAElasticity : public CNumerics {
   su2double Mu        = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
   su2double Lambda    = 0.0;              /*!< \brief Aux. variable, Lame's coeficient. */
   su2double Kappa     = 0.0;              /*!< \brief Aux. variable, Compressibility constant. */
+  su2double ThermalStressTerm = 0.0;      /*!< \brief Aux. variable, Relationship between stress and delta T. */
 
   su2double *E_i      = nullptr;          /*!< \brief Young's modulus of elasticity. */
   su2double *Nu_i     = nullptr;          /*!< \brief Poisson's ratio. */
@@ -66,6 +67,8 @@ class CFEAElasticity : public CNumerics {
   su2double *Rho_s_DL_i = nullptr;        /*!< \brief Structural density (for dead loads). */
   su2double *Alpha_i  = nullptr;          /*!< \brief Thermal expansion coefficient. */
 
+  su2double ReferenceTemperature = 0.0;   /*!< \brief Reference temperature for thermal expansion. */
+
   su2double **Ba_Mat = nullptr;           /*!< \brief Matrix B for node a - Auxiliary. */
   su2double **Bb_Mat = nullptr;           /*!< \brief Matrix B for node b - Auxiliary. */
   su2double *Ni_Vec  = nullptr;           /*!< \brief Vector of shape functions - Auxiliary. */
@@ -74,8 +77,6 @@ class CFEAElasticity : public CNumerics {
   su2double **GradNi_Ref_Mat = nullptr;   /*!< \brief Gradients of Ni - Auxiliary. */
   su2double **GradNi_Curr_Mat = nullptr;  /*!< \brief Gradients of Ni - Auxiliary. */
 
-  su2double *FAux_Dead_Load = nullptr;    /*!< \brief Auxiliar vector for the dead loads */
-
   su2double *DV_Val = nullptr;            /*!< \brief For optimization cases, value of the design variables. */
   unsigned short n_DV = 0;                /*!< \brief For optimization cases, number of design variables. */
 
@@ -232,6 +233,8 @@ class CFEAElasticity : public CNumerics {
     Mu     = E / (2.0*(1.0 + Nu));
     Lambda = Nu*E/((1.0+Nu)*(1.0-2.0*Nu));
     Kappa  = Lambda + (2/3)*Mu;
+    /*--- https://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php ---*/
+    ThermalStressTerm = -Alpha * E / (1 - (plane_stress ? 1 : 2) * Nu);
   }
 
   /*!
@@ -240,8 +243,8 @@ class CFEAElasticity : public CNumerics {
    * \param[in] jVar - Index j.
    * \return 1 if i=j, 0 otherwise.
    */
-  inline static su2double deltaij(unsigned short iVar, unsigned short jVar) {
-    return su2double(iVar==jVar);
+  inline static passivedouble deltaij(unsigned short iVar, unsigned short jVar) {
+    return static_cast<passivedouble>(iVar == jVar);
   }
 
 };
diff --git a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp
index 0e57f54e529..49ba78c09a2 100644
--- a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp
+++ b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp
@@ -141,8 +141,9 @@ class CFEANonlinearElasticity : public CFEAElasticity {
    * \brief Compute the stress tensor.
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
+   * \param[in] iGauss - Index of Gaussian integration point.
    */
-  virtual void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) = 0;
+  virtual void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) = 0;
 
   /*!
    * \brief Update an element with Maxwell's stress.
diff --git a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp b/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp
index a36eee0a67d..dbd7b4db422 100644
--- a/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp
+++ b/SU2_CFD/include/numerics/elasticity/nonlinear_models.hpp
@@ -73,7 +73,7 @@ class CFEM_NeoHookean_Comp final : public CFEANonlinearElasticity {
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
    */
-  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override;
+  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
 };
 
@@ -124,7 +124,7 @@ class CFEM_Knowles_NearInc final : public CFEANonlinearElasticity {
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
    */
-  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override;
+  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
 };
 
@@ -172,7 +172,7 @@ class CFEM_DielectricElastomer final : public CFEANonlinearElasticity {
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
    */
-  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override;
+  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
 };
 
@@ -222,6 +222,6 @@ class CFEM_IdealDE final : public CFEANonlinearElasticity {
    * \param[in,out] element_container - The finite element.
    * \param[in] config - Definition of the problem.
    */
-  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config) override;
+  void Compute_Stress_Tensor(CElement *element_container, const CConfig *config, unsigned short iGauss) override;
 
 };
diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
index b91be37223e..73df569ff26 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
@@ -35,7 +35,6 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   nDim = val_nDim;
   nVar = val_nVar;
 
-  bool body_forces = config->GetDeadLoad();  // Body forces (dead loads).
   bool pseudo_static = config->GetPseudoStatic();
 
   unsigned short iVar;
@@ -62,14 +61,11 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   Rho_s = Rho_s_i[0];
   Rho_s_DL = Rho_s_DL_i[0];
   Alpha = Alpha_i[0];
+  ReferenceTemperature = config->GetMaterialReferenceTemperature();
 
   Compute_Lame_Parameters();
 
-  // Auxiliary vector for body forces (dead load)
-  FAux_Dead_Load = nullptr;
-  if (body_forces) FAux_Dead_Load = new su2double [nDim];
-
-  plane_stress = (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRESS);
+  plane_stress = (nDim == 2) && (config->GetElas2D_Formulation() == STRUCT_2DFORM::PLANE_STRESS);
 
   KAux_ab = new su2double* [nDim];
   for (iVar = 0; iVar < nDim; iVar++) {
@@ -162,12 +158,11 @@ CFEAElasticity::~CFEAElasticity() {
 
   delete[] DV_Val;
 
-  delete [] FAux_Dead_Load;
-
   delete [] E_i;
   delete [] Nu_i;
   delete [] Rho_s_i;
   delete [] Rho_s_DL_i;
+  delete [] Alpha_i;
   delete [] Ni_Vec;
 }
 
@@ -242,43 +237,32 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config)
   AD::SetPreaccIn(Rho_s_DL);
   element->SetPreaccIn_Coords(false);
 
-  unsigned short iGauss, nGauss;
-  unsigned short iNode, iDim, nNode;
-
-  su2double Weight, Jac_X;
-
   /* -- Gravity directionality:
    * -- For 2D problems, we assume the direction for gravity is -y
    * -- For 3D problems, we assume the direction for gravity is -z
    */
   su2double g_force[3] = {0.0,0.0,0.0};
-
-  if (nDim == 2) g_force[1] = -1*STANDARD_GRAVITY;
-  else if (nDim == 3) g_force[2] = -1*STANDARD_GRAVITY;
+  g_force[nDim - 1] = -STANDARD_GRAVITY;
 
   element->ClearElement();       /*--- Restart the element to avoid adding over previous results. --*/
   element->ComputeGrad_Linear(); /*--- Need to compute the gradients to obtain the Jacobian. ---*/
 
-  nNode = element->GetnNodes();
-  nGauss = element->GetnGaussPoints();
-
-  for (iGauss = 0; iGauss < nGauss; iGauss++) {
+  const auto nNode = element->GetnNodes();
+  const auto nGauss = element->GetnGaussPoints();
 
-    Weight = element->GetWeight(iGauss);
-    Jac_X = element->GetJ_X(iGauss);      /*--- The dead load is computed in the reference configuration ---*/
+  for (auto iGauss = 0u; iGauss < nGauss; iGauss++) {
 
-    /*--- Retrieve the values of the shape functions for each node ---*/
-    /*--- This avoids repeated operations ---*/
-    for (iNode = 0; iNode < nNode; iNode++) {
-      Ni_Vec[iNode] = element->GetNi(iNode,iGauss);
-    }
+    const auto Weight = element->GetWeight(iGauss);
+    /*--- The dead load is computed in the reference configuration ---*/
+    const auto Jac_X = element->GetJ_X(iGauss);
 
-    for (iNode = 0; iNode < nNode; iNode++) {
+    for (auto iNode = 0; iNode < nNode; iNode++) {
+      const auto Ni = element->GetNi(iNode,iGauss);
 
-      for (iDim = 0; iDim < nDim; iDim++) {
-        FAux_Dead_Load[iDim] = Weight * Ni_Vec[iNode] * Jac_X * Rho_s_DL * g_force[iDim];
+      su2double FAux_Dead_Load[3] = {0.0,0.0,0.0};
+      for (auto iDim = 0u; iDim < nDim; iDim++) {
+        FAux_Dead_Load[iDim] = Weight * Ni * Jac_X * Rho_s_DL * g_force[iDim];
       }
-
       element->Add_FDL_a(iNode, FAux_Dead_Load);
 
     }
diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
index de7ad86219f..8e15ed0cf6b 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp
@@ -97,10 +97,11 @@ void CFEALinearElasticity::Compute_Tangent_Matrix(CElement *element, const CConf
       }
     }
 
+    const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature);
+
     for (iNode = 0; iNode < nNode; iNode++) {
 
       su2double KAux_t_a[3] = {0.0};
-      const su2double thermalStress = -E * 2e-4 / (1 - 2 * Nu);
       for (iVar = 0; iVar < nDim; iVar++) {
         KAux_t_a[iVar] += Weight * thermalStress * GradNi_Ref_Mat[iNode][iVar] * Jac_X;
       }
@@ -321,14 +322,17 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element,
 
     }
 
-    /*--- Compute the Stress Vector as D*epsilon ---*/
+    /*--- Compute the Stress Vector as D*epsilon + thermal stress ---*/
 
     su2double Stress[DIM_STRAIN_3D] = {0.0};
 
+    const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature);
+
     for (iVar = 0; iVar < bDim; iVar++) {
       for (jVar = 0; jVar < bDim; jVar++) {
         Stress[iVar] += D_Mat[iVar][jVar]*Strain[jVar];
       }
+      if (iVar < nDim) Stress[iVar] += thermalStress;
       avgStress[iVar] += Stress[iVar] / nGauss;
     }
 
@@ -370,10 +374,10 @@ CFEAMeshElasticity::CFEAMeshElasticity(unsigned short val_nDim, unsigned short v
                                        unsigned long val_nElem, const CConfig *config) :
                                        CFEALinearElasticity() {
   DV_Val         = nullptr;
-  FAux_Dead_Load = nullptr;
   Rho_s_i        = nullptr;
   Rho_s_DL_i     = nullptr;
   Nu_i           = nullptr;
+  Alpha_i        = nullptr;
 
   nDim = val_nDim;
   nVar = val_nVar;
diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp
index 0b8ed7cb855..8f74c164dd6 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp
@@ -349,7 +349,7 @@ void CFEANonlinearElasticity::Compute_Tangent_Matrix(CElement *element, const CC
 
     /*--- Compute the constitutive matrix ---*/
 
-    Compute_Stress_Tensor(element, config);
+    Compute_Stress_Tensor(element, config, iGauss);
 //    if (maxwell_stress) Add_MaxwellStress(element, config);
     Compute_Constitutive_Matrix(element, config);
 
@@ -571,7 +571,7 @@ void CFEANonlinearElasticity::Compute_NodalStress_Term(CElement *element, const
 
     /*--- Compute the stress tensor ---*/
 
-    Compute_Stress_Tensor(element, config);
+    Compute_Stress_Tensor(element, config, iGauss);
 //    if (maxwell_stress) Add_MaxwellStress(element, config);
 
     for (iNode = 0; iNode < nNode; iNode++) {
@@ -850,7 +850,7 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen
 
     /*--- Compute the stress tensor ---*/
 
-    Compute_Stress_Tensor(element, config);
+    Compute_Stress_Tensor(element, config, iGauss);
     if (maxwell_stress) Add_MaxwellStress(element, config);
 
     avgStress[0] += Stress_Tensor[0][0] / nGauss;
diff --git a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
index 55d48a29010..4b787db3416 100644
--- a/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
+++ b/SU2_CFD/src/numerics/elasticity/nonlinear_models.cpp
@@ -98,7 +98,7 @@ void CFEM_NeoHookean_Comp::Compute_Constitutive_Matrix(CElement *element, const
 
 }
 
-void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfig *config) {
+void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) {
 
   unsigned short iVar,jVar;
   su2double Mu_J = 0.0, Lambda_J = 0.0;
@@ -109,10 +109,12 @@ void CFEM_NeoHookean_Comp::Compute_Stress_Tensor(CElement *element, const CConfi
     Lambda_J = Lambda/J_F;
   }
 
+  const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature);
+
   for (iVar = 0; iVar < 3; iVar++) {
     for (jVar = 0; jVar < 3; jVar++) {
       su2double dij = deltaij(iVar,jVar);
-      Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) - E * 2e-4 / (1 - 2 * Nu)) * dij;
+      Stress_Tensor[iVar][jVar] = Mu_J * (b_Mat[iVar][jVar] - dij) + (Lambda_J * log(J_F) + thermalStress) * dij;
     }
   }
 
@@ -182,7 +184,7 @@ void CFEM_Knowles_NearInc::Compute_Constitutive_Matrix(CElement *element, const
 
 }
 
-void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfig *config) {
+void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) {
 
   /* -- Suchocki (2011) (full reference in class constructor). ---*/
 
@@ -199,10 +201,12 @@ void CFEM_Knowles_NearInc::Compute_Stress_Tensor(CElement *element, const CConfi
   Ek = Kappa * (2.0 * J_F - 1.0);
   Pr = Kappa * (J_F - 1.0);
 
+  const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature);
+
   for (iVar = 0; iVar < 3; iVar++){
     for (jVar = 0; jVar < 3; jVar++){
       Stress_Tensor[iVar][jVar] = term1 * (b_Mat_Iso[iVar][jVar] - deltaij(iVar,jVar)*trbbar) +
-                                  deltaij(iVar,jVar) * Pr;
+                                  deltaij(iVar,jVar) * (Pr + thermalStress);
     }
   }
 
@@ -234,7 +238,7 @@ void CFEM_DielectricElastomer::Compute_Constitutive_Matrix(CElement *element, co
 
 }
 
-void CFEM_DielectricElastomer::Compute_Stress_Tensor(CElement *element, const CConfig *config) {
+void CFEM_DielectricElastomer::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) {
 
   unsigned short iDim, jDim;
 
@@ -315,12 +319,11 @@ void CFEM_IdealDE::Compute_Constitutive_Matrix(CElement *element, const CConfig
 
 }
 
-void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *config) {
+void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *config, unsigned short iGauss) {
 
   /* -- Zhao, X. and Suo, Z. (2008) (full reference in class constructor). ---*/
 
   unsigned short iVar, jVar;
-  su2double dij = 0.0;
 
   /*--- Compute the isochoric deformation gradient Fbar and left Cauchy-Green tensor bbar ---*/
   Compute_Isochoric_F_b();
@@ -333,13 +336,12 @@ void CFEM_IdealDE::Compute_Stress_Tensor(CElement *element, const CConfig *confi
   Pr = Kappa * (J_F - 1.0);
   Eg23 = 2.0 * Eg / 3.0;
 
-  // Stress tensor
+  const su2double thermalStress = ThermalStressTerm * (element->GetTemperature(iGauss) - ReferenceTemperature);
 
   for (iVar = 0; iVar < 3; iVar++){
-    for (jVar = 0; jVar < 3; jVar++){
-      if (iVar == jVar) dij = 1.0;
-      else if (iVar != jVar) dij = 0.0;
-      Stress_Tensor[iVar][jVar] = Eg * ( b_Mat_Iso[iVar][jVar] - dij * trbbar) + dij * Pr ;
+    for (jVar = 0; jVar < 3; jVar++) {
+      const su2double dij = deltaij(iVar, jVar);
+      Stress_Tensor[iVar][jVar] = Eg * ( b_Mat_Iso[iVar][jVar] - dij * trbbar) + dij * (Pr + thermalStress);
     }
   }
 
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index 4cd718d6ce1..44a19e1421a 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -71,6 +71,10 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge
       element_container[FEA_TERM][EL_TRIA+offset] = new CTRIA1();
       element_container[FEA_TERM][EL_QUAD+offset] = new CQUAD4();
 
+      /*--- Initialize temperature ---*/
+      element_container[FEA_TERM][EL_TRIA+offset]->SetTemperature(config->GetTemperature_FreeStream());
+      element_container[FEA_TERM][EL_QUAD+offset]->SetTemperature(config->GetTemperature_FreeStream());
+
       if (de_effects) {
         element_container[DE_TERM][EL_TRIA+offset] = new CTRIA1();
         element_container[DE_TERM][EL_QUAD+offset] = new CQUAD4();
@@ -82,6 +86,10 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge
       element_container[FEA_TERM][EL_PYRAM+offset] = new CPYRAM5();
       element_container[FEA_TERM][EL_PRISM+offset] = new CPRISM6();
 
+      for (const auto el : {EL_TETRA, EL_HEXA, EL_PYRAM, EL_PRISM}) {
+        element_container[FEA_TERM][el+offset]->SetTemperature(config->GetTemperature_FreeStream());
+      }
+
       if (de_effects) {
         element_container[DE_TERM][EL_TETRA+offset] = new CTETRA1();
         element_container[DE_TERM][EL_HEXA +offset] = new CHEXA8 ();
diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
index 128958ecb19..b08a38181d6 100644
--- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
+++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
@@ -14,13 +14,15 @@ MATERIAL_MODEL= LINEAR_ELASTIC
 MESH_FILENAME= meshBeam_3d.su2
 ELASTICITY_MODULUS=3E7
 POISSON_RATIO=0.3
+MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5
+MATERIAL_REFERENCE_TEMPERATURE= 288.15
+FREESTREAM_TEMPERATURE= 350
 MATERIAL_DENSITY=7854
 MARKER_CLAMPED = ( left )
 MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0)
-%MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0)
 LINEAR_SOLVER= CONJUGATE_GRADIENT
 LINEAR_SOLVER_PREC= ILU
-LINEAR_SOLVER_ERROR= 1E-9
+LINEAR_SOLVER_ERROR= 1E-8
 LINEAR_SOLVER_ITER= 1000
 MESH_FORMAT= SU2
 TABULAR_FORMAT= CSV
@@ -29,4 +31,4 @@ VOLUME_FILENAME= beam
 RESTART_FILENAME= restart_beam.dat
 SOLUTION_FILENAME= restart_beam.dat
 OUTPUT_WRT_FREQ= 1
-INNER_ITER=10
+INNER_ITER=1

From bac0e067fe3169496b1484346be42d0d2177497a Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Sun, 22 Dec 2024 18:15:30 -0800
Subject: [PATCH 17/21] update testcase

---
 TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg | 5 +++--
 config_template.cfg                             | 6 ++++++
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
index b08a38181d6..56b57777a55 100644
--- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
+++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg
@@ -18,8 +18,9 @@ MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5
 MATERIAL_REFERENCE_TEMPERATURE= 288.15
 FREESTREAM_TEMPERATURE= 350
 MATERIAL_DENSITY=7854
-MARKER_CLAMPED = ( left )
-MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0, upper, 0, right, 0)
+MARKER_CLAMPED = ( left, right )
+MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0 )
+MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 )
 LINEAR_SOLVER= CONJUGATE_GRADIENT
 LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ERROR= 1E-8
diff --git a/config_template.cfg b/config_template.cfg
index 8a953b7be78..2ecb8ca03b2 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1727,6 +1727,12 @@ ELASTICITY_MODULUS= 1000.0
 % Poisson ratio
 POISSON_RATIO= 0.35
 %
+% Thermal expansion coefficient
+MATERIAL_THERMAL_EXPANSION_COEFF= 0
+%
+% Freestream temperature at which there is no stress from thermal expansion
+MATERIAL_REFERENCE_TEMPERATURE= 288.15
+%
 % Knowles B constant
 KNOWLES_B= 1.0
 %

From f0808bd7dc0a7109b6f823ce1b854b0c18d73bad Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Mon, 23 Dec 2024 18:28:26 -0800
Subject: [PATCH 18/21] update cases

---
 Common/src/CConfig.cpp           | 5 ++++-
 TestCases/hybrid_regression.py   | 4 ++--
 TestCases/parallel_regression.py | 4 ++--
 TestCases/serial_regression.py   | 4 ++--
 4 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 5c220686271..ba2f975c699 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -914,7 +914,10 @@ void CConfig::SetPointersNull() {
   Inlet_Velocity              = nullptr;   Inflow_Mach             = nullptr;   Inflow_Pressure       = nullptr;
   Outlet_Pressure             = nullptr;   Isothermal_Temperature  = nullptr;
 
-  ElasticityMod             = nullptr;     PoissonRatio                = nullptr;     MaterialDensity       = nullptr;
+  ElasticityMod = nullptr;
+  PoissonRatio = nullptr;
+  MaterialDensity = nullptr;
+  MaterialThermalExpansion = nullptr;
 
   Load_Dir = nullptr;            Load_Dir_Value = nullptr;          Load_Dir_Multiplier = nullptr;
   Disp_Dir = nullptr;            Disp_Dir_Value = nullptr;          Disp_Dir_Multiplier = nullptr;
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index e5b0a847e61..fb05c2513db 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -680,8 +680,8 @@ def main():
     statbeam3d.cfg_dir   = "fea_fsi/StatBeam_3d"
     statbeam3d.cfg_file  = "configBeam_3d.cfg"
     statbeam3d.test_iter = 0
-    statbeam3d.test_vals         = [-2.378370, -1.585252, -2.028505, 64359.000000]
-    statbeam3d.test_vals_aarch64 = [-2.382650, -1.561882, -2.045083, 64433.000000]
+    statbeam3d.test_vals         = [-2.787802, -1.721974, -2.438436, 110350]
+    statbeam3d.test_vals_aarch64 = [-2.787802, -1.721974, -2.438436, 110350]
     test_list.append(statbeam3d)
 
     # Dynamic beam, 2d
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index f19b323f761..4a86a6fcc00 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1220,8 +1220,8 @@ def main():
     statbeam3d.cfg_dir   = "fea_fsi/StatBeam_3d"
     statbeam3d.cfg_file  = "configBeam_3d.cfg"
     statbeam3d.test_iter = 0
-    statbeam3d.test_vals         = [-8.396797, -8.162206, -8.156102, 64095.000000]
-    statbeam3d.test_vals_aarch64 = [-8.396793, -8.162255, -8.156118, 64095.0] #last 4 columns
+    statbeam3d.test_vals         = [-6.058758, -5.750933, -5.892188, 110190]
+    statbeam3d.test_vals_aarch64 = [-6.058758, -5.750933, -5.892188, 110190]
     statbeam3d.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
     test_list.append(statbeam3d)
 
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index 4aa5f883481..a0d531353f6 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -1012,8 +1012,8 @@ def main():
     statbeam3d.cfg_dir   = "fea_fsi/StatBeam_3d"
     statbeam3d.cfg_file  = "configBeam_3d.cfg"
     statbeam3d.test_iter = 0
-    statbeam3d.test_vals         = [-8.498245, -8.230816, -8.123810, 64095.000000]
-    statbeam3d.test_vals_aarch64 = [-8.498254, -8.230683, -8.123819, 64095.0] #last 4 columns
+    statbeam3d.test_vals         = [-6.168640, -5.939035, -6.071159, 110190]
+    statbeam3d.test_vals_aarch64 = [-6.168640, -5.939035, -6.071159, 110190] #last 4 columns
     test_list.append(statbeam3d)
 
     # Mix elem, 3d beam, Knowles

From 2cda5ec8ff45ea99466ceb59187f031cba654b41 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Tue, 24 Dec 2024 14:10:28 -0800
Subject: [PATCH 19/21] generalize body force for FEA solver and add
 centrifugal loads

---
 Common/include/CConfig.hpp                    |  7 ++-
 Common/src/CConfig.cpp                        | 10 ++--
 SU2_CFD/include/numerics/CNumerics.hpp        |  4 +-
 .../numerics/elasticity/CFEAElasticity.hpp    |  6 +--
 SU2_CFD/include/solvers/CFEASolver.hpp        |  9 ++--
 SU2_CFD/include/solvers/CSolver.hpp           |  6 +--
 .../numerics/elasticity/CFEAElasticity.cpp    | 51 ++++++++++++++-----
 SU2_CFD/src/output/CAdjElasticityOutput.cpp   |  4 +-
 SU2_CFD/src/solvers/CFEASolver.cpp            | 26 +++-------
 SU2_CFD/src/variables/CFEAVariable.cpp        |  2 +-
 TestCases/disc_adj_fea/configAD_fem.cfg       |  3 +-
 .../flow_load_sens/configAD_fem.cfg           |  1 -
 config_template.cfg                           |  5 +-
 13 files changed, 73 insertions(+), 61 deletions(-)

diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index d91a6011611..b5407e71eeb 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -1004,7 +1004,7 @@ class CConfig {
   bool ExtraOutput;           /*!< \brief Check if extra output need. */
   bool Wall_Functions;           /*!< \brief Use wall functions with the turbulence model */
   long ExtraHeatOutputZone;      /*!< \brief Heat solver zone with extra screen output */
-  bool DeadLoad;                 /*!< \brief Application of dead loads to the FE analysis */
+  bool CentrifugalForce;         /*!< \brief Application of centrifugal forces to the FE analysis */
   bool PseudoStatic;             /*!< \brief Application of dead loads to the FE analysis */
   bool SteadyRestart;            /*!< \brief Restart from a steady state for FSI problems. */
   su2double Newmark_beta,        /*!< \brief Parameter alpha for Newmark method. */
@@ -8942,10 +8942,9 @@ class CConfig {
   su2double GetAitkenDynMinInit(void) const { return AitkenDynMinInit; }
 
   /*!
-   * \brief Decide whether to apply dead loads to the model.
-   * \return <code>TRUE</code> if the dead loads are to be applied, <code>FALSE</code> otherwise.
+   * \brief Decide whether to apply centrifugal forces to the model.
    */
-  bool GetDeadLoad(void) const { return DeadLoad; }
+  bool GetCentrifugalForce(void) const { return CentrifugalForce; }
 
   /*!
    * \brief Identifies if the mesh is matching or not (temporary, while implementing interpolation procedures).
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index ba2f975c699..757388c9cc8 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -2507,11 +2507,11 @@ void CConfig::SetConfig_Options() {
   addEnumOption("NONLINEAR_FEM_SOLUTION_METHOD", Kind_SpaceIteScheme_FEA, Space_Ite_Map_FEA, STRUCT_SPACE_ITE::NEWTON);
   /* DESCRIPTION: Formulation for bidimensional elasticity solver */
   addEnumOption("FORMULATION_ELASTICITY_2D", Kind_2DElasForm, ElasForm_2D, STRUCT_2DFORM::PLANE_STRAIN);
-  /*  DESCRIPTION: Apply dead loads
-  *  Options: NO, YES \ingroup Config */
-  addBoolOption("DEAD_LOAD", DeadLoad, false);
+  /*  DESCRIPTION: Apply centrifugal forces
+   *  Options: NO, YES \ingroup Config */
+  addBoolOption("CENTRIFUGAL_FORCE", CentrifugalForce, false);
   /*  DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis)
-  *  Options: NO, YES \ingroup Config */
+   *  Options: NO, YES \ingroup Config */
   addBoolOption("PSEUDO_STATIC", PseudoStatic, false);
   /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */
   addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25);
@@ -3070,6 +3070,8 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){
             newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN instead.\n\n");
           else if (!option_name.compare("SPECIES_USE_STRONG_BC"))
             newString.append("SPECIES_USE_STRONG_BC is deprecated. Use MARKER_SPECIES_STRONG_BC= (marker1, ...) instead.\n\n");
+          else if (!option_name.compare("DEAD_LOAD"))
+            newString.append("DEAD_LOAD is deprecated. Use GRAVITY_FORCE or BODY_FORCE instead.\n\n");
           else {
             /*--- Find the most likely candidate for the unrecognized option, based on the length
              of start and end character sequences shared by candidates and the option. ---*/
diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index 535335947ca..2e45916edb1 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -1483,10 +1483,10 @@ class CNumerics {
   inline virtual void Compute_Mass_Matrix(CElement *element_container, const CConfig* config) { }
 
   /*!
-   * \brief A virtual member to compute the residual component due to dead loads
+   * \brief A virtual member to compute the residual component due to body forces.
    * \param[in] element_container - Element structure for the particular element integrated.
    */
-  inline virtual void Compute_Dead_Load(CElement *element_container, const CConfig* config) { }
+  inline virtual void Compute_Body_Forces(CElement *element_container, const CConfig* config) { }
 
   /*!
    * \brief A virtual member to compute the averaged nodal stresses
diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
index bbaa805e046..a1395fdf130 100644
--- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
+++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp
@@ -159,11 +159,11 @@ class CFEAElasticity : public CNumerics {
   void Compute_Mass_Matrix(CElement *element_container, const CConfig *config) final;
 
   /*!
-   * \brief Compute the nodal gravity loads for an element.
-   * \param[in,out] element_container - The element for which the dead loads are computed.
+   * \brief Compute the nodal inertial loads for an element.
+   * \param[in,out] element_container - The element for which the inertial loads are computed.
    * \param[in] config - Definition of the problem.
    */
-  void Compute_Dead_Load(CElement *element_container, const CConfig *config) final;
+  void Compute_Body_Forces(CElement *element_container, const CConfig *config) final;
 
   /*!
    * \brief Build the tangent stiffness matrix of an element.
diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp
index 4d8bf74be8e..9442fa190d7 100644
--- a/SU2_CFD/include/solvers/CFEASolver.hpp
+++ b/SU2_CFD/include/solvers/CFEASolver.hpp
@@ -98,6 +98,7 @@ class CFEASolver : public CFEASolverBase {
   bool element_based;          /*!< \brief Bool to determine if an element-based file is used. */
   bool topol_filter_applied;   /*!< \brief True if density filtering has been performed. */
   bool initial_calc = true;    /*!< \brief Becomes false after first call to Preprocessing. */
+  bool body_forces = false;    /*!< \brief Whether any body force is active. */
 
   /*!
    * \brief The highest level in the variable hierarchy this solver can safely use,
@@ -335,14 +336,14 @@ class CFEASolver : public CFEASolverBase {
                            const CConfig *config);
 
   /*!
-   * \brief Compute the dead loads.
+   * \brief Compute the inertial loads.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  void Compute_DeadLoad(CGeometry *geometry,
-                        CNumerics **numerics,
-                        const CConfig *config) final;
+  void Compute_BodyForces(CGeometry *geometry,
+                            CNumerics **numerics,
+                            const CConfig *config) final;
 
   /*!
    * \brief Clamped boundary conditions.
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index b6c0d45d219..2ed78a78a9b 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/SU2_CFD/include/solvers/CSolver.hpp
@@ -3645,9 +3645,9 @@ class CSolver {
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void Compute_DeadLoad(CGeometry *geometry,
-                                       CNumerics **numerics,
-                                       const CConfig *config) { }
+  inline virtual void Compute_BodyForces(CGeometry *geometry,
+                                           CNumerics **numerics,
+                                           const CConfig *config) { }
 
   /*!
    * \brief A virtual member. Set the volumetric heat source
diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
index 73df569ff26..ca3da2b26ae 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
@@ -27,6 +27,7 @@
 
 #include "../../../include/numerics/elasticity/CFEAElasticity.hpp"
 #include "../../../../Common/include/parallelization/omp_structure.hpp"
+#include "../../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
 
 CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
@@ -226,7 +227,11 @@ void CFEAElasticity::Compute_Mass_Matrix(CElement *element, const CConfig *confi
 }
 
 
-void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) {
+void CFEAElasticity::Compute_Body_Forces(CElement *element, const CConfig *config) {
+
+  const bool gravity = config->GetGravityForce();
+  const bool body_force = config->GetBody_Force();
+  const bool centrifugal = config->GetCentrifugalForce();
 
   /*--- Initialize values for the material model considered ---*/
   SetElement_Properties(element, config);
@@ -237,13 +242,24 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config)
   AD::SetPreaccIn(Rho_s_DL);
   element->SetPreaccIn_Coords(false);
 
-  /* -- Gravity directionality:
-   * -- For 2D problems, we assume the direction for gravity is -y
-   * -- For 3D problems, we assume the direction for gravity is -z
-   */
-  su2double g_force[3] = {0.0,0.0,0.0};
-  g_force[nDim - 1] = -STANDARD_GRAVITY;
+  std::array<su2double, 3> acceleration{};
+  if (gravity) {
+    /*--- For 2D problems, we assume gravity is in the -y direction,
+    * and for 3D problems in the -z direction. ---*/
+    acceleration[nDim - 1] = -STANDARD_GRAVITY;
+  } else if (body_force) {
+    for (auto iDim = 0u; iDim < nDim; iDim++) {
+      acceleration[iDim] = config->GetBody_Force_Vector()[iDim];
+    }
+  }
 
+  std::array<su2double, 3> center{}, omega{};
+  if (centrifugal) {
+    for (auto iDim = 0u; iDim < nDim; iDim++) {
+      center[iDim] = config->GetMotion_Origin(iDim);
+      omega[iDim] = config->GetRotation_Rate(iDim);
+    }
+  }
   element->ClearElement();       /*--- Restart the element to avoid adding over previous results. --*/
   element->ComputeGrad_Linear(); /*--- Need to compute the gradients to obtain the Jacobian. ---*/
 
@@ -256,15 +272,26 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config)
     /*--- The dead load is computed in the reference configuration ---*/
     const auto Jac_X = element->GetJ_X(iGauss);
 
-    for (auto iNode = 0; iNode < nNode; iNode++) {
+    for (auto iNode = 0u; iNode < nNode; iNode++) {
       const auto Ni = element->GetNi(iNode,iGauss);
 
-      su2double FAux_Dead_Load[3] = {0.0,0.0,0.0};
+      auto total_accel = acceleration;
+      if (centrifugal) {
+        /*--- For nonlinear this should probably use the current (deformed)
+         * coordinates, but it should not make a big difference. ---*/
+        std::array<su2double, 3> r{}, wr{}, w2r{};
+        for (auto iDim = 0u; iDim < nDim; iDim++) {
+          r[iDim] = element->GetRef_Coord(iNode, iDim) - center[iDim];
+        }
+        GeometryToolbox::CrossProduct(omega.data(), r.data(), wr.data());
+        GeometryToolbox::CrossProduct(omega.data(), wr.data(), w2r.data());
+        for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] += w2r[iDim];
+      }
+      std::array<su2double, 3> force{};
       for (auto iDim = 0u; iDim < nDim; iDim++) {
-        FAux_Dead_Load[iDim] = Weight * Ni * Jac_X * Rho_s_DL * g_force[iDim];
+        force[iDim] = Weight * Ni * Jac_X * Rho_s_DL * total_accel[iDim];
       }
-      element->Add_FDL_a(iNode, FAux_Dead_Load);
-
+      element->Add_FDL_a(iNode, force.data());
     }
 
   }
diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
index 7466b506d60..add084006ea 100644
--- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp
+++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
@@ -113,7 +113,7 @@ void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){
     if (config->GetTime_Domain() && !config->GetPseudoStatic()) {
       AddHistoryOutput("SENS_RHO_" + iVarS, "Sens[Rho" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Material density");
     }
-    if (config->GetDeadLoad()) {
+    if (config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce()) {
       AddHistoryOutput("SENS_RHO_DL_" + iVarS, "Sens[RhoDL" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Dead load density");
     }
   }
@@ -151,7 +151,7 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge
     if (config->GetTime_Domain() && !config->GetPseudoStatic()) {
       SetHistoryOutputValue("SENS_RHO_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar));
     }
-    if (config->GetDeadLoad()) {
+    if (config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce()) {
       SetHistoryOutputValue("SENS_RHO_DL_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho_DL(iVar));
     }
   }
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index 44a19e1421a..23bad36a8a4 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -49,7 +49,7 @@ CFEASolver::CFEASolver(LINEAR_SOLVER_MODE mesh_deform_mode) : CFEASolverBase(mes
 
 CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(geometry, config) {
 
-  bool dynamic = (config->GetTime_Domain());
+  bool dynamic = config->GetTime_Domain();
   config->SetDelta_UnstTimeND(config->GetDelta_UnstTime());
 
   /*--- Test whether we consider dielectric elastomers ---*/
@@ -59,6 +59,7 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge
   element_based = false;
   topol_filter_applied = false;
   initial_calc = true;
+  body_forces = config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce();
 
   /*--- Here is where we assign the kind of each element ---*/
 
@@ -121,22 +122,15 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(ge
 
   /*--- The length of the solution vector depends on whether the problem is static or dynamic ---*/
 
-  unsigned short nSolVar;
   string text_line, filename;
   ifstream restart_file;
 
-  if (dynamic) nSolVar = 3 * nVar;
-  else nSolVar = nVar;
-
-  auto* SolInit = new su2double[nSolVar]();
-
   /*--- Initialize from zero everywhere ---*/
 
-  nodes = new CFEABoundVariable(SolInit, nPoint, nDim, nVar, config);
+  std::array<su2double, 3 * MAXNVAR> zeros{};
+  nodes = new CFEABoundVariable(zeros.data(), nPoint, nDim, nVar, config);
   SetBaseClassPointerToNodes();
 
-  delete [] SolInit;
-
   /*--- Set which points are vertices and allocate boundary data. ---*/
 
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++)
@@ -567,7 +561,6 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container,
 
   const bool dynamic = config->GetTime_Domain();
   const bool disc_adj_fem = (config->GetKind_Solver() == MAIN_SOLVER::DISC_ADJ_FEM);
-  const bool body_forces = config->GetDeadLoad();
   const bool topology_mode = config->GetTopology_Optimization();
 
   /*
@@ -602,7 +595,7 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container,
    * Only initialized once, at the first iteration of the first time step.
    */
   if (body_forces && (initial_calc || disc_adj_fem))
-    Compute_DeadLoad(geometry, numerics, config);
+    Compute_BodyForces(geometry, numerics, config);
 
   /*--- Clear the linear system solution. ---*/
   SU2_OMP_PARALLEL
@@ -1415,7 +1408,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
 
 }
 
-void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
+void CFEASolver::Compute_BodyForces(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   /*--- Start OpenMP parallel region. ---*/
 
@@ -1465,7 +1458,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con
         /*--- Set the properties of the element and compute its mass matrix. ---*/
         element->Set_ElProperties(element_properties[iElem]);
 
-        numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Dead_Load(element, config);
+        numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Body_Forces(element, config);
 
         /*--- Add contributions of this element to the mass matrix. ---*/
         for (iNode = 0; iNode < nNodes; iNode++) {
@@ -2113,7 +2106,6 @@ void CFEASolver::ImplicitNewmark_Iteration(const CGeometry *geometry, CNumerics
   const bool linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL);
   const bool nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE);
   const bool newton_raphson = (config->GetKind_SpaceIteScheme_FEA() == STRUCT_SPACE_ITE::NEWTON);
-  const bool body_forces = config->GetDeadLoad();
 
   /*--- For simplicity, no incremental loading is handled with increment of 1. ---*/
   const su2double loadIncr = config->GetIncrementalLoad()? loadIncrement : su2double(1.0);
@@ -2301,7 +2293,6 @@ void CFEASolver::GeneralizedAlpha_Iteration(const CGeometry *geometry, CNumerics
   const bool linear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL);
   const bool nonlinear_analysis = (config->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE);
   const bool newton_raphson = (config->GetKind_SpaceIteScheme_FEA() == STRUCT_SPACE_ITE::NEWTON);
-  const bool body_forces = config->GetDeadLoad();
 
   /*--- Blend between previous and current timestep. ---*/
   const su2double alpha_f = config->Get_Int_Coeffs(2);
@@ -2925,9 +2916,6 @@ void CFEASolver::Compute_OFVolFrac(CGeometry *geometry, const CConfig *config)
 
 void CFEASolver::Compute_OFCompliance(CGeometry *geometry, const CConfig *config)
 {
-  /*--- Types of loads to consider ---*/
-  const bool body_forces = config->GetDeadLoad();
-
   /*--- If the loads are being applied incrementaly ---*/
   const bool incremental_load = config->GetIncrementalLoad();
 
diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp
index 7e486af3a21..4760706a814 100644
--- a/SU2_CFD/src/variables/CFEAVariable.cpp
+++ b/SU2_CFD/src/variables/CFEAVariable.cpp
@@ -46,7 +46,7 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig
    * still work as expected for the primal solver). ---*/
 
   const bool dynamic_analysis   = config->GetTime_Domain();
-  const bool body_forces        = config->GetDeadLoad();
+  const bool body_forces        = config->GetGravityForce() || config->GetBody_Force() || config->GetCentrifugalForce();
   const bool prestretch_fem     = config->GetPrestretch();  // Structure is prestretched
   const bool discrete_adjoint   = config->GetDiscrete_Adjoint();
   const bool refgeom            = config->GetRefGeom(); // Reference geometry needs to be stored
diff --git a/TestCases/disc_adj_fea/configAD_fem.cfg b/TestCases/disc_adj_fea/configAD_fem.cfg
index 188c7e451bf..fe1d2c6c5b7 100644
--- a/TestCases/disc_adj_fea/configAD_fem.cfg
+++ b/TestCases/disc_adj_fea/configAD_fem.cfg
@@ -16,7 +16,7 @@ RESTART_SOL= NO
 
 ITER = 10
 
-OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY 
+OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY
 REFERENCE_GEOMETRY = YES
 
 REFERENCE_GEOMETRY_FILENAME = reference_geometry.dat
@@ -39,7 +39,6 @@ MATERIAL_COMPRESSIBILITY= COMPRESSIBLE
 ELASTICITY_MODULUS=21000
 POISSON_RATIO=0.4
 MATERIAL_DENSITY=100
-DEAD_LOAD=NO
 FORMULATION_ELASTICITY_2D = PLANE_STRAIN
 
 NONLINEAR_FEM_SOLUTION_METHOD = NEWTON_RAPHSON
diff --git a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg
index 63edf6869cf..b5d0f15d944 100644
--- a/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg
+++ b/TestCases/py_wrapper/disc_adj_fea/flow_load_sens/configAD_fem.cfg
@@ -35,7 +35,6 @@ MATERIAL_COMPRESSIBILITY= COMPRESSIBLE
 ELASTICITY_MODULUS=21000
 POISSON_RATIO=0.4
 MATERIAL_DENSITY=100
-DEAD_LOAD=NO
 FORMULATION_ELASTICITY_2D = PLANE_STRAIN
 
 NONLINEAR_FEM_SOLUTION_METHOD = NEWTON_RAPHSON
diff --git a/config_template.cfg b/config_template.cfg
index 2ecb8ca03b2..17891c34373 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -2023,10 +2023,7 @@ NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON
 % Formulation for bidimensional elasticity solver
 FORMULATION_ELASTICITY_2D= PLANE_STRAIN
 %
-% Apply dead loads
-DEAD_LOAD= NO
-%
-% pseudo static analysis (no density in dynamic analysis)
+% Pseudo static analysis (no density in dynamic analysis)
 PSEUDO_STATIC= NO
 %
 % Dynamic or static structural analysis (deprecated -> use TIME_DOMAIN)

From bee055e4f89e2f00f3c32d186ceb6a904a1287ad Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Tue, 24 Dec 2024 19:12:01 -0800
Subject: [PATCH 20/21] add verification case

---
 .../numerics/elasticity/CFEAElasticity.cpp    |  2 +-
 .../fea_fsi/rotating_cylinder/config.cfg      | 34 +++++++++++++++++++
 TestCases/parallel_regression.py              | 12 +++++++
 3 files changed, 47 insertions(+), 1 deletion(-)
 create mode 100644 TestCases/fea_fsi/rotating_cylinder/config.cfg

diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
index ca3da2b26ae..efe0e19a0b2 100644
--- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
+++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp
@@ -285,7 +285,7 @@ void CFEAElasticity::Compute_Body_Forces(CElement *element, const CConfig *confi
         }
         GeometryToolbox::CrossProduct(omega.data(), r.data(), wr.data());
         GeometryToolbox::CrossProduct(omega.data(), wr.data(), w2r.data());
-        for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] += w2r[iDim];
+        for (auto iDim = 0; iDim < 3; ++iDim) total_accel[iDim] -= w2r[iDim];
       }
       std::array<su2double, 3> force{};
       for (auto iDim = 0u; iDim < nDim; iDim++) {
diff --git a/TestCases/fea_fsi/rotating_cylinder/config.cfg b/TestCases/fea_fsi/rotating_cylinder/config.cfg
new file mode 100644
index 00000000000..7a9d124da76
--- /dev/null
+++ b/TestCases/fea_fsi/rotating_cylinder/config.cfg
@@ -0,0 +1,34 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SU2 configuration file                                                       %
+% Case description: Spinning Cylinder                                          %
+% Author: Pedro Gomes                                                          %
+% Institution: SU2 Foundation                                                  %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+SOLVER= ELASTICITY
+MATH_PROBLEM= DIRECT
+GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS
+MATERIAL_MODEL= LINEAR_ELASTIC
+
+ELASTICITY_MODULUS= 2E11
+POISSON_RATIO= 0.3
+MATERIAL_DENSITY= 7854
+
+MARKER_SYM= ( x_minus, per_1, per_2 )
+MARKER_PRESSURE= ( inner, 0, outer, 0, x_plus, 0 )
+
+CENTRIFUGAL_FORCE= YES
+ROTATION_RATE= ( 1500, 0, 0 )
+MOTION_ORIGIN= ( 0, 0, 0 )
+
+LINEAR_SOLVER= CONJUGATE_GRADIENT
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-8
+LINEAR_SOLVER_ITER= 1000
+
+SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS
+OUTPUT_WRT_FREQ= 1
+INNER_ITER= 1
+
+MESH_FILENAME= cylinder.su2
+
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index d60bd28e07e..b0c18f755b3 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1225,6 +1225,18 @@ def main():
     statbeam3d.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
     test_list.append(statbeam3d)
 
+    # Rotating cylinder, 3d
+    rotating_cylinder_fea           = TestCase('rotating_cylinder_fea')
+    rotating_cylinder_fea.cfg_dir   = "fea_fsi/rotating_cylinder"
+    rotating_cylinder_fea.cfg_file  = "config.cfg"
+    rotating_cylinder_fea.test_iter = 0
+    # For a thin disk with the inner and outer radius of this geometry, from
+    # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4,
+    # the maximum stress is 165.6MPa, we get a Von Misses stress very close to that.
+    rotating_cylinder_fea.test_vals = [-1.986097, -1.023250, -1.022700, 47, -8.101266, 1.6458e+08]
+    rotating_cylinder_fea.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
+    test_list.append(rotating_cylinder_fea)
+
     # Dynamic beam, 2d
     dynbeam2d           = TestCase('dynbeam2d')
     dynbeam2d.cfg_dir   = "fea_fsi/DynBeam_2d"

From 28fb23ddef1d624dbbe485f122e8cab3123cfa67 Mon Sep 17 00:00:00 2001
From: Pedro Gomes <pcarruscag@gmail.com>
Date: Wed, 25 Dec 2024 10:01:45 -0800
Subject: [PATCH 21/21] update values and config template

---
 TestCases/parallel_regression.py |  3 +-
 config_template.cfg              | 49 +++++++++++++++-----------------
 2 files changed, 24 insertions(+), 28 deletions(-)

diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index b0c18f755b3..8bc45cb6a8b 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1233,8 +1233,7 @@ def main():
     # For a thin disk with the inner and outer radius of this geometry, from
     # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4,
     # the maximum stress is 165.6MPa, we get a Von Misses stress very close to that.
-    rotating_cylinder_fea.test_vals = [-1.986097, -1.023250, -1.022700, 47, -8.101266, 1.6458e+08]
-    rotating_cylinder_fea.command   = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
+    rotating_cylinder_fea.test_vals = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8]
     test_list.append(rotating_cylinder_fea)
 
     # Dynamic beam, 2d
diff --git a/config_template.cfg b/config_template.cfg
index 17891c34373..49cca43332b 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1742,12 +1742,17 @@ KNOWLES_N= 1.0
 % ID of the region we want to compute the sensitivities using direct differentiation
 FEA_ID_DIRECTDIFF= 0
 %
-%
 RESTART_STEADY_STATE= NO
 %
 % Time discretization
 TIME_DISCRE_FEA= NEWMARK_IMPLICIT
 %
+% Iterative method for non-linear structural analysis
+NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON
+%
+% Pseudo static analysis (no density in dynamic analysis)
+PSEUDO_STATIC= NO
+%
 % Parameter alpha for Newmark scheme (s)
 NEWMARK_BETA= 0.25
 %
@@ -1781,6 +1786,9 @@ MATERIAL_MODEL= LINEAR_ELASTIC
 % Compressibility of the material
 MATERIAL_COMPRESSIBILITY= COMPRESSIBLE
 %
+% Formulation for 2-dimensional elasticity solver
+FORMULATION_ELASTICITY_2D= PLANE_STRAIN
+%
 % -------------------- Dielectric effects  ------------------%
 %
 % Include DE effects
@@ -1794,6 +1802,20 @@ ELECTRIC_FIELD_MOD= 20e5
 %
 % Direction of the Electic Fields
 ELECTRIC_FIELD_DIR= (0.0, 1.0)
+%
+% ------------------------ Prestretch -----------------------%
+%
+% Consider a prestretch in the structural domain
+PRESTRETCH= NO
+%
+% Filename to input for prestretching membranes
+PRESTRETCH_FILENAME= prestretch_file.dat
+%
+% ----------------------- Body Forces -----------------------%
+%
+% Centrifugal forces due to ROTATION_RATE around MOTION_ORIGIN.
+% GRAVITY_FORCE and BODY_FORCE can also be used with the FEA solver.
+CENTRIFUGAL_FORCE= NO
 
 % -------------------- Weakly Coupled Heat ------------------%
 %
@@ -2010,31 +2032,6 @@ STRESS_PENALTY_PARAM= (1.0, 10.0)
 % Preaccumulation in the AD mode.
 PREACC= YES
 
-% ---------------- PRESTRETCH FOR STRUCTURES -------------------%
-% Consider a prestretch in the structural domain
-PRESTRETCH= NO
-%
-% Filename to input for prestretching membranes
-PRESTRETCH_FILENAME= prestretch_file.dat
-%
-% Iterative method for non-linear structural analysis
-NONLINEAR_FEM_SOLUTION_METHOD= NEWTON_RAPHSON
-%
-% Formulation for bidimensional elasticity solver
-FORMULATION_ELASTICITY_2D= PLANE_STRAIN
-%
-% Pseudo static analysis (no density in dynamic analysis)
-PSEUDO_STATIC= NO
-%
-% Dynamic or static structural analysis (deprecated -> use TIME_DOMAIN)
-DYNAMIC_ANALYSIS= NO
-%
-% Time Step for dynamic analysis (s) (deprecated -> use TIME_STEP)
-DYN_TIMESTEP= 0.0
-%
-% Total Physical Time for dual time stepping simulations (s) (deprecated -> use MAX_TIME)
-DYN_TIME= 1.0
-
 % ---------------- MESH DEFORMATION PARAMETERS (NEW SOLVER) -------------------%
 %
 % Use the reformatted pseudo-elastic solver for grid deformation