Add simple custom outputs (functions of other outputs) to all solvers

Common/include/CConfig.hpp

@@ -8658,7 +8658,6 @@ class CConfig {
   unsigned short GetKind_Average(void) const { return Kind_Average; }
 
   /*!
-   *
    * \brief Get the direct differentation method.
    * \return direct differentiation method.
    */

SU2_CFD/include/output/CFlowOutput.hpp

@@ -215,7 +215,8 @@ class CFlowOutput : public CFVMOutput{
    * \brief Helper for custom outputs, converts variable names to indices and pointers which are then used
    * to evaluate the custom expressions.
    */
-  void ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, CustomOutput& output) const;
+  void ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, bool allowSkip,
+                                       CustomOutput& output) const;
 
   /*!
    * \brief Compute value of the Q criteration for vortex idenfitication
@@ -296,7 +297,7 @@ class CFlowOutput : public CFVMOutput{
    * \param[in] force_writing - boolean that forces writing of volume output
    * \param[in] iFile - index to the file that we need to consider for volume output
    */
-  bool WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile) override; 
+  bool WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile) override;
   /*!
    * \brief Write the forces breakdown file
    * \param[in] config - Definition of the particular problem per zone.

SU2_CFD/include/output/COutput.hpp

@@ -223,6 +223,10 @@ class COutput {
      We store pointers to the required outputs to speed-up access. ---*/
     std::vector<const su2double*> otherOutputs;
 
+    /*--- For discrete adjoint we may need to skip some expressions because there is one output class
+     for the primal solver and one for the discrete adjoint (each with different variables). ---*/
+    bool skip = false;
+
     /*--- For evaluation, "vars" is a functor (i.e. has operator()) that returns the value of a variable at a given
      point. For example, it can be a wrapper to the primitives pointer, in which case varIndices needs to be setup
      with primitive indices. ---*/
