checkpoint fixing AD index issues in the flow tractions

pcarruscag · pcarruscag · commit 1549d45be67b · 2023-04-22T17:22:11.000-07:00
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
@@ -1439,7 +1439,7 @@ enum class STRUCT_TIME_INT {
 };
 static const MapType<std::string, STRUCT_TIME_INT> Time_Int_Map_FEA = {
   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.
 };
 
 /*!
diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/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.
diff --git a/SU2_CFD/include/variables/CDiscAdjVariable.hpp b/SU2_CFD/include/variables/CDiscAdjVariable.hpp
@@ -72,6 +72,7 @@ class CDiscAdjVariable : public CVariable {
   inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const final {
     return Sensitivity(iPoint,iDim);
   }
+  inline const MatrixType& GetSensitivity() const final { return Sensitivity; }
 
   /*!
    * \brief Set/store the dual time contributions to the adjoint variable.
diff --git a/SU2_CFD/include/variables/CFEABoundVariable.hpp b/SU2_CFD/include/variables/CFEABoundVariable.hpp
@@ -159,10 +159,15 @@ class CFEABoundVariable final : public CFEAVariable {
    */
   void RegisterFlowTraction() override;
 
+  /*!
+   * \brief Resets the AD indices of the flow tractions.
+   */
+  void ResetInputFlowTraction() override;
+
   /*!
    * \brief Extract the flow traction derivatives.
    */
