Merge pull request #2008 from su2code/fea_time_params

Unify Flow and FEA unsteady options + fix unsteady FEA adjoints + add FEA python wrapper examples

Author: pcarruscag

Common/include/CConfig.hpp

@@ -150,7 +150,6 @@ class CConfig {
   su2double CL_Target;         /*!< \brief Fixed Cl mode Target Cl. */
   su2double Confinement_Param; /*!< \brief Confinement paramenter for Vorticity Confinement method. */
   TIME_MARCHING TimeMarching;        /*!< \brief Steady or unsteady (time stepping or dual time stepping) computation. */
-  unsigned short Dynamic_Analysis;   /*!< \brief Static or dynamic structural analysis. */
   su2double FixAzimuthalLine;        /*!< \brief Fix an azimuthal line due to misalignments of the nearfield. */
   su2double **DV_Value;              /*!< \brief Previous value of the design variable. */
   su2double Venkat_LimiterCoeff;     /*!< \brief Limiter coefficient */
@@ -168,9 +167,6 @@ class CConfig {
   su2double HarmonicBalance_Period;  /*!< \brief Period of oscillation to be used with harmonic balance computations. */
   su2double Delta_UnstTime,          /*!< \brief Time step for unsteady computations. */
   Delta_UnstTimeND;                  /*!< \brief Time step for unsteady computations (non dimensional). */
-  su2double Delta_DynTime,        /*!< \brief Time step for dynamic structural computations. */
-  Total_DynTime,                  /*!< \brief Total time for dynamic structural computations. */
-  Current_DynTime;                /*!< \brief Global time of the dynamic structural computations. */
   su2double Total_UnstTime,       /*!< \brief Total time for unsteady computations. */
   Total_UnstTimeND;               /*!< \brief Total time for unsteady computations (non dimensional). */
   su2double Current_UnstTime,     /*!< \brief Global time of the unsteady simulation. */
@@ -2112,17 +2108,11 @@ class CConfig {
   su2double GetElasticyMod(unsigned short id_val) const { return ElasticityMod[id_val]; }
 
   /*!
-    * \brief Decide whether to apply DE effects to the model.
-    * \return <code>TRUE</code> if the DE effects are to be applied, <code>FALSE</code> otherwise.
-    */
+   * \brief Decide whether to apply DE effects to the model.
+   * \return <code>TRUE</code> if the DE effects are to be applied, <code>FALSE</code> otherwise.
+   */
   bool GetDE_Effects(void) const { return DE_Effects; }
 
-  /*!
-    * \brief Decide whether to predict the DE effects for the next time step.
-    * \return <code>TRUE</code> if the DE effects are to be applied, <code>FALSE</code> otherwise.
-    */
-   bool GetDE_Predicted(void);
-
   /*!
    * \brief Get the number of different electric constants.
    * \return Value of the DE modulus.
@@ -8722,34 +8712,6 @@ class CConfig {
    */
   bool GetSteadyRestart(void) const { return SteadyRestart; }
 
-  /*!
-   * \brief Provides information about the time integration of the structural analysis, and change the write in the output
-   *        files information about the iteration.
-   * \return The kind of time integration: Static or dynamic analysis
-   */
-  unsigned short GetDynamic_Analysis(void) const { return Dynamic_Analysis; }
-
-  /*!
-   * \brief If we are prforming an unsteady simulation, there is only
-   *        one value of the time step for the complete simulation.
-   * \return Value of the time step in an unsteady simulation (non dimensional).
-   */
-  su2double GetDelta_DynTime(void) const { return Delta_DynTime; }
-
-  /*!
-   * \brief If we are prforming an unsteady simulation, there is only
-   *        one value of the time step for the complete simulation.
-   * \return Value of the time step in an unsteady simulation (non dimensional).
-   */
-  su2double GetTotal_DynTime(void) const { return Total_DynTime; }
-
-  /*!
-   * \brief If we are prforming an unsteady simulation, there is only
-   *        one value of the time step for the complete simulation.
-   * \return Value of the time step in an unsteady simulation (non dimensional).
-   */
-  su2double GetCurrent_DynTime(void) const { return Current_DynTime; }
-
   /*!
    * \brief Get the current instance.
    * \return Current instance identifier.
@@ -8781,22 +8743,10 @@ class CConfig {
   unsigned short GetnIntCoeffs(void) const { return nIntCoeffs; }
 
   /*!
-   * \brief Get the number of different values for the elasticity modulus.
-   * \return Number of different values for the elasticity modulus.
+   * \brief Get the number of different materials for the elasticity solver.
+   * \return Number of different materials.
    */
-  unsigned short GetnElasticityMod(void) const { return nElasticityMod; }
-
-  /*!
-   * \brief Get the number of different values for the Poisson ratio.
-   * \return Number of different values for the Poisson ratio.
-   */
-  unsigned short GetnPoissonRatio(void) const { return nPoissonRatio; }
-
-  /*!
-   * \brief Get the number of different values for the Material density.
-   * \return Number of different values for the Material density.
-   */
-  unsigned short GetnMaterialDensity(void) const { return nMaterialDensity; }
+  unsigned short GetnElasticityMat(void) const { return nElasticityMod; }
 
   /*!
    * \brief Get the integration coefficients for the Generalized Alpha - Newmark integration integration scheme.
@@ -9282,12 +9232,6 @@ class CConfig {
    */
   void SetnTime_Iter(unsigned long val_iter) { nTimeIter = val_iter; }
 
-  /*!
-   * \brief Get the number of pseudo-time iterations
-   * \return Number of pseudo-time steps run for the single-zone problem
-   */
-  unsigned long GetnIter(void) const { return nIter; }
-
   /*!
    * \brief Get the restart iteration
    * \return Iteration for the restart of multizone problems

Common/include/geometry/elements/CElementProperty.hpp

@@ -93,7 +93,7 @@ class CProperty {
   /*!
    * \brief Extract the derivative of the Design density.
    */
-  inline virtual su2double GetAdjointDensity(void) const { return 0.0; }
+  inline virtual su2double GetAdjointDensity(void) { return 0.0; }
 
   /*!
    * \brief Register the Design density as an AD input variable.
@@ -177,7 +177,11 @@ class CElementProperty final : public CProperty {
   /*!
    * \brief Extract the derivative of the Design density.
    */
-  inline su2double GetAdjointDensity(void) const override { return SU2_TYPE::GetDerivative(design_rho); }
+  inline su2double GetAdjointDensity(void) override {
+    su2double der = SU2_TYPE::GetDerivative(design_rho);
+    AD::ResetInput(design_rho);
+    return der;
+  }
 
   /*!
    * \brief Register the Design density as an AD input variable.

Common/include/option_structure.hpp

@@ -1434,14 +1434,12 @@ static const MapType<std::string, ENUM_HEAT_TIMESTEP> Heat_TimeStep_Map = {
  * \brief Type of time integration schemes
  */
 enum class STRUCT_TIME_INT {
-  CD_EXPLICIT,       /*!< \brief Support for implementing an explicit method. */
   NEWMARK_IMPLICIT,  /*!< \brief Implicit Newmark integration definition. */
   GENERALIZED_ALPHA, /*!< \brief Support for implementing another implicit method. */
 };
 static const MapType<std::string, STRUCT_TIME_INT> Time_Int_Map_FEA = {
-  MakePair("CD_EXPLICIT", STRUCT_TIME_INT::CD_EXPLICIT)
   MakePair("NEWMARK_IMPLICIT", STRUCT_TIME_INT::NEWMARK_IMPLICIT)
-  MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA)
+  // MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) Not fully implemented.
 };
 
 /*!
@@ -2368,18 +2366,6 @@ enum class RECORDING {
   SOLUTION_AND_MESH,
 };
 
-/*!
- * \brief Types of schemes for dynamic structural computations
- */
-enum ENUM_DYNAMIC {
-  STATIC = 0,     /*!< \brief A static structural computation. */
-  DYNAMIC = 1     /*!< \brief Use a time stepping strategy for dynamic computations. */
-};
-static const MapType<std::string, ENUM_DYNAMIC> Dynamic_Map = {
-  MakePair("NO", STATIC)
-  MakePair("YES", DYNAMIC)
-};
-
 /*!
  * \brief Types of input file formats
  */

Common/src/CConfig.cpp

@@ -2435,12 +2435,6 @@ void CConfig::SetConfig_Options() {
   /*  DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis)
   *  Options: NO, YES \ingroup Config */
   addBoolOption("PSEUDO_STATIC", PseudoStatic, false);
-  /* DESCRIPTION: Dynamic or static structural analysis */
-  addEnumOption("DYNAMIC_ANALYSIS", Dynamic_Analysis, Dynamic_Map, STATIC);
-  /* DESCRIPTION: Time Step for dynamic analysis (s) */
-  addDoubleOption("DYN_TIMESTEP", Delta_DynTime, 0.0);
-  /* DESCRIPTION: Total Physical Time for dual time stepping simulations (s) */
-  addDoubleOption("DYN_TIME", Total_DynTime, 1.0);
   /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */
   addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25);
   /* DESCRIPTION: Parameter delta for Newmark scheme (s) */
@@ -2987,6 +2981,12 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){
           newString.append("\n");
           if (!option_name.compare("SINGLEZONE_DRIVER"))
             newString.append("Option SINGLEZONE_DRIVER is deprecated, it does not have a replacement.\n\n");
+          else if (!option_name.compare("DYN_TIMESTEP"))
+            newString.append("DYN_TIMESTEP is deprecated. Use TIME_STEP instead.\n\n");
+          else if (!option_name.compare("DYN_TIME"))
+            newString.append("DYN_TIME is deprecated. Use MAX_TIME instead.\n\n");
+          else if (!option_name.compare("DYNAMIC_ANALYSIS"))
+            newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN 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. ---*/
@@ -3660,7 +3660,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   if (Time_Domain){
     Delta_UnstTime = Time_Step;
-    Delta_DynTime  = Time_Step;
 
     if (TimeMarching == TIME_MARCHING::TIME_STEPPING){ InnerIter = 1; }
 
@@ -4649,6 +4648,11 @@ 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 (nElectric_Constant == 0) {
     nElectric_Constant = 1;
     Electric_Constant = new su2double[1]; Electric_Constant[0] = 0.0;
@@ -6803,7 +6807,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
   else {
     if (Time_Domain) {
       cout << "Dynamic structural analysis."<< endl;
-      cout << "Time step provided by the user for the dynamic analysis(s): "<< Delta_DynTime << "." << endl;
+      cout << "Time step provided by the user for the dynamic analysis(s): "<< Time_Step << "." << endl;
     } else {
       cout << "Static structural analysis." << endl;
     }
@@ -6871,14 +6875,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
 
     if (fea) {
       switch (Kind_TimeIntScheme_FEA) {
-        case STRUCT_TIME_INT::CD_EXPLICIT:
-          cout << "Explicit time integration (NOT IMPLEMENTED YET)." << endl;
-          break;
         case STRUCT_TIME_INT::GENERALIZED_ALPHA:
           cout << "Generalized-alpha method." << endl;
           break;
         case STRUCT_TIME_INT::NEWMARK_IMPLICIT:
-          if (Dynamic_Analysis) cout << "Newmark implicit method for the structural time integration." << endl;
+          if (Time_Domain) cout << "Newmark implicit method for the structural time integration." << endl;
           switch (Kind_Linear_Solver) {
             case BCGSTAB:
               cout << "BCGSTAB is used for solving the linear system." << endl;
@@ -8416,9 +8417,6 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver,
 
     case MAIN_SOLVER::FEM_ELASTICITY:
     case MAIN_SOLVER::DISC_ADJ_FEM:
-
-      Current_DynTime = static_cast<su2double>(TimeIter)*Delta_DynTime;
-
       if (val_system == RUNTIME_FEA_SYS) {
         SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE);
         SetKind_TimeIntScheme(NONE);
@@ -9994,7 +9992,6 @@ void CConfig::SetMultizone(const CConfig *driver_config, const CConfig* const* c
   /*--- Fix the Time Step for all subdomains, for the case of time-dependent problems ---*/
   if (driver_config->GetTime_Domain()){
     Delta_UnstTime = driver_config->GetTime_Step();
-    Delta_DynTime  = driver_config->GetTime_Step();
 
     Time_Domain = true;
   }

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -499,6 +499,15 @@ class CDriverBase {
         main_geometry->GetnVertex(iMarker), "MarkerPrimitives", false);
   }
 
+  /*!
+   * \brief Get a read-only view of the geometry sensitivity of a discrete adjoint solver.
+   */
+  inline CPyWrapperMatrixView Sensitivity(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    auto& sensitivity = const_cast<su2activematrix&>(solver->GetNodes()->GetSensitivity());
+    return CPyWrapperMatrixView(sensitivity, "Sensitivity", true);
+  }
+
   /*!
    * \brief Set the temperature of a vertex on a specified marker (MARKER_PYTHON_CUSTOM).
    * \note This can be the input of a heat or flow solver in a CHT setting.

SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp

@@ -50,7 +50,8 @@ class CDiscAdjFEASolver final : public CSolver {
     su2double* val = nullptr;         /*!< \brief Value of the variable. */
     su2double* LocalSens = nullptr;   /*!< \brief Local sensitivity (domain). */
     su2double* GlobalSens = nullptr;  /*!< \brief Global sensitivity (mpi). */
-    su2double* TotalSens = nullptr;   /*!< \brief Total sensitivity (time domain). */
+    su2double* OldSens = nullptr;     /*!< \brief Previous global sensitivity, used to update the total. */
+    su2double* TotalSens = nullptr;   /*!< \brief Total sensitivity (integrated over time). */
 
     su2double& operator[] (unsigned short i) { return val[i]; }
     const su2double& operator[] (unsigned short i) const { return val[i]; }
@@ -61,6 +62,7 @@ class CDiscAdjFEASolver final : public CSolver {
       val = new su2double[n]();
       LocalSens = new su2double[n]();
       GlobalSens = new su2double[n]();
+      OldSens = new su2double[n]();
       TotalSens = new su2double[n]();
     }
 
@@ -69,6 +71,7 @@ class CDiscAdjFEASolver final : public CSolver {
       delete [] val;
       delete [] LocalSens;
       delete [] GlobalSens;
+      delete [] OldSens;
       delete [] TotalSens;
     }
 
@@ -80,10 +83,19 @@ class CDiscAdjFEASolver final : public CSolver {
       for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]);
 
       SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
+
+      for (auto i = 0u; i < size; ++i) {
+        /*--- Update the total by subtracting the old and adding the new value.
+         * Then update the old value for the next call to this function. ---*/
+        TotalSens[i] += GlobalSens[i] - OldSens[i];
+        OldSens[i] = GlobalSens[i];
+      }
     }
 
-    void UpdateTotal() {
-      for (auto i = 0u; i < size; ++i) TotalSens[i] += GlobalSens[i];
+    void Store() {
+      /*--- Clears the old values such that on the next time step the total is
+       * incremented instead of updated. ---*/
+      for (auto i = 0u; i < size; ++i) OldSens[i] = 0.0;
     }
 
     ~SensData() { clear(); }
@@ -213,42 +225,6 @@ class CDiscAdjFEASolver final : public CSolver {
    */
   inline su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.TotalSens[iDVFEA]; }
 
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the Young Modulus E
-   */
-  inline su2double GetGlobal_Sens_E(unsigned short iVal) const override { return E.GlobalSens[iVal]; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the Mach sensitivity for the Poisson's ratio Nu
-   */
-  inline su2double GetGlobal_Sens_Nu(unsigned short iVal) const override { return Nu.GlobalSens[iVal]; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the Electric Field in the region iEField
-   */
-  inline su2double GetGlobal_Sens_EField(unsigned short iEField) const override { return EField.GlobalSens[iEField]; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA
-   */
-  inline su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.GlobalSens[iDVFEA]; }
-
-  /*!
-   * \brief Get the total sensitivity for the structural density
-   * \return Value of the structural density sensitivity
-   */
-  inline su2double GetGlobal_Sens_Rho(unsigned short iVal) const override { return Rho.GlobalSens[iVal]; }
-
-  /*!
-   * \brief Get the total sensitivity for the structural weight
-   * \return Value of the structural weight sensitivity
-   */
-  inline su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const override { return Rho_DL.GlobalSens[iVal]; }
-
   /*!
    * \brief Get the value of the Young modulus from the adjoint solver
    * \return Value of the Young modulus from the adjoint solver

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -682,28 +682,9 @@ class CFEASolver : public CFEASolverBase {
   inline su2double GetFSI_ConvValue(unsigned short val_index) const final { return FSI_Conv[val_index]; }
 
   /*!
-   * \brief Retrieve the value of the dynamic Aitken relaxation factor.
-   * \return Value of the dynamic Aitken relaxation factor.
+   * \brief Store the value of the last Aitken relaxation factor in the current time step.
    */
-  inline su2double GetWAitken_Dyn(void) const final { return WAitken_Dyn; }
-
-  /*!
-   * \brief Retrieve the value of the last Aitken relaxation factor in the previous time step.
-   * \return Value of the last Aitken relaxation factor in the previous time step.
-   */
-  inline su2double GetWAitken_Dyn_tn1(void) const final { return WAitken_Dyn_tn1; }
-
-  /*!
-   * \brief Set the value of the dynamic Aitken relaxation factor
-   * \param[in] Value of the dynamic Aitken relaxation factor
-   */
-  inline void SetWAitken_Dyn(su2double waitk) final { WAitken_Dyn = waitk; }
-
-  /*!
-   * \brief Set the value of the last Aitken relaxation factor in the current time step.
-   * \param[in] Value of the last Aitken relaxation factor in the current time step.
-   */
-  inline void SetWAitken_Dyn_tn1(su2double waitk_tn1) final { WAitken_Dyn_tn1 = waitk_tn1; }
+  inline void SetWAitken_Dyn_tn1() final { WAitken_Dyn_tn1 = WAitken_Dyn; }
 
   /*!
    * \brief Set the value of the load increment for nonlinear structural analysis