Merge pull request #1311 from su2code/fix_cp_inv_design

Fix inverse design Cp function

Author: pcarruscag

SU2_CFD/include/output/CFlowOutput.hpp

@@ -77,17 +77,16 @@ class CFlowOutput : public CFVMOutput{
 
   /*!
    * \brief Add CP inverse design output as history fields
-   * \param[in] config - Definition of the particular problem.
    */
-  void Add_CpInverseDesignOutput(CConfig *config);
+  void Add_CpInverseDesignOutput();
 
   /*!
-   * \brief Set CP inverse design output field values
-   * \param[in] solver - The container holding all solution data.
+   * \brief Set CP inverse design output field values (and also into the solver).
+   * \param[in,out] solver - The container holding all solution data.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void Set_CpInverseDesign(CSolver *solver, CGeometry *geometry, CConfig *config);
+  void Set_CpInverseDesign(CSolver *solver, const CGeometry *geometry, const CConfig *config);
 
   /*!
    * \brief Compute value of the Q criteration for vortex idenfitication

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -279,7 +279,7 @@ void CFlowCompOutput::SetHistoryOutputFields(CConfig *config){
 
   /*--- Add Cp diff fields ---*/
 
-  Add_CpInverseDesignOutput(config);
+  Add_CpInverseDesignOutput();
 
 }
 

SU2_CFD/src/output/CFlowOutput.cpp

@@ -695,123 +695,96 @@ void CFlowOutput::SetRotatingFrameCoefficients(CConfig *config, CSolver *flow_so
 }
 
 