-  inline su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const override {
+  inline su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const override {
     if (!fsi_analysis) return 0.0;
     if (!VertexMap.GetVertexIndex(iPoint)) return 0.0;
     return SU2_TYPE::GetDerivative(FlowTraction(iPoint,iDim));
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
@@ -2165,7 +2165,12 @@ class CVariable {
   /*!
    * \brief A virtual member.
    */
-  inline virtual su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
+  inline virtual void ResetInputFlowTraction() { }
+
+  /*!
+   * \brief A virtual member.
+   */
+  inline virtual su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
 
   /*!
    * \brief Register the variables in the solution array as input/output variable.
@@ -2224,6 +2229,7 @@ class CVariable {
    * \return value of the Sensitivity
    */
   inline virtual su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
+  inline virtual const MatrixType& GetSensitivity() const { AssertOverride(); return Solution; }
 
   inline virtual void SetTau_Wall(unsigned long iPoint, su2double tau_wall) {}
 
diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
@@ -193,6 +193,23 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry,
   SetVolumeOutputValue("SENSITIVITY-Y", iPoint, Node_Struc->GetSensitivity(iPoint, 1));
   if (nDim == 3)
     SetVolumeOutputValue("SENSITIVITY-Z", iPoint, Node_Struc->GetSensitivity(iPoint, 2));
+
+  if (!config->GetTime_Domain()) return;
+
+  SetVolumeOutputValue("SENS_DISP-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 0));
+  SetVolumeOutputValue("SENS_DISP-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 1));
+  if (nDim == 3)
+    SetVolumeOutputValue("SENS_DISP-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2));
+
+  SetVolumeOutputValue("SENS_VEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim));
+  SetVolumeOutputValue("SENS_VEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim + 1));
+  if (nDim == 3)
+    SetVolumeOutputValue("SENS_VEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 5));
+
+  SetVolumeOutputValue("SENS_ACCEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim));
+  SetVolumeOutputValue("SENS_ACCEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim + 1));
+  if (nDim == 3)
+    SetVolumeOutputValue("SENS_ACCEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 8));
 }
 
 void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
@@ -223,4 +240,23 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("SENSITIVITY-Z", "Sensitivity_z", "SENSITIVITY", "geometric sensitivity  in the z direction");
   /// END_GROUP
 
+  if (!config->GetTime_Domain()) return;
+
+  /*--- Sensitivities with respect to initial conditions. ---*/
+
+  AddVolumeOutput("SENS_DISP-X", "SensInitialDisp_x", "SENSITIVITY_T0", "sensitivity to the initial x displacement");
+  AddVolumeOutput("SENS_DISP-Y", "SensInitialDisp_y", "SENSITIVITY_T0", "sensitivity to the initial y displacement");
+  if (nDim == 3)
+    AddVolumeOutput("SENS_DISP-Z", "SensInitialDisp_z", "SENSITIVITY_T0", "sensitivity to the initial z displacement");
+
+  AddVolumeOutput("SENS_VEL-X", "SensInitialVel_x", "SENSITIVITY_T0", "sensitivity to the initial x velocity");
+  AddVolumeOutput("SENS_VEL-Y", "SensInitialVel_y", "SENSITIVITY_T0", "sensitivity to the initial y velocity");
+  if (nDim == 3)
+    AddVolumeOutput("SENS_VEL-Z", "SensInitialVel_z", "SENSITIVITY_T0", "sensitivity to the initial z velocity");
+
+  AddVolumeOutput("SENS_ACCEL-X", "SensInitialAccel_x", "SENSITIVITY_T0", "sensitivity to the initial x acceleration");
+  AddVolumeOutput("SENS_ACCEL-Y", "SensInitialAccel_y", "SENSITIVITY_T0", "sensitivity to the initial y acceleration");
+  if (nDim == 3)
+    AddVolumeOutput("SENS_ACCEL-Z", "SensInitialAccel_z", "SENSITIVITY_T0", "sensitivity to the initial z acceleration");
+
 }
diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp
@@ -873,8 +873,8 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry
           const auto varIdx = i % CustomOutput::MAX_VARS_PER_SOLVER;
           if (solIdx == FLOW_SOL) {
             return flowNodes->GetPrimitive(iPoint, varIdx);
-          }             return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx);
-         
+          }
+          return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx);
         } else {
           return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE];
         }
diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp
@@ -125,14 +125,18 @@ CDiscAdjFEASolver::~CDiscAdjFEASolver() { delete nodes; }
 void CDiscAdjFEASolver::SetRecording(CGeometry* geometry, CConfig *config){
 
   unsigned long iPoint;
-  unsigned short iVar;
+  unsigned short iDim, iVar;
 
   /*--- Reset the solution to the initial (converged) solution ---*/
 
-  for (iPoint = 0; iPoint < nPoint; iPoint++)
+  for (iPoint = 0; iPoint < nPoint; iPoint++) {
     for (iVar = 0; iVar < nVar; iVar++)
       direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]);
 
+    auto Coord = geometry->nodes->GetCoord(iPoint);
+    for (iDim = 0; iDim < nDim; iDim++) AD::ResetInput(Coord[iDim]);
+  }
+
   /*--- Reset the input for time n ---*/
 
   if (config->GetTime_Domain()) {
@@ -198,6 +202,8 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config,
       for (iVar = 0; iVar < nDV; iVar++) AD::ResetInput(DV[iVar]);
     }
 
+    direct_solver->GetNodes()->ResetInputFlowTraction();
+
     if (!reset) {
       E.Register();
       Nu.Register();
@@ -294,7 +300,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *c
   if (config->GetnMarker_Fluid_Load() > 0) {
     for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
       for (unsigned short iDim = 0; iDim < nDim; iDim++){
-        su2double val_sens = direct_solver->GetNodes()->ExtractFlowTraction_Sensitivity(iPoint,iDim);
+        su2double val_sens = direct_solver->GetNodes()->ExtractFlowTractionSensitivity(iPoint,iDim);
         nodes->SetFlowTractionSensitivity(iPoint, iDim, val_sens);
       }
     }
@@ -366,12 +372,12 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo
 
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
 
-    auto Coord = geometry->nodes->GetCoord(iPoint);
+    //auto Coord = geometry->nodes->GetCoord(iPoint);
 
     for (unsigned short iDim = 0; iDim < nDim; iDim++) {
 
       su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim);
-      AD::ResetInput(Coord[iDim]);
+      //AD::ResetInput(Coord[iDim]);
 
       if (!time_domain) {
         nodes->SetSensitivity(iPoint, iDim, Sensitivity);
diff --git a/SU2_CFD/src/variables/CFEABoundVariable.cpp b/SU2_CFD/src/variables/CFEABoundVariable.cpp
@@ -75,3 +75,10 @@ void CFEABoundVariable::RegisterFlowTraction() {
     for (unsigned long iVar = 0; iVar < nVar; iVar++)
       AD::RegisterInput(FlowTraction(iVertex,iVar));
 }
+
+void CFEABoundVariable::ResetInputFlowTraction() {
+  if (!fsi_analysis) return;
+  for (unsigned long iVertex = 0; iVertex < FlowTraction.rows(); iVertex++)
+    for (unsigned long iVar = 0; iVar < nVar; iVar++)
+      AD::ResetInput(FlowTraction(iVertex,iVar));
+}
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
@@ -276,7 +276,7 @@ def main():
     discadj_fsi2.cfg_dir   = "disc_adj_fsi/Airfoil_2d"
     discadj_fsi2.cfg_file  = "config.cfg"
     discadj_fsi2.test_iter = 8
-    discadj_fsi2.test_vals         = [-4.349377, 0.128475, -1.303589, 7.5407e-09, 2.3244]
+    discadj_fsi2.test_vals         = [-4.349377, 0.189968, -1.303589, 7.5407e-09, 2.3244]
     discadj_fsi2.test_vals_aarch64 = [-3.479505, 0.127953, -1.303589, 7.5407e-09, 2.3244]
     discadj_fsi2.tol       = 0.00001
     test_list.append(discadj_fsi2)
@@ -413,7 +413,7 @@ def main():
     pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea"
     pywrapper_Unst_FEA_AD.cfg_file = "config.cfg"
     pywrapper_Unst_FEA_AD.test_iter = 100
-    pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015]
+    pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.031988, 0.032015, -18.44906, -18.45085]
     pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py")
     pywrapper_Unst_FEA_AD.timeout = 1600
     pywrapper_Unst_FEA_AD.tol = 0.00001
diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py
@@ -55,7 +55,7 @@
 
 MESH_FORMAT= RECTANGLE
 MESH_BOX_SIZE= ( 17, 5, 0 )
-MESH_BOX_LENGTH= ( 0.5, 0.05, 0 )
+MESH_BOX_LENGTH= ( 0.5, __HEIGHT__, 0 )
 
 OUTPUT_FILES= RESTART, PARAVIEW
 OUTPUT_WRT_FREQ= 1
@@ -84,6 +84,7 @@
 UNST_ADJOINT_ITER= 21
 ITER_AVERAGE_OBJ= 0
 SOLUTION_FILENAME= restart.dat
+VOLUME_OUTPUT= SENSITIVITY_T0
 """
 
 
@@ -110,12 +111,16 @@ def ApplyLoad(driver, marker_id, peak_load):
   return derivatives
 
 
-def RunPrimal(density, peak_load):
-  """Runs the primal solver for a given density and peak load and returns the time average objective function."""
+def RunPrimal(density, peak_load, height):
+  """
+  Runs the primal solver for a given density, peak load, and beam height.
+  Returns the time average objective function.
+  """
   comm = MPI.COMM_WORLD
 
   with open('config_unsteady.cfg', 'w') as f:
-    f.write(common_settings.replace('__DENSITY__', str(density)) + primal_settings)
+    f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) +
+            primal_settings)
 
   # Initialize the primal driver of SU2, this includes solver preprocessing.
   try:
@@ -144,15 +149,16 @@ def RunPrimal(density, peak_load):
   return tavg_stress_penalty
 
 
-def RunAdjoint(density, peak_load):
+def RunAdjoint(density, peak_load, height):
   """
   Runs the adjoint solver and returns the sensitivity of the objective function to the peak
-  load and to the material density.
+  load, to the material density, and to the beam height.
   """
   comm = MPI.COMM_WORLD
 
   with open('config_unsteady_ad.cfg', 'w') as f:
-    f.write(common_settings.replace('__DENSITY__', str(density)) + adjoint_settings)
+    f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) +
+            adjoint_settings)
 
   # Initialize the adjoint driver of SU2, this includes solver preprocessing.
   try:
@@ -182,7 +188,8 @@ def RunAdjoint(density, peak_load):
     driver.Postprocess()
     driver.Update()
 
-    # Accumulate sensitivies.
+    # Accumulate load sensitivies (the solver doesn't accumulate
+    # these for when they are used for FSI adjoints).
     for i_vertex in range(n_vertex):
       load_sens += derivatives[i_vertex] * driver.GetMarkerFEALoadSensitivity(marker_id, i_vertex)[1]
 
@@ -192,37 +199,51 @@ def RunAdjoint(density, peak_load):
 
   rho_sens = driver.GetOutputValue('SENS_RHO_0')
 
+  height_sens = 0.0
+  sensitivity = driver.Sensitivity(driver.GetSolverIndices()['ADJ.FEA'])
+  coords = driver.Coordinates()
+
+  for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()):
+    y = coords(i_node, 1)
+    dy_dh = y / height
+    height_sens += dy_dh * sensitivity(i_node, 1)
+
   # Finalize the solver and exit cleanly.
   driver.Finalize()
 
-  return comm.allreduce(load_sens), rho_sens
+  return comm.allreduce(load_sens), rho_sens, comm.allreduce(height_sens)
 
 
 def main():
   comm = MPI.COMM_WORLD
   rank = comm.Get_rank()
 
   # Run the primal with 2 loads to compute the sensitivity via finite differences.
-  obj_pert_load = RunPrimal(1, 2.002)
-  obj_pert_rho = RunPrimal(1.0001, 2)
-  obj = RunPrimal(1, 2)
+  obj_pert_height = RunPrimal(1, 2, 0.05001)
+  obj_pert_load = RunPrimal(1, 2.002, 0.05)
+  obj_pert_rho = RunPrimal(1.0001, 2, 0.05)
+  # Run the un-perturbed last to use the restarts for the adjoint.
+  obj = RunPrimal(1, 2, 0.05)
+  sens_height_fd = (obj_pert_height - obj) / 0.00001
   sens_load_fd = (obj_pert_load - obj) / 0.002
   sens_rho_fd = (obj_pert_rho - obj) / 0.0001
 
-  sens_load, sens_rho = RunAdjoint(1, 2)
+  sens_load, sens_rho, sens_height = RunAdjoint(1, 2, 0.05)
 
   if rank == 0:
-    print("      Finite Differences\tDiscrete Adjoint")
-    print(f"Load  {sens_load_fd}\t{sens_load}")
-    print(f"Rho   {sens_rho_fd}\t{sens_rho}")
+    print("       Finite Differences\tDiscrete Adjoint")
+    print(f"Height {sens_height_fd}\t{sens_height}")
+    print(f"Load   {sens_load_fd}\t{sens_load}")
+    print(f"Rho    {sens_rho_fd}\t{sens_rho}")
 
+  assert abs(sens_height / sens_height_fd - 1) < 1e-4, "Error in geometric derivatives."
   assert abs(sens_load / sens_load_fd - 1) < 1e-4, "Error in load derivative."
-  assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in density derivative."
+  assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in material derivative."
 
   # Print results for the regression script to check.
   if rank == 0:
     print("\n")
-    print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho)
+    print(100, 100, sens_load_fd, sens_load, sens_rho_fd, sens_rho, sens_height_fd, sens_height)
 
 
 if __name__ == '__main__':

Original file line number	Diff line number	Diff line change
`@@ -72,6 +72,7 @@ class CDiscAdjVariable : public CVariable {`
`72`	`72`	`inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const final {`
`73`	`73`	`return Sensitivity(iPoint,iDim);`
`74`	`74`	`}`
	`75`	`+ inline const MatrixType& GetSensitivity() const final { return Sensitivity; }`
`75`	`76`
`76`	`77`	`/*!`
`77`	`78`	`* \brief Set/store the dual time contributions to the adjoint variable.`