Merge branch 'coupled_thermoelasticity' of github.com:Vaish-W/SU2 int…

…o coupled_thermoelasticity

Add nodal temperature integration in Compute_StiffMatrix

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

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. ---*/

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

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.

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -99,6 +99,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 Pointer to the heat solver for coupled simulations.
@@ -352,14 +353,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.

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

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

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

SU2_CFD/src/solvers/CFEASolver.cpp

@@ -50,7 +50,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 ---*/
@@ -60,6 +60,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 ---*/
 
@@ -122,22 +123,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++)
@@ -568,7 +562,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();
 
   if (config->GetWeakly_Coupled_Heat()) {
@@ -607,7 +600,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
@@ -1424,7 +1417,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. ---*/
 
@@ -1474,7 +1467,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++) {
@@ -2122,7 +2115,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);
@@ -2310,7 +2302,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);
@@ -2934,9 +2925,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();
 

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

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

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
+

TestCases/parallel_regression.py

@@ -1225,6 +1225,17 @@ 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 = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8]
+    test_list.append(rotating_cylinder_fea)
+
     # Dynamic beam, 2d
     dynbeam2d           = TestCase('dynbeam2d')
     dynbeam2d.cfg_dir   = "fea_fsi/DynBeam_2d"