Merge pull request #1705 from su2code/feature_custom_functions

User defined functions for history outputs

Author: pcarruscag

Common/include/CConfig.hpp

@@ -445,6 +445,7 @@ class CConfig {
 
   bool ReorientElements;       /*!< \brief Flag for enabling element reorientation. */
   string CustomObjFunc;        /*!< \brief User-defined objective function. */
+  string CustomOutputs;        /*!< \brief User-defined functions for outputs. */
   unsigned short nDV,                  /*!< \brief Number of design variables. */
   nObj, nObjW;                         /*! \brief Number of objective functions. */
   unsigned short* nDV_Value;           /*!< \brief Number of values for each design variable (might be different than 1 if we allow arbitrary movement). */
@@ -5224,6 +5225,11 @@ class CConfig {
    */
   const string& GetCustomObjFunc() const { return CustomObjFunc; }
 
+  /*!
+   * \brief Get the user expressions for custom outputs.
+   */
+  const string& GetCustomOutputs() const { return CustomOutputs; }
+
   /*!
    * \brief Get the kind of sensitivity smoothing technique.
    * \return Kind of sensitivity smoothing technique.

Common/src/CConfig.cpp

@@ -1932,6 +1932,8 @@ void CConfig::SetConfig_Options() {
 
   /*!\brief CUSTOM_OBJFUNC \n DESCRIPTION: User-provided definition of a custom objective function. \ingroup Config*/
   addStringOption("CUSTOM_OBJFUNC", CustomObjFunc, "");
+  /*!\brief CUSTOM_OUTPUTS \n DESCRIPTION: User-provided definitions for custom output. \ingroup Config*/
+  addStringOption("CUSTOM_OUTPUTS", CustomOutputs, "");
 
   /* DESCRIPTION: parameter for the definition of a complex objective function */
   addDoubleOption("DCD_DCL_VALUE", dCD_dCL, 0.0);

SU2_CFD/include/output/CFlowOutput.hpp

@@ -199,6 +199,84 @@ class CFlowOutput : public CFVMOutput{
    */
   void Set_NearfieldInverseDesign(CSolver *solver, const CGeometry *geometry, const CConfig *config);
 
+  /*!
+   * \brief Compute the custom outputs.
+   * \param[in] solver - The container holding all solution data.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetCustomOutputs(const CSolver* const* solver, const CGeometry *geometry, const CConfig *config);
+
+  /*!
+   * \brief Helper for custom outputs, converts variable names to indices and pointers which are then used
+   * to evaluate the custom expressions.
+   */
+  template <class FlowIndices>
+  void ConvertVariableSymbolsToIndices(const FlowIndices& idx, CustomOutput& output) const {
+
+    static const auto knownVariables =
+        "TEMPERATURE, TEMPERATURE_VE, VELOCITY_X, VELOCITY_Y, VELOCITY_Z, PRESSURE,\n"
+        "DENSITY, ENTHALPY, SOUND_SPEED, LAMINAR_VISCOSITY, EDDY_VISCOSITY, THERMAL_CONDUCTIVITY\n"
+        "TURB[0,1,...], RAD[0,1,...], SPECIES[0,1,...]";
+
+    auto IndexOfVariable = [&](const FlowIndices& idx, const std::string& var) {
+      /*--- Primitives of the flow solver. ---*/
+      const auto flow_offset = FLOW_SOL * CustomOutput::MAX_VARS_PER_SOLVER;
+
+      if ("TEMPERATURE" == var) return flow_offset + idx.Temperature();
+      if ("TEMPERATURE_VE" == var) return flow_offset + idx.Temperature_ve();
+      if ("VELOCITY_X" == var) return flow_offset + idx.Velocity();
+      if ("VELOCITY_Y" == var) return flow_offset + idx.Velocity() + 1;
+      if ("VELOCITY_Z" == var) return flow_offset + idx.Velocity() + 2;
+      if ("PRESSURE" == var) return flow_offset + idx.Pressure();
+      if ("DENSITY" == var) return flow_offset + idx.Density();
+      if ("ENTHALPY" == var) return flow_offset + idx.Enthalpy();
+      if ("SOUND_SPEED" == var) return flow_offset + idx.SoundSpeed();
+      if ("LAMINAR_VISCOSITY" == var) return flow_offset + idx.LaminarViscosity();
+      if ("EDDY_VISCOSITY" == var) return flow_offset + idx.EddyViscosity();
+      if ("THERMAL_CONDUCTIVITY" == var) return flow_offset + idx.ThermalConductivity();
+
+      /*--- Index-based (no name) access to variables of other solvers. ---*/
+      auto GetIndex = [](const std::string& s, int nameLen) {
+        /*--- Extract an int from "name[int]", nameLen is the length of "name". ---*/
+        return std::stoi(std::string(s.begin() + nameLen + 1, s.end() - 1));
+      };
+      if (var.rfind("SPECIES", 0) == 0) return SPECIES_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 7);
+      if (var.rfind("TURB", 0) == 0) return TURB_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 4);
+      if (var.rfind("RAD", 0) == 0) return RAD_SOL * CustomOutput::MAX_VARS_PER_SOLVER + GetIndex(var, 3);
+
+      return CustomOutput::NOT_A_VARIABLE;
+    };
+
+    output.otherOutputs.clear();
+    output.varIndices.clear();
+    output.varIndices.reserve(output.varSymbols.size());
+
+    for (const auto& var : output.varSymbols) {
+      output.varIndices.push_back(IndexOfVariable(idx, var));
+
+      if (output.type == OperationType::FUNCTION && output.varIndices.back() != CustomOutput::NOT_A_VARIABLE) {
+        SU2_MPI::Error("Custom outputs of type 'Function' cannot reference solver variables.", CURRENT_FUNCTION);
+      }
+      /*--- Symbol is a valid solver variable. ---*/
+      if (output.varIndices.back() < CustomOutput::NOT_A_VARIABLE) continue;
+
+      /*--- An index above NOT_A_VARIABLE is not valid with current solver settings. ---*/
+      if (output.varIndices.back() > CustomOutput::NOT_A_VARIABLE) {
+        SU2_MPI::Error("Inactive solver variable (" + var + ") used in function " + output.name + "\n"
+                       "E.g. this may only be a variable of the compressible solver.", CURRENT_FUNCTION);
+      }
+
+      /*--- An index equal to NOT_A_VARIABLE may refer to a history output. ---*/
+      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: " + knownVariables, CURRENT_FUNCTION);
+      }
+    }
+  }
+
   /*!
    * \brief Compute value of the Q criteration for vortex idenfitication
    * \param[in] VelocityGradient - Velocity gradients
@@ -238,6 +316,25 @@ class CFlowOutput : public CFVMOutput{
     return Q;
   }
 
+  /*!
+   * \brief Returns the axisymmetric factor for a point on a marker.
+   */
+  template <class GeoNodes>
+  inline su2double GetAxiFactor(bool axisymmetric, const GeoNodes& nodes, unsigned long iPoint,
+                                unsigned short iMarker) {
+    if (!axisymmetric) return 1.0;
+
+    if (nodes.GetCoord(iPoint, 1) > EPS) return 2 * PI_NUMBER * nodes.GetCoord(iPoint, 1);
+
+    for (const auto jPoint : nodes.GetPoints(iPoint)) {
+      if (nodes.GetVertex(jPoint, iMarker) >= 0) {
+        /*--- Not multiplied by two since we need to half the y coordinate. ---*/
+        return PI_NUMBER * nodes.GetCoord(jPoint, 1);
+      }
+    }
+    return 0.0;
+  }
+
   /*!
    * \brief Write information to meta data file
    * \param[in] config - Definition of the particular problem per zone.

SU2_CFD/include/output/COutput.hpp

@@ -172,7 +172,7 @@ class COutput {
   //! Structure to store the value initial residuals for relative residual computation
   std::map<string, su2double> initialResiduals;
 
-  /** \brief Struct to hold a parsed user-defined expression. */
+  /*! \brief Struct to hold a parsed user-defined expression. */
   struct CustomHistoryOutput {
     mel::ExpressionTree<passivedouble> expression;
     /*--- Pointers to values in the history output maps, to avoid key lookup every time. ---*/
@@ -186,17 +186,62 @@ class COutput {
 
   CustomHistoryOutput customObjFunc;  /*!< \brief User-defined expression for a custom objective. */
 
-   /*----------------------------- Volume output ----------------------------*/
+  /*! \brief Type of operation for custom outputs. */
+  enum class OperationType { MACRO, FUNCTION, AREA_AVG, AREA_INT, MASSFLOW_AVG, MASSFLOW_INT };
 
-   CParallelDataSorter* volumeDataSorter;    //!< Volume data sorter
-   CParallelDataSorter* surfaceDataSorter;   //!< Surface data sorter
+  /*! \brief Struct to hold a parsed custom output function. */
+  struct CustomOutput {
+    /*--- First level of parsing the syntax "name : type{func}[markers];". ---*/
+    std::string name;
+    OperationType type;
+    std::string func;
+    std::vector<std::string> markers;
 
-   vector<string> volumeFieldNames;     //!< Vector containing the volume field names
-   unsigned short nVolumeFields;        /*!< \brief Number of fields in the volume output */
+    /*--- Second level, func into expression, and acceleration structures. ---*/
+    mel::ExpressionTree<passivedouble> expression;
+    std::vector<std::string> varSymbols;
+    std::vector<unsigned short> markerIndices;
+
+    /*--- The symbols (strings) are associated with an integer index for efficiency. For evaluation this index
+     is passed to a functor that returns the value associated with the symbol. This functor is an input to "eval()"
+     and needs to be generated on-the-fly for each point. The functor approach is more generic than a pointer, for
+     example it allows wrapping the access to multiple solvers. The interpretation of these indices is dictated by
+     the functor used in eval, for example indices may be established as 32 * solver_idx + variable_idx.
+     The parts of the code that assign and interpret indices need to be in sync. ---*/
+    std::vector<unsigned long> varIndices;
+
+    /*--- Offset between varIndices of different solvers (see above). Power of 2 to make decoding faster. ---*/
+    static constexpr unsigned long MAX_VARS_PER_SOLVER = 32;
+
+    /*--- Arbitrary number to indicate that a string did not match a variable. ---*/
+    static constexpr unsigned long NOT_A_VARIABLE = MAX_SOLS * MAX_VARS_PER_SOLVER;
+
+    /*--- Other outputs can be referenced in expressions, e.g. to compute variance.
+     We store pointers to the required outputs to speed-up access. ---*/
+    std::vector<const su2double*> otherOutputs;
+
+    /*--- 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. ---*/
+    template <class Variables>
+    su2double eval(const Variables& vars) const {
+      return mel::Eval<su2double>(expression, [&](int iSymbol) {return vars(varIndices[iSymbol]);});
+    }
+  };
+
+  std::vector<CustomOutput> customOutputs;  /*!< \brief User-defined outputs. */
+
+  /*----------------------------- Volume output ----------------------------*/
 
-   string volumeFilename,               //!< Volume output filename
-   surfaceFilename,                     //!< Surface output filename
-   restartFilename;                     //!< Restart output filename
+  CParallelDataSorter* volumeDataSorter;    //!< Volume data sorter
+  CParallelDataSorter* surfaceDataSorter;   //!< Surface data sorter
+
+  vector<string> volumeFieldNames;     //!< Vector containing the volume field names
+  unsigned short nVolumeFields;        //!< Number of fields in the volume output
+
+  string volumeFilename,               //!< Volume output filename
+  surfaceFilename,                     //!< Surface output filename
+  restartFilename;                     //!< Restart output filename
 
   /** \brief Structure to store information for a volume output field.
    *
@@ -628,6 +673,29 @@ class COutput {
     }
   }
 
+  /*!
+   * \brief Returns a pointer to the value of an history output.
+   * \note For per-surface outputs the marker index is specified as "name[index]".
+   */
+  inline const su2double* GetPtrToHistoryOutput(const string& name) const {
+    /*--- Decide if it should be per surface. ---*/
+    const auto pos = name.find('[');
+    const su2double* ptr = nullptr;
+    if (pos == std::string::npos) {
+      const auto it = historyOutput_Map.find(name);
+      if (it != historyOutput_Map.end()) {
+        ptr = &(it->second.value);
+      }
+    } else {
+      const auto idx = std::stoi(std::string(name.begin()+pos+1, name.end()-1));
+      const auto it = historyOutputPerSurface_Map.find(std::string(name, 0, pos));
+      if (it != historyOutputPerSurface_Map.end()) {
+        ptr = &(it->second[idx].value);
+      }
+    }
+    return ptr;
+  }
+
   /*!
    * \brief Setup a custom history output object for a given expression.
    * \param[in] expression - Some user-defined math with the history field names as variables.
@@ -726,6 +794,11 @@ class COutput {
    */
   void SetCommonHistoryFields();
 
+  /*!
+   * \brief Parses user-defined outputs.
+   */
+  void SetCustomOutputs(const CConfig *config);
+
   /*!
    * \brief Load values of the history fields common for all solvers.
    * \param[in] config - Definition of the particular problem.

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -51,6 +51,7 @@ class CIncEulerVariable : public CFlowVariable {
     inline IndexType Temperature() const { return nDim+1; }
     inline IndexType Density() const { return nDim+2; }
     inline IndexType Beta() const { return nDim+3; }
+    inline IndexType SoundSpeed() const { return Beta(); }
     inline IndexType LaminarViscosity() const { return nDim+4; }
     inline IndexType EddyViscosity() const { return nDim+5; }
     inline IndexType ThermalConductivity() const { return nDim+6; }
@@ -59,6 +60,7 @@ class CIncEulerVariable : public CFlowVariable {
 
     /*--- For compatible interface with NEMO. ---*/
     inline IndexType Temperature_ve() const { return std::numeric_limits<IndexType>::max(); }
+    inline IndexType Enthalpy() const { return std::numeric_limits<IndexType>::max(); }
   };
 
  protected:

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -58,6 +58,8 @@ class CNEMOEulerVariable : public CFlowVariable {
     inline IndexType RhoCvve() const {return nSpecies+nDim+7;}
     inline IndexType LaminarViscosity() const {return nSpecies+nDim+8;}
     inline IndexType EddyViscosity() const {return nSpecies+nDim+9;}
+
+    inline IndexType ThermalConductivity() const {return std::numeric_limits<IndexType>::max();}
   };
 
  protected:

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -87,10 +87,7 @@ CFlowCompOutput::CFlowCompOutput(const CConfig *config, unsigned short nDim) : C
   if (convFields.empty() ) convFields.emplace_back("RMS_DENSITY");
 
   if (config->GetFixed_CL_Mode()) {
-    bool found = false;
-    for (unsigned short iField = 0; iField < convFields.size(); iField++)
-      if (convFields[iField] == "LIFT") found = true;
-    if (!found) {
+    if (std::find(convFields.begin(), convFields.end(), "LIFT") != convFields.end()) {
       if (rank == MASTER_NODE)
         cout<<"  Fixed CL: Adding LIFT as Convergence Field to ensure convergence to target CL"<<endl;
       convFields.emplace_back("LIFT");
@@ -446,6 +443,8 @@ void CFlowCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSol
 
   /*--- Keep this as last, since it uses the history values that were set. ---*/
 
+  SetCustomOutputs(solver, geometry, config);
+
   SetCustomAndComboObjectives(FLOW_SOL, config, solver);
 
 }

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -278,6 +278,8 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
 
   /*--- Keep this as last, since it uses the history values that were set. ---*/
 
+  SetCustomOutputs(solver, geometry, config);
+
   SetCustomAndComboObjectives(FLOW_SOL, config, solver);
 
 }