@@ -807,6 +811,14 @@ class COutput {
    */
   void SetCustomOutputs(const CConfig *config);
 
+  /*!
+   * \brief Evaluates function-type custom outputs.
+   * Derived classes can use this to compute simple expressions of other outputs if they
+   * do not implement surface averages. This should be called just before evaluating the
+   * custom objective function.
+   */
+  void ComputeSimpleCustomOutputs(const CConfig *config);
+
   /*!
    * \brief Load values of the history fields common for all solvers.
    * \param[in] config - Definition of the particular problem.

SU2_CFD/src/output/CAdjElasticityOutput.cpp

@@ -172,6 +172,8 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge
       SetHistoryOutputValue("BGS_ADJ_DISP_Z", log10(solver[ADJFEA_SOL]->GetRes_BGS(2)));
     }
   }
+
+  ComputeSimpleCustomOutputs(config);
 }
 
 void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){

SU2_CFD/src/output/CAdjFlowCompOutput.cpp

@@ -225,6 +225,8 @@ void CAdjFlowCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, C
   }
 
   LoadHistoryDataAdjScalar(config, solver);
+
+  ComputeSimpleCustomOutputs(config);
 }
 
 void CAdjFlowCompOutput::SetVolumeOutputFields(CConfig *config) {

SU2_CFD/src/output/CAdjFlowIncOutput.cpp

@@ -275,6 +275,8 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
   }
 
   LoadHistoryDataAdjScalar(config, solver);
+
+  ComputeSimpleCustomOutputs(config);
 }
 
 void CAdjFlowIncOutput::SetVolumeOutputFields(CConfig *config) {

SU2_CFD/src/output/CAdjFlowOutput.cpp

@@ -172,6 +172,8 @@ void CAdjFlowOutput::LoadHistoryDataAdjScalar(const CConfig* config, const CSolv
     SetHistoryOutputValue("LINSOL_ITER_SPECIES", adjspecies_solver->GetIterLinSolver());
     SetHistoryOutputValue("LINSOL_RESIDUAL_SPECIES", log10(adjspecies_solver->GetResLinSolver()));
   }
+
+  ComputeSimpleCustomOutputs(config);
 }
 
 void CAdjFlowOutput::SetVolumeOutputFieldsAdjScalarSolution(const CConfig* config) {

SU2_CFD/src/output/CAdjHeatOutput.cpp

@@ -142,6 +142,7 @@ void CAdjHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
     SetHistoryOutputValue("DEFORM_RESIDUAL", log10(solver[MESH_SOL]->System.GetResidual()));
   }
 
+  ComputeSimpleCustomOutputs(config);
 }
 
 void CAdjHeatOutput::SetVolumeOutputFields(CConfig *config){

SU2_CFD/src/output/CElasticityOutput.cpp

@@ -146,6 +146,9 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
     SetHistoryOutputValue("VOLUME_FRACTION", fea_solver->GetTotal_OFVolFrac());
     SetHistoryOutputValue("TOPOL_DISCRETENESS", fea_solver->GetTotal_OFDiscreteness());
   }
+
+  ComputeSimpleCustomOutputs(config);
+
   /*--- Keep this as last, since it uses the history values that were set. ---*/
   SetCustomAndComboObjectives(FEA_SOL, config, solver);
 

SU2_CFD/src/output/CFlowCompFEMOutput.cpp

@@ -263,6 +263,7 @@ void CFlowCompFEMOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, C
 
   SetAerodynamicCoefficients(config, flow_solver);
 
+  ComputeSimpleCustomOutputs(config);
 }
 
 bool CFlowCompFEMOutput::SetInitResiduals(const CConfig *config){

SU2_CFD/src/output/CFlowOutput.cpp

@@ -739,13 +739,13 @@ void CFlowOutput::SetAnalyzeSurfaceSpeciesVariance(const CSolver* const*solver,
   for (unsigned short iMarker_Analyze = 0; iMarker_Analyze < nMarker_Analyze; iMarker_Analyze++) {
     su2double SpeciesVariance = Surface_SpeciesVariance_Total[iMarker_Analyze];
     SetHistoryOutputPerSurfaceValue("SURFACE_SPECIES_VARIANCE", SpeciesVariance, iMarker_Analyze);
+    config->SetSurface_Species_Variance(iMarker_Analyze, SpeciesVariance);
     Tot_Surface_SpeciesVariance += SpeciesVariance;
-    config->SetSurface_Species_Variance(iMarker_Analyze, Tot_Surface_SpeciesVariance);
   }
   SetHistoryOutputValue("SURFACE_SPECIES_VARIANCE", Tot_Surface_SpeciesVariance);
 }
 
-void CFlowOutput::ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx,
+void CFlowOutput::ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, const bool allowSkip,
                                                   CustomOutput& output) const {
   const auto nameToIndex = PrimitiveNameToIndexMap(idx);
 
@@ -796,23 +796,38 @@ void CFlowOutput::ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsign
     output.varIndices.back() += output.otherOutputs.size();
     output.otherOutputs.push_back(GetPtrToHistoryOutput(var));
     if (output.otherOutputs.back() == nullptr) {
-      SU2_MPI::Error("Invalid history output or solver variable (" + var + ") used in function " + output.name +
-                     "\nValid solvers variables:\n" + knownVariables.str(), CURRENT_FUNCTION);
+      if (!allowSkip) {
+        SU2_MPI::Error("Invalid history output or solver variable (" + var + ") used in function " + output.name +
+                       "\nValid solvers variables:\n" + knownVariables.str(), CURRENT_FUNCTION);
+      } else {
+        if (rank == MASTER_NODE) {
+          std::cout << "Info: Ignoring function " + output.name + " because it may be used by the primal/adjoint "
+                       "solver.\n      If the function is ignored twice it is invalid." << std::endl;
+        }
+        output.skip = true;
+        break;
+      }
     }
   }
 }
 
 void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry *geometry, const CConfig *config) {
 
+  const bool adjoint = config->GetDiscrete_Adjoint();
   const bool axisymmetric = config->GetAxisymmetric();
   const auto* flowNodes = su2staticcast_p<const CFlowVariable*>(solver[FLOW_SOL]->GetNodes());
 
   for (auto& output : customOutputs) {
+    if (output.skip) continue;
+
     if (output.varIndices.empty()) {
+      const bool allowSkip = adjoint && (output.type == OperationType::FUNCTION);
+
       /*--- Setup indices for the symbols in the expression. ---*/
       const auto primIdx = CPrimitiveIndices<unsigned long>(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE,
           config->GetNEMOProblem(), nDim, config->GetnSpecies());
-      ConvertVariableSymbolsToIndices(primIdx, output);
+      ConvertVariableSymbolsToIndices(primIdx, allowSkip, output);
+      if (output.skip) continue;
 
       /*--- Convert marker names to their index (if any) in this rank. Or probe locations to nearest points. ---*/
 

SU2_CFD/src/output/CHeatOutput.cpp

@@ -89,6 +89,8 @@ void CHeatOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver
   SetHistoryOutputValue("LINSOL_RESIDUAL", log10(heat_solver->GetResLinSolver()));
   SetHistoryOutputValue("CFL_NUMBER", config->GetCFL(MESH_0));
 
+  ComputeSimpleCustomOutputs(config);
+
   /*--- Keep this as last, since it uses the history values that were set. ---*/
   SetCustomAndComboObjectives(HEAT_SOL, config, solver);
 }

SU2_CFD/src/output/COutput.cpp

@@ -2135,7 +2135,7 @@ bool COutput::WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_
     return ((Iter % config->GetVolumeOutputFrequency(iFile) == 0)) || force_writing;
   }
       return ((Iter > 0) && (Iter % config->GetVolumeOutputFrequency(iFile) == 0)) || force_writing;
- 
+
 }
 
 void COutput::SetCommonHistoryFields() {
@@ -2287,6 +2287,46 @@ void COutput::SetCustomOutputs(const CConfig* config) {
 
 }
 
+void COutput::ComputeSimpleCustomOutputs(const CConfig *config) {
+  const bool adjoint = config->GetDiscrete_Adjoint();
+
+  for (auto& output : customOutputs) {
+    if (output.type != OperationType::FUNCTION) {
+      if (adjoint) continue;
+      SU2_MPI::Error("The current solver can only use 'Function' custom outputs.", CURRENT_FUNCTION);
+    }
+
+    if (output.varIndices.empty()) {
+      output.varIndices.reserve(output.varSymbols.size());
+      output.otherOutputs.reserve(output.varSymbols.size());
+
+      for (const auto& var : output.varSymbols) {
+        output.varIndices.push_back(output.varIndices.size());
+        output.otherOutputs.push_back(GetPtrToHistoryOutput(var));
+        if (output.otherOutputs.back() == nullptr) {
+          if (!adjoint) {
+            // In primal mode all functions must be valid.
+            SU2_MPI::Error("Invalid history output (" + var + ") used in function " + output.name, CURRENT_FUNCTION);
+          } else {
+            if (rank == MASTER_NODE) {
+              std::cout << "Info: Ignoring function " + output.name + " because it may be used by the primal/adjoint "
+                           "solver.\n      If the function is ignored twice it is invalid." << std::endl;
+            }
+            output.skip = true;
+            break;
+          }
+        }
+      }
+    }
+    if (output.skip) continue;
+
+    auto Functor = [&](unsigned long i) {
+      return *output.otherOutputs[i];
+    };
+    SetHistoryOutputValue(output.name, output.Eval(Functor));
+  }
+}
+
 void COutput::LoadCommonHistoryData(const CConfig *config) {
 
   SetHistoryOutputValue("TIME_STEP", config->GetDelta_UnstTimeND()*config->GetTime_Ref());

TestCases/disc_adj_fsi/Airfoil_2d/config.cfg

@@ -13,6 +13,6 @@ MESH_FILENAME= mesh.su2
 
 OBJECTIVE_FUNCTION= CUSTOM_OBJFUNC
 
-SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], SENS_E_0[1], SENS_NU_0[1]
+SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], sens_e[1], SENS_NU_0[1]
 
 %WRT_ZONE_CONV=YES

