Merge pull request #2017 from su2code/py_wrapper_heat_example

pcarruscag · web-flow · commit f2b75c5b08d9 · 2023-04-30T17:20:21.000-07:00
Python wrapper heat solver example + disc adj minor fixes and cleanup
diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp
@@ -473,6 +473,25 @@ class CDriverBase {
         "MarkerSolutionTimeN of " + solver->GetSolverName(), false);
   }
 
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on all mesh nodes of a solver.
+   */
+  inline CPyWrapperMatrixView SolutionTimeN1(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    return CPyWrapperMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), "SolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on the mesh nodes of a marker.
+   */
+  inline CPyWrapperMarkerMatrixView MarkerSolutionTimeN1(unsigned short iSolver, unsigned short iMarker) {
+    auto* solver = GetSolverAndCheckMarker(iSolver, iMarker);
+    return CPyWrapperMarkerMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker),
+        "MarkerSolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
   /*!
    * \brief Get the flow solver primitive variable names with their associated indices.
    * These correspond to the column indices in the matrix returned by Primitives.
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
@@ -500,6 +500,7 @@ class CVariable {
    * \return Pointer to the solution (at time n-1) vector.
    */
   inline su2double *GetSolution_time_n1(unsigned long iPoint) { return Solution_time_n1[iPoint]; }
+  inline MatrixType& GetSolution_time_n1() { return Solution_time_n1; }
 
   /*!
    * \brief Set the value of the old residual.
diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp
@@ -129,6 +129,10 @@ CDiscAdjSinglezoneDriver::~CDiscAdjSinglezoneDriver() {
 
 void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
+
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Preprocess the adjoint iteration ---*/
diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp
@@ -110,8 +110,10 @@ void CSinglezoneDriver::StartSolver() {
 
 void CSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
-  /*--- Set the current time iteration in the config ---*/
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
 
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Store the current physical time in the config container, as
diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp
diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp
@@ -43,79 +43,64 @@ void CDiscAdjHeatIteration::Preprocess(COutput* output, CIntegration**** integra
   /*--- For the unsteady adjoint, load direct solutions from restart files. ---*/
 
   if (config[val_iZone]->GetTime_Marching() != TIME_MARCHING::STEADY) {
-    const int Direct_Iter = static_cast<int>(config[val_iZone]->GetUnst_AdjointIter()) - static_cast<int>(TimeIter) - 2 + dual_time;
+    const int Direct_Iter = static_cast<int>(config[val_iZone]->GetUnst_AdjointIter()) -
+                            static_cast<int>(TimeIter) - 2 + dual_time;
 
-    /*--- For dual-time stepping we want to load the already converged solution at timestep n ---*/
+    /*--- For dual-time stepping we want to load the already converged solution at previous timesteps.
+     * In general we only load one file and shift the previously loaded solutions, on the first we
+     * load one or two more (depending on dual time order). ---*/
 
-    if (TimeIter == 0) {
-      if (dual_time_2nd) {
-        /*--- Load solution at timestep n-2 ---*/
-
-        LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
-
-        /*--- Push solution back to correct array ---*/
+    if (dual_time_2nd) {
+      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
+    } else if (dual_time_1st) {
+      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 1);
+    }
 
+    if (TimeIter == 0) {
+      /*--- Push solution back one level. ---*/
+      if (dual_time) {
         for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
           solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n();
-          solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1();
         }
       }
-      if (dual_time) {
-        /*--- Load solution at timestep n-1 ---*/
 
+      /*--- If required load another time step. Push the previous time step to n-1 and the
+       loaded time step to n. ---*/
+      if (dual_time_2nd) {
         LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 1);
 
-        /*--- Push solution back to correct array ---*/
-
         for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
+          solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1();
           solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n();
         }
       }
 
-      /*--- Load solution timestep n ---*/
-
+      /*--- Load current solution. ---*/
       LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter);
     }
     if ((TimeIter > 0) && dual_time) {
-      /*--- Load solution timestep n - 2 ---*/
-
-      LoadUnsteady_Solution(geometry, solver, config, val_iZone, val_iInst, Direct_Iter - 2);
-
       /*--- Temporarily store the loaded solution in the Solution_Old array ---*/
-
       for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++)
         solvers[iMesh][HEAT_SOL]->Set_OldSolution();
 
-      /*--- Set Solution at timestep n to solution at n-1 ---*/
-
+      /*--- Move timestep n to current solution. ---*/
       for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
         for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
           solvers[iMesh][HEAT_SOL]->GetNodes()->SetSolution(
               iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n(iPoint));
         }
       }
-      if (dual_time_1st) {
-        /*--- Set Solution at timestep n-1 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-          }
-        }
-      }
-      if (dual_time_2nd) {
-        /*--- Set Solution at timestep n-1 to solution at n-2 ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_time_n1(iPoint));
-          }
-        }
-        /*--- Set Solution at timestep n-2 to the previously loaded solution ---*/
-        for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-          for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
-            solvers[iMesh][HEAT_SOL]->GetNodes()->Set_Solution_time_n1(
-                iPoint, solvers[iMesh][HEAT_SOL]->GetNodes()->GetSolution_Old(iPoint));
+
+      /*--- Finally, place the loaded solution in the correct place (n or n-1 depending on order). ---*/
+      for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
+        auto* heatSol = solvers[iMesh][HEAT_SOL];
+        for (auto iPoint = 0ul; iPoint < geometry[val_iZone][val_iInst][iMesh]->GetnPoint(); iPoint++) {
+          if (dual_time_2nd) {
+            /*--- If required also move timestep n-1 to timestep n. ---*/
+            heatSol->GetNodes()->Set_Solution_time_n(iPoint, heatSol->GetNodes()->GetSolution_time_n1(iPoint));
+            heatSol->GetNodes()->Set_Solution_time_n1(iPoint, heatSol->GetNodes()->GetSolution_Old(iPoint));
+          } else {
+            heatSol->GetNodes()->Set_Solution_time_n(iPoint, heatSol->GetNodes()->GetSolution_Old(iPoint));
           }
         }
       }
@@ -140,23 +125,23 @@ void CDiscAdjHeatIteration::Preprocess(COutput* output, CIntegration**** integra
 void CDiscAdjHeatIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config,
                                                   unsigned short val_iZone, unsigned short val_iInst,
                                                   int val_DirectIter) {
+  auto solvers = solver[val_iZone][val_iInst];
+  auto geometries = geometry[val_iZone][val_iInst];
 
   if (val_DirectIter >= 0) {
     if (rank == MASTER_NODE)
       cout << " Loading heat solution from direct iteration " << val_DirectIter << " for zone " << val_iZone << "." << endl;
 
-    solver[val_iZone][val_iInst][MESH_0][HEAT_SOL]->LoadRestart(
-        geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone], val_DirectIter, false);
+    solvers[MESH_0][HEAT_SOL]->LoadRestart(geometries, solvers, config[val_iZone], val_DirectIter, false);
   }
   else {
     /*--- If there is no solution file we set the freestream condition ---*/
     if (rank == MASTER_NODE)
       cout << " Setting freestream conditions at direct iteration " << val_DirectIter << " for zone " << val_iZone << "." << endl;
 
     for (auto iMesh = 0u; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++) {
-      solver[val_iZone][val_iInst][iMesh][HEAT_SOL]->SetFreeStream_Solution(config[val_iZone]);
-      solver[val_iZone][val_iInst][iMesh][HEAT_SOL]->Postprocessing(
-          geometry[val_iZone][val_iInst][iMesh], solver[val_iZone][val_iInst][iMesh], config[val_iZone], iMesh);
+      solvers[iMesh][HEAT_SOL]->SetFreeStream_Solution(config[val_iZone]);
+      solvers[iMesh][HEAT_SOL]->Postprocessing(geometries[iMesh], solvers[iMesh], config[val_iZone], iMesh);
     }
   }
 }
@@ -177,26 +162,28 @@ void CDiscAdjHeatIteration::InitializeAdjoint(CSolver***** solver, CGeometry****
 
 void CDiscAdjHeatIteration::RegisterInput(CSolver***** solver, CGeometry**** geometry, CConfig** config,
                                           unsigned short iZone, unsigned short iInst, RECORDING kind_recording) {
+  auto solvers0 = solver[iZone][iInst][MESH_0];
+  auto geometry0 = geometry[iZone][iInst][MESH_0];
 
   if (kind_recording == RECORDING::SOLUTION_VARIABLES || kind_recording == RECORDING::SOLUTION_AND_MESH) {
     /*--- Register flow and turbulent variables as input ---*/
 
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->RegisterSolution(geometry0, config[iZone]);
 
-    solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJHEAT_SOL]->RegisterVariables(geometry0, config[iZone]);
   }
   else if (kind_recording == RECORDING::MESH_COORDS) {
     /*--- Register node coordinates as input ---*/
 
-    geometry[iZone][iInst][MESH_0]->RegisterCoordinates();
+    geometry0->RegisterCoordinates();
   }
   else if (kind_recording == RECORDING::MESH_DEFORM) {
     /*--- Register the variables of the mesh deformation ---*/
     /*--- Undeformed mesh coordinates ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterSolution(geometry0, config[iZone]);
 
     /*--- Boundary displacements ---*/
-    solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]);
+    solvers0[ADJMESH_SOL]->RegisterVariables(geometry0, config[iZone]);
   }
 }
 
diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp
@@ -246,19 +246,19 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){
 
   /*--- 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");
+  AddVolumeOutput("SENS_DISP-X", "SensitivityDispN_x", "SENSITIVITY_N", "sensitivity to the previous x displacement");
+  AddVolumeOutput("SENS_DISP-Y", "SensitivityDispN_y", "SENSITIVITY_N", "sensitivity to the previous y displacement");
   if (nDim == 3)
-    AddVolumeOutput("SENS_DISP-Z", "SensInitialDisp_z", "SENSITIVITY_T0", "sensitivity to the initial z displacement");
+    AddVolumeOutput("SENS_DISP-Z", "SensitivityDispN_z", "SENSITIVITY_N", "sensitivity to the previous 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");
+  AddVolumeOutput("SENS_VEL-X", "SensitivityVelN_x", "SENSITIVITY_N", "sensitivity to the previous x velocity");
+  AddVolumeOutput("SENS_VEL-Y", "SensitivityVelN_y", "SENSITIVITY_N", "sensitivity to the previous y velocity");
   if (nDim == 3)
-    AddVolumeOutput("SENS_VEL-Z", "SensInitialVel_z", "SENSITIVITY_T0", "sensitivity to the initial z velocity");
+    AddVolumeOutput("SENS_VEL-Z", "SensitivityVelN_z", "SENSITIVITY_N", "sensitivity to the previous 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");
+  AddVolumeOutput("SENS_ACCEL-X", "SensitivityAccelN_x", "SENSITIVITY_N", "sensitivity to the previous x acceleration");
+  AddVolumeOutput("SENS_ACCEL-Y", "SensitivityAccelN_y", "SENSITIVITY_N", "sensitivity to the previous y acceleration");
   if (nDim == 3)
-    AddVolumeOutput("SENS_ACCEL-Z", "SensInitialAccel_z", "SENSITIVITY_T0", "sensitivity to the initial z acceleration");
+    AddVolumeOutput("SENS_ACCEL-Z", "SensitivityAccelN_z", "SENSITIVITY_N", "sensitivity to the previous z acceleration");
 
 }
diff --git a/SU2_CFD/src/output/CAdjHeatOutput.cpp b/SU2_CFD/src/output/CAdjHeatOutput.cpp
@@ -177,6 +177,14 @@ void CAdjHeatOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("SENSITIVITY", "Surface_Sensitivity", "SENSITIVITY", "sensitivity in normal direction");
   /// END_GROUP
 
+  if (!config->GetTime_Domain()) return;
+
+  /*--- Sensitivities with respect to initial conditions. ---*/
+
+  AddVolumeOutput("SENS_TEMP_N", "SensitivityTempN", "SENSITIVITY_N", "sensitivity to the previous temperature");
+  if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+    AddVolumeOutput("SENS_TEMP_N1", "SensitivityTempN1", "SENSITIVITY_N", "sensitivity to the previous-1 temperature");
+  }
 }
 
 void CAdjHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
@@ -200,6 +208,13 @@ void CAdjHeatOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   if (nDim == 3)
     SetVolumeOutputValue("SENSITIVITY-Z", iPoint, Node_AdjHeat->GetSensitivity(iPoint, 2));
 
+  if (!config->GetTime_Domain()) return;
+
+  SetVolumeOutputValue("SENS_TEMP_N", iPoint, Node_AdjHeat->GetSolution_time_n(iPoint, 0));
+  if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+    SetVolumeOutputValue("SENS_TEMP_N1", iPoint, Node_AdjHeat->GetSolution_time_n1(iPoint, 0));
+  }
+
 }
 
 void CAdjHeatOutput::LoadSurfaceData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint, unsigned short iMarker, unsigned long iVertex){
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
@@ -421,6 +421,19 @@ def main():
     test_list.append(pywrapper_Unst_FEA_AD)
     pass_list.append(pywrapper_Unst_FEA_AD.run_test())
 
+    # Heat solver unsteady AD
+    pywrapper_Unst_Heat_AD = TestCase('pywrapper_Unst_Heat_AD')
+    pywrapper_Unst_Heat_AD.cfg_dir = "py_wrapper/custom_heat_flux"
+    pywrapper_Unst_Heat_AD.cfg_file = "run_ad.py"
+    pywrapper_Unst_Heat_AD.test_iter = 100
+    pywrapper_Unst_Heat_AD.test_vals = [0.776365, 0.776430, 1.000003]
+    pywrapper_Unst_Heat_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py")
+    pywrapper_Unst_Heat_AD.timeout = 1600
+    pywrapper_Unst_Heat_AD.tol = 0.00001
+    pywrapper_Unst_Heat_AD.new_output = False
+    test_list.append(pywrapper_Unst_Heat_AD)
+    pass_list.append(pywrapper_Unst_Heat_AD.run_test())
+
     # Flow AD Mesh Displacement Sensitivity
     pywrapper_CFD_AD_MeshDisp               = TestCase('pywrapper_CFD_AD_MeshDisp')
     pywrapper_CFD_AD_MeshDisp.cfg_dir       = "py_wrapper/disc_adj_flow/mesh_disp_sens"
diff --git a/TestCases/py_wrapper/custom_heat_flux/run_ad.py b/TestCases/py_wrapper/custom_heat_flux/run_ad.py
diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py
