Merge pull request #1909 from su2code/feature_probes

Point probes

Author: pcarruscag

SU2_CFD/include/output/COutput.hpp

@@ -187,7 +187,7 @@ class COutput {
   CustomHistoryOutput customObjFunc;  /*!< \brief User-defined expression for a custom objective. */
 
   /*! \brief Type of operation for custom outputs. */
-  enum class OperationType { MACRO, FUNCTION, AREA_AVG, AREA_INT, MASSFLOW_AVG, MASSFLOW_INT };
+  enum class OperationType { MACRO, FUNCTION, AREA_AVG, AREA_INT, MASSFLOW_AVG, MASSFLOW_INT, PROBE };
 
   /*! \brief Struct to hold a parsed custom output function. */
   struct CustomOutput {
@@ -201,6 +201,9 @@ class COutput {
     mel::ExpressionTree<passivedouble> expression;
     std::vector<std::string> varSymbols;
     std::vector<unsigned short> markerIndices;
+    static constexpr long PROBE_NOT_SETUP = -2;
+    static constexpr long PROBE_NOT_OWNED = -1;
+    long iPoint = PROBE_NOT_SETUP;
 
     /*--- 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()"

SU2_CFD/src/output/CFlowOutput.cpp

@@ -25,6 +25,7 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <sstream>
 #include <string>
 
 #include "../../include/output/CFlowOutput.hpp"
@@ -765,16 +766,42 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry
       } else {
         SU2_MPI::Error("Unknown flow solver type.", CURRENT_FUNCTION);
       }
-      /*--- Convert marker names to their index (if any) in this rank. ---*/
-
-      output.markerIndices.clear();
-      for (const auto& marker : output.markers) {
-        for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); ++iMarker) {
-          if (config->GetMarker_All_TagBound(iMarker) == marker) {
-            output.markerIndices.push_back(iMarker);
-            continue;
+      /*--- Convert marker names to their index (if any) in this rank. Or probe locations to nearest points. ---*/
+
+      if (output.type != OperationType::PROBE) {
+        output.markerIndices.clear();
+        for (const auto& marker : output.markers) {
+          for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); ++iMarker) {
+            if (config->GetMarker_All_TagBound(iMarker) == marker) {
+              output.markerIndices.push_back(iMarker);
+              continue;
+            }
+          }
+        }
+      } else {
+        if (output.markers.size() != nDim) {
+          SU2_MPI::Error("Wrong number of coordinates to specify probe " + output.name, CURRENT_FUNCTION);
+        }
+        su2double coord[3] = {};
+        for (auto iDim = 0u; iDim < nDim; ++iDim) coord[iDim] = std::stod(output.markers[iDim]);
+        su2double minDist = std::numeric_limits<su2double>::max();
+        unsigned long minPoint = 0;
+        for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); ++iPoint) {
+          const su2double dist = GeometryToolbox::SquaredDistance(nDim, coord, geometry->nodes->GetCoord(iPoint));
+          if (dist < minDist) {
+            minDist = dist;
+            minPoint = iPoint;
           }
         }
+        /*--- Decide which rank owns the probe. ---*/
+        su2double globMinDist;
+        SU2_MPI::Allreduce(&minDist, &globMinDist, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
+        output.iPoint = fabs(minDist - globMinDist) < EPS ? minPoint : CustomOutput::PROBE_NOT_OWNED;
+        if (output.iPoint != CustomOutput::PROBE_NOT_OWNED) {
+          std::cout << "Probe " << output.name << " is using global point "
+                    << geometry->nodes->GetGlobalIndex(output.iPoint)
+                    << ", distance from target location is " << sqrt(minDist) << std::endl;
+        }
       }
     }
 
@@ -787,6 +814,37 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry
       continue;
     }
 
+    /*--- Prepares the functor that maps symbol indices to values at a given point
+     * (see ConvertVariableSymbolsToIndices). ---*/
+
+    auto MakeFunctor = [&](unsigned long iPoint) {
+      /*--- This returns another lambda that captures iPoint by value. ---*/
+      return [&, iPoint](unsigned long i) {
+        if (i < CustomOutput::NOT_A_VARIABLE) {
+          const auto solIdx = i / CustomOutput::MAX_VARS_PER_SOLVER;
+          const auto varIdx = i % CustomOutput::MAX_VARS_PER_SOLVER;
+          if (solIdx == FLOW_SOL) {
+            return flowNodes->GetPrimitive(iPoint, varIdx);
+          } else {
+            return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx);
+          }
+        } else {
+          return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE];
+        }
+      };
+    };
+
+    if (output.type == OperationType::PROBE) {
+      su2double value = std::numeric_limits<su2double>::max();
+      if (output.iPoint != CustomOutput::PROBE_NOT_OWNED) {
+        value = output.eval(MakeFunctor(output.iPoint));
+      }
+      su2double tmp = value;
+      SU2_MPI::Allreduce(&tmp, &value, 1, MPI_DOUBLE, MPI_MIN, SU2_MPI::GetComm());
+      SetHistoryOutputValue(output.name, value);
+      continue;
+    }
+
     /*--- Surface integral of the expression. ---*/
 
     std::array<su2double, 2> integral = {0.0, 0.0};