TestCases/disc_adj_fsi/Airfoil_2d/configFEA.cfg

@@ -7,7 +7,7 @@ REFERENCE_NODE= 234
 REFERENCE_NODE_DISPLACEMENT= (0.0, 0.0)
 REFERENCE_NODE_PENALTY= 1.0
 DESIGN_VARIABLE_FEA= YOUNG_MODULUS
-CUSTOM_OBJFUNC= '1e2 * (1e3 * REFERENCE_NODE + TOPOL_COMPLIANCE)'
+CUSTOM_OBJFUNC= '1e2 * (ref_node + TOPOL_COMPLIANCE)'
 
 % Solid properties ----------------------------------------------------- %
 MATERIAL_MODEL= NEO_HOOKEAN
@@ -41,7 +41,10 @@ INCREMENTAL_LOAD= YES
 % In\Out --------------------------------------------------------------- %
 MESH_FILENAME= mesh.su2
 MESH_FORMAT= SU2
-
+% Scale SENS_E_0 to avoid a very small value in the regression test.
+% Test using custom outputs for the adjoint and primal solver.
+CUSTOM_OUTPUTS= 'sens_e : Function{1e8 * SENS_E_0};\
+                 ref_node : Function{1e3 * REFERENCE_NODE}'
 RESTART_SOL= NO
 SOLUTION_FILENAME= solution_solid.dat
 SOLUTION_ADJ_FILENAME= adjoint_solid.dat

TestCases/parallel_regression_AD.py

@@ -276,8 +276,8 @@ 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.192713, -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.test_vals         = [-4.349377, 0.192713, -1.303589, 0.75407, 2.3244]
+    discadj_fsi2.test_vals_aarch64 = [-3.479505, 0.127953, -1.303589, 0.75407, 2.3244]
     discadj_fsi2.tol       = 0.00001
     test_list.append(discadj_fsi2)