cleanup the time integration of fea material sensitivities

Author: pcarruscag

Common/include/CConfig.hpp

@@ -2108,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.
@@ -8749,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.
-   */
-  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.
+   * \brief Get the number of different materials for the elasticity solver.
+   * \return Number of different materials.
    */
-  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.

Common/src/CConfig.cpp

@@ -4648,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;

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];
+      /*--- 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/CSolver.hpp

@@ -3201,42 +3201,6 @@ class CSolver {
    */
   inline virtual su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; }
 
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the Young Modulus E
-   */
-  inline virtual su2double GetGlobal_Sens_E(unsigned short iVal) const { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the Poisson's ratio Nu
-   */
-  inline virtual su2double GetGlobal_Sens_Nu(unsigned short iVal) const { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the structural density sensitivity
-   */
-  inline virtual su2double GetGlobal_Sens_Rho(unsigned short iVal) const { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the structural weight sensitivity
-   */
-  inline virtual su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the Electric Field in the region iEField
-   */
-  inline virtual su2double GetGlobal_Sens_EField(unsigned short iEField) const { return 0.0; }
-
-  /*!
-   * \brief A virtual member.
-   * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA
-   */
-  inline virtual su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; }
-
   /*!
    * \brief A virtual member.
    * \return Value of the Young modulus from the adjoint solver

SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp

@@ -32,32 +32,6 @@
 
 CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) {
   fem_iteration = new CFEAIteration(config);
-
-  // TEMPORARY output only for standalone structural problems
-  if (config->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) {
-    bool de_effects = config->GetDE_Effects();
-    unsigned short iVar;
-
-    /*--- Header of the temporary output file ---*/
-    ofstream myfile_res;
-    myfile_res.open("Results_Reverse_Adjoint.txt");
-
-    myfile_res << "Obj_Func"
-               << " ";
-    for (iVar = 0; iVar < config->GetnElasticityMod(); iVar++) myfile_res << "Sens_E_" << iVar << "\t";
-
-    for (iVar = 0; iVar < config->GetnPoissonRatio(); iVar++) myfile_res << "Sens_Nu_" << iVar << "\t";
-
-    if (config->GetTime_Domain()) {
-      for (iVar = 0; iVar < config->GetnMaterialDensity(); iVar++) myfile_res << "Sens_Rho_" << iVar << "\t";
-    }
-
-    if (de_effects) {
-      for (iVar = 0; iVar < config->GetnElectric_Field(); iVar++) myfile_res << "Sens_EField_" << iVar << "\t";
-    }
-
-    myfile_res << "\n";
-  }
 }
 
 CDiscAdjFEAIteration::~CDiscAdjFEAIteration() = default;
@@ -71,31 +45,49 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat
   auto dirNodes = solvers0[FEA_SOL]->GetNodes();
   auto adjNodes = solvers0[ADJFEA_SOL]->GetNodes();
 
-  /*--- For the dynamic adjoint, load direct solutions from restart files. ---*/
+  auto StoreDirectSolution = [&]() {
+    for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) {
+      adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint));
+    }
+  };
+
+  /*--- For the dynamic adjoint, load direct solutions from restart files.
+   * For steady, store the direct solution to be able to reset it later. ---*/
 
   if (config[iZone]->GetTime_Domain()) {
     const int TimeIter = config[iZone]->GetTimeIter();
     const int Direct_Iter = SU2_TYPE::Int(config[iZone]->GetUnst_AdjointIter()) - TimeIter - 1;
 
-    /*--- We want to load the already converged solution at timesteps n and n-1 ---*/
-
-    /*--- Load solution at timestep n-1 ---*/
+    /*--- We need the already converged solution at timesteps n and n-1. On the first
+     * adjoint time step we load both, then n-1 becomes "n" and we only load n-2. ---*/
 
-    LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
+    if (TimeIter > 0) {
+      /*--- Save n-1 to become n. ---*/
+      for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) {
+        adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution_time_n(iPoint));
+      }
+    }
 
-    /*--- Push solution back to correct array ---*/
+    /*--- Load solution at timestep n-1 and push back to correct array. ---*/
 
+    LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1);
     dirNodes->Set_Solution_time_n();
 