-void CFlowOutput::Add_CpInverseDesignOutput(CConfig *config){
+void CFlowOutput::Add_CpInverseDesignOutput(){
 
   AddHistoryOutput("INVERSE_DESIGN_PRESSURE", "Cp_Diff", ScreenOutputFormat::FIXED, "CP_DIFF", "Cp difference for inverse design");
 
 }
 
-void CFlowOutput::Set_CpInverseDesign(CSolver *solver, CGeometry *geometry, CConfig *config){
-
-  unsigned short iMarker, icommas, Boundary;
-  unsigned long iVertex, iPoint, (*Point2Vertex)[2], nPointLocal = 0, nPointGlobal = 0;
-  su2double XCoord, YCoord, ZCoord, Pressure, PressureCoeff = 0, Cp, CpTarget, *Normal = nullptr, Area, PressDiff = 0.0;
-  bool *PointInDomain;
-  string text_line, surfCp_filename;
-  ifstream Surface_file;
+void CFlowOutput::Set_CpInverseDesign(CSolver *solver, const CGeometry *geometry, const CConfig *config){
 
   /*--- Prepare to read the surface pressure files (CSV) ---*/
 
-  surfCp_filename = config->GetUnsteady_FileName("TargetCp", (int)curTimeIter, ".dat");
-
-  /*--- Read the surface pressure file ---*/
+  const auto surfCp_filename = config->GetUnsteady_FileName("TargetCp", curTimeIter, ".dat");
 
-  string::size_type position;
+  /*--- Read the surface pressure file, on the first inner iteration. ---*/
 
+  ifstream Surface_file;
   Surface_file.open(surfCp_filename);
 
-  if (!(Surface_file.fail())) {
-
-    nPointLocal = geometry->GetnPoint();
-    SU2_MPI::Allreduce(&nPointLocal, &nPointGlobal, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
-
-    Point2Vertex = new unsigned long[nPointGlobal][2];
-    PointInDomain = new bool[nPointGlobal];
-
-    for (iPoint = 0; iPoint < nPointGlobal; iPoint ++)
-      PointInDomain[iPoint] = false;
-
-    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-      Boundary   = config->GetMarker_All_KindBC(iMarker);
-
-      if (config->GetSolid_Wall(iMarker) || (Boundary == NEARFIELD_BOUNDARY)) {
-        for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
-
-          /*--- The Pressure file uses the global numbering ---*/
-
-          iPoint = geometry->nodes->GetGlobalIndex(geometry->vertex[iMarker][iVertex]->GetNode());
-
-          if (geometry->vertex[iMarker][iVertex]->GetNode() < geometry->GetnPointDomain()) {
-            Point2Vertex[iPoint][0] = iMarker;
-            Point2Vertex[iPoint][1] = iVertex;
-            PointInDomain[iPoint] = true;
-            solver->SetCPressureTarget(iMarker, iVertex, 0.0);
-          }
-
-        }
-      }
-    }
+  if (!Surface_file.good()) {
+    solver->SetTotal_CpDiff(0.0);
+    SetHistoryOutputValue("INVERSE_DESIGN_PRESSURE", 0.0);
+    return;
+  }
 
+  if ((config->GetInnerIter() == 0) || config->GetDiscrete_Adjoint()) {
+    string text_line;
     getline(Surface_file, text_line);
 
     while (getline(Surface_file, text_line)) {
-      for (icommas = 0; icommas < 50; icommas++) {
-        position = text_line.find( ",", 0 );
-        if (position!=string::npos) text_line.erase (position,1);
-      }
-      stringstream  point_line(text_line);
-
-      if (geometry->GetnDim() == 2) point_line >> iPoint >> XCoord >> YCoord >> Pressure >> PressureCoeff;
-      if (geometry->GetnDim() == 3) point_line >> iPoint >> XCoord >> YCoord >> ZCoord >> Pressure >> PressureCoeff;
-
-      if (PointInDomain[iPoint]) {
-
-        /*--- Find the vertex for the Point and Marker ---*/
-
-        iMarker = Point2Vertex[iPoint][0];
-        iVertex = Point2Vertex[iPoint][1];
-
-        solver->SetCPressureTarget(iMarker, iVertex, PressureCoeff);
-
+      /*--- remove commas ---*/
+      for (auto& c : text_line) if (c == ',') c = ' ';
+      stringstream point_line(text_line);
+
+      /*--- parse line ---*/
+      unsigned long iPointGlobal;
+      su2double XCoord, YCoord, ZCoord=0, Pressure, PressureCoeff;
+
+      point_line >> iPointGlobal >> XCoord >> YCoord;
+      if (nDim == 3) point_line >> ZCoord;
+      point_line >> Pressure >> PressureCoeff;
+
+      const auto iPoint = geometry->GetGlobal_to_Local_Point(iPointGlobal);
+
+      /*--- If the point is on this rank set the Cp to associated vertices
+       *    (one point may be shared by multiple vertices). ---*/
+      if (iPoint >= 0) {
+        bool set = false;
+        for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); ++iMarker) {
+          const auto iVertex = geometry->nodes->GetVertex(iPoint, iMarker);
+
+          if (iVertex >= 0) {
+            solver->SetCPressureTarget(iMarker, iVertex, PressureCoeff);
+            set = true;
+          }
+        }
+        if (!set)
+          cout << "WARNING: In file " << surfCp_filename << ", point " << iPointGlobal << " is not a vertex." << endl;
       }
-
     }
+  }
 
-    Surface_file.close();
-
-    delete [] Point2Vertex;
-    delete [] PointInDomain;
+  /*--- Compute the pressure difference. ---*/
 
-    /*--- Compute the pressure difference ---*/
+  su2double PressDiff = 0.0;
 
-    PressDiff = 0.0;
-    for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-      Boundary   = config->GetMarker_All_KindBC(iMarker);
+  for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); ++iMarker) {
 
-      if (config->GetSolid_Wall(iMarker) || (Boundary == NEARFIELD_BOUNDARY)) {
-        for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
+    const auto Boundary = config->GetMarker_All_KindBC(iMarker);
 
-          Normal = geometry->vertex[iMarker][iVertex]->GetNormal();
+    if (config->GetSolid_Wall(iMarker) || (Boundary == NEARFIELD_BOUNDARY)) {
+      for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) {
 
-          Cp = solver->GetCPressure(iMarker, iVertex);
-          CpTarget = solver->GetCPressureTarget(iMarker, iVertex);
+        const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode();
+        if (!geometry->nodes->GetDomain(iPoint)) continue;
 
-          Area = GeometryToolbox::Norm(nDim, Normal);
+        const auto Cp = solver->GetCPressure(iMarker, iVertex);
+        const auto CpTarget = solver->GetCPressureTarget(iMarker, iVertex);
 
-          PressDiff += Area * (CpTarget - Cp) * (CpTarget - Cp);
-        }
+        const auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal();
+        const auto Area = GeometryToolbox::Norm(nDim, Normal);
 
+        PressDiff += Area * pow(CpTarget-Cp, 2);
       }
     }
-
-    su2double MyPressDiff = PressDiff;
-    SU2_MPI::Allreduce(&MyPressDiff, &PressDiff, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
   }
+  su2double tmp = PressDiff;
+  SU2_MPI::Allreduce(&tmp, &PressDiff, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
 
-  /*--- Update the total Cp difference coeffient ---*/
+  /*--- Update the total Cp difference coeffient. ---*/
 
   solver->SetTotal_CpDiff(PressDiff);
-
   SetHistoryOutputValue("INVERSE_DESIGN_PRESSURE", PressDiff);
 
 }
@@ -820,7 +793,7 @@ su2double CFlowOutput::GetQ_Criterion(su2double** VelocityGradient) const {
 
   /*--- Make a 3D copy of the gradient so we do not have worry about nDim ---*/
 
-  su2double Grad_Vel[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}};
+  su2double Grad_Vel[3][3] = {{0.0}};
 
   for (unsigned short iDim = 0; iDim < nDim; iDim++)
     for (unsigned short jDim = 0 ; jDim < nDim; jDim++)

SU2_CFD/src/output/CNEMOCompOutput.cpp

@@ -284,7 +284,7 @@ void CNEMOCompOutput::SetHistoryOutputFields(CConfig *config){
 
   /*--- Add Cp diff fields ---*/
 
-  Add_CpInverseDesignOutput(config);
+  Add_CpInverseDesignOutput();
 
 }