@@ -812,23 +870,7 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry
           }
           weight *= GetAxiFactor(axisymmetric, *geometry->nodes, iPoint, iMarker);
           local_integral[1] += weight;
-
-          /*--- Prepare the functor that maps symbol indices to values (see ConvertVariableSymbolsToIndices). ---*/
-
-          auto Functor = [&](unsigned long i) {
-            if (i < CustomOutput::NOT_A_VARIABLE) {
-              const auto solIdx = i / CustomOutput::MAX_VARS_PER_SOLVER;
-              const auto varIdx = i % CustomOutput::MAX_VARS_PER_SOLVER;
-              if (solIdx == FLOW_SOL) {
-                return flowNodes->GetPrimitive(iPoint, varIdx);
-              } else {
-                return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx);
-              }
-            } else {
-              return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE];
-            }
-          };
-          local_integral[0] += weight * output.eval(Functor);
+          local_integral[0] += weight * output.eval(MakeFunctor(iPoint));
         }
         END_SU2_OMP_FOR
       }

SU2_CFD/src/output/COutput.cpp

@@ -2186,6 +2186,7 @@ void COutput::SetCustomOutputs(const CConfig* config) {
     {"AreaInt", OperationType::AREA_INT},
     {"MassFlowAvg", OperationType::MASSFLOW_AVG},
     {"MassFlowInt", OperationType::MASSFLOW_INT},
+    {"Probe", OperationType::PROBE},
   };
   std::stringstream knownOps;
   for (const auto& item : opMap) knownOps << item.first << ", ";

TestCases/parallel_regression.py

@@ -272,7 +272,7 @@ def main():
     flatplate_udobj.cfg_dir   = "user_defined_functions"
     flatplate_udobj.cfg_file  = "lam_flatplate.cfg"
     flatplate_udobj.test_iter = 20
-    flatplate_udobj.test_vals = [-6.653802, -1.181430, -0.794887, 0.000611, -0.000369, 0.000736, -0.001104, 596.690000, 299.800000, 296.890000, 21.492000, 0.563990, 2.278700]
+    flatplate_udobj.test_vals = [-6.653802, -1.181430, -0.794887, 0.000611, -0.000369, 0.000736, -0.001104, 596.690000, 299.800000, 296.890000, 21.492000, 0.563990, 37.148, 2.278700]
     test_list.append(flatplate_udobj)
 
     # Laminar cylinder (steady)

TestCases/user_defined_functions/lam_flatplate.cfg

@@ -24,19 +24,22 @@ RESTART_SOL= NO
 % - AreaAvg and AreaInt: Computes an area average or integral of a field (the
 %   expression) over the list of markers.
 % - MassFlowAvg and MassFlowInt: Computes a mass flow average or integral.
+% - Probe: Evaluates the expression using the values of the mesh point closest
+%   to the coordinates specified inside "[]", [x, y] or [x, y, z] (2 or 3D).
 % NOTE: Each custom output can only use one type, e.g. it is not possible to
 % write 'p_drop : AreaAvg{PRESSURE}[inlet] - AreaAvg{PRESSURE}[outlet]'. This
 % would need to be separated into two AreaAvg outputs and one Function to
 % compute their difference.
 CUSTOM_OUTPUTS= 'velocity : Macro{sqrt(pow(VELOCITY_X, 2) + pow(VELOCITY_Y, 2) + pow(VELOCITY_Z, 2))};\
                  avg_vel : AreaAvg{$velocity}[z_minus, z_plus];\
                  var_vel : AreaAvg{pow($velocity - avg_vel, 2)}[z_minus, z_plus];\
-                 dev_vel : Function{sqrt(var_vel) / avg_vel}'
+                 dev_vel : Function{sqrt(var_vel) / avg_vel};\
+                 probe1 : Probe{$velocity}[0.005, 0.005, 0.05]'
 %
 % "COMBO" is the name and group of the output for the objective function
 % (regardless of definition). "CUSTOM" is the group for all custom outpus.
 SCREEN_OUTPUT= INNER_ITER, RMS_DENSITY, RMS_ENERGY, LINSOL_RESIDUAL, FORCE_Z,\
-               SURFACE_MASSFLOW, SURFACE_TOTAL_TEMPERATURE, avg_vel, dev_vel, COMBO
+               SURFACE_MASSFLOW, SURFACE_TOTAL_TEMPERATURE, avg_vel, dev_vel, probe1, COMBO
 HISTORY_OUTPUT = ITER, AERO_COEFF, FLOW_COEFF, FLOW_COEFF_SURF, CUSTOM, COMBO
 OBJECTIVE_FUNCTION= CUSTOM_OBJFUNC
 % Here we define how the custom objective is computed from other outputs. For