-    /*--- Load solution timestep n ---*/
+    /*--- Load or set solution at timestep n. ---*/
 
-    LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter);
-  }
-
-  /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later. ---*/
+    if (TimeIter == 0) {
+      LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter);
+      StoreDirectSolution();
+    } else {
+      /*--- Set n-1 as n. ---*/
+      for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++)
+        for (auto iVar = 0u; iVar < solvers0[ADJFEA_SOL]->GetnVar(); iVar++)
+          dirNodes->SetSolution(iPoint, iVar, adjNodes->GetSolution_Direct(iPoint)[iVar]);
+    }
 
-  for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) {
-    adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint));
+  } else {
+    /*--- Steady. ---*/
+    StoreDirectSolution();
   }
 
   solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false);
@@ -178,7 +170,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge
   const int mat_knowles = MAT_KNOWLES+offset;
   const int de_term = DE_TERM+offset;
 
-  for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMod(); iProp++) {
+  for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMat(); iProp++) {
     su2double E = adj_solver->GetVal_Young(iProp);
     su2double nu = adj_solver->GetVal_Poisson(iProp);
     su2double rho = adj_solver->GetVal_Rho(iProp);
@@ -302,45 +294,8 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra
                                        CSolver***** solver, CNumerics****** numerics, CConfig** config,
                                        CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement,
                                        CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) {
-  const bool dynamic = (config[iZone]->GetTime_Domain());
   auto solvers0 = solver[iZone][iInst][MESH_0];
 
-  // TEMPORARY output only for standalone structural problems
-  if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) {
-    unsigned short iVar;
-
-    const bool de_effects = config[iZone]->GetDE_Effects();
-
-    /*--- Header of the temporary output file ---*/
-    ofstream myfile_res;
-    myfile_res.open("Results_Reverse_Adjoint.txt", ios::app);
-
-    myfile_res.precision(15);
-
-    myfile_res << config[iZone]->GetTimeIter() << "\t";
-
-    solvers0[FEA_SOL]->Evaluate_ObjFunc(config[iZone], solvers0);
-    myfile_res << scientific << solvers0[FEA_SOL]->GetTotal_ComboObj() << "\t";
-
-    for (iVar = 0; iVar < config[iZone]->GetnElasticityMod(); iVar++)
-      myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_E(iVar) << "\t";
-    for (iVar = 0; iVar < config[iZone]->GetnPoissonRatio(); iVar++)
-      myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Nu(iVar) << "\t";
-    if (dynamic) {
-      for (iVar = 0; iVar < config[iZone]->GetnMaterialDensity(); iVar++)
-        myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar) << "\t";
-    }
-    if (de_effects) {
-      for (iVar = 0; iVar < config[iZone]->GetnElectric_Field(); iVar++)
-        myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_EField(iVar) << "\t";
-    }
-    for (iVar = 0; iVar < solvers0[ADJFEA_SOL]->GetnDVFEA(); iVar++) {
-      myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iVar) << "\t";
-    }
-
-    myfile_res << "\n";
-  }
-
   // TEST: for implementation of python framework in standalone structural problems
   if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) {
     /*--- Header of the temporary output file ---*/

SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp

@@ -38,23 +38,19 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar,
   bool body_forces = config->GetDeadLoad();  // Body forces (dead loads).
   bool pseudo_static = config->GetPseudoStatic();
 
-  unsigned short iVar, nProp;
+  unsigned short iVar;
 
   /*--- Initialize vector structures for multiple material definition ---*/
-  nProp = config->GetnElasticityMod();
+  const auto nProp = config->GetnElasticityMat();
 
   E_i = new su2double[nProp];
   for (iVar = 0; iVar < nProp; iVar++)
     E_i[iVar] = config->GetElasticyMod(iVar);
 
-  nProp = config->GetnPoissonRatio();
-
   Nu_i = new su2double[nProp];
   for (iVar = 0; iVar < nProp; iVar++)
     Nu_i[iVar] = config->GetPoissonRatio(iVar);
 
-  nProp = config->GetnMaterialDensity();
-
   Rho_s_i = new su2double[nProp];     // For inertial effects
   Rho_s_DL_i = new su2double[nProp];  // For dead loads