test, place custom_residual in CSolver

Author: bigfooted

Common/include/CConfig.hpp

@@ -707,6 +707,7 @@ class CConfig {
   Wrt_Restart_Overwrite,              /*!< \brief Overwrite restart files or append iteration number.*/
   Wrt_Surface_Overwrite,              /*!< \brief Overwrite surface output files or append iteration number.*/
   Wrt_Volume_Overwrite,               /*!< \brief Overwrite volume output files or append iteration number.*/
+  PyCustom_Source,                    /*!< \brief Use a user defined custom source term .*/
   Restart_Flow;                       /*!< \brief Restart flow solution for adjoint and linearized problems. */
   unsigned short nMarker_Monitoring,  /*!< \brief Number of markers to monitor. */
   nMarker_Designing,                  /*!< \brief Number of markers for the objective function. */
@@ -3087,6 +3088,12 @@ class CConfig {
    */
   unsigned short GetnMarker_PyCustom(void) const { return nMarker_PyCustom; }
 
+  /*!
+   * \brief Get the Python custom source term activation.
+   * \return Custom source term is active or not.
+   */
+  bool GetPyCustom_Source(void) const { return PyCustom_Source; }
+
   /*!
    * \brief Get the total number of moving markers.
    * \return Total number of moving markers.

Common/src/CConfig.cpp

@@ -1548,6 +1548,9 @@ void CConfig::SetConfig_Options() {
   /*!\brief MARKER_PYTHON_CUSTOM\n DESCRIPTION: Python customizable marker(s) \ingroup Config*/
   addStringListOption("MARKER_PYTHON_CUSTOM", nMarker_PyCustom, Marker_PyCustom);
 
+  /*!\brief PYTHON_CUSTOM_SOURCE\n DESCRIPTION: Python custom source \ingroup Config*/
+  addBoolOption("PYTHON_CUSTOM_SOURCE", PyCustom_Source, false);
+
   /*!\brief MARKER_WALL_FUNCTIONS\n DESCRIPTION: Viscous wall markers for which wall functions must be applied.
    Format: (Wall function marker, wall function type, ...) \ingroup Config*/
   addWallFunctionOption("MARKER_WALL_FUNCTIONS", nMarker_WallFunctions, Marker_WallFunctions,

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -316,36 +316,36 @@ class CFVMFlowSolverBase : public CSolver {
   }
 
 
-inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_container,
-                                   CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
+// inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_container,
+//                                    CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
 
-  /*--- Pick one numerics object per thread. ---*/
-  CNumerics* numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num()*MAX_TERMS];
+//   /*--- Pick one numerics object per thread. ---*/
+//   CNumerics* numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num()*MAX_TERMS];
 
-  unsigned short iVar;
-  unsigned long iPoint;
-  AD::StartNoSharedReading();
+//   unsigned short iVar;
+//   unsigned long iPoint;
+//   AD::StartNoSharedReading();
 
-  SU2_OMP_FOR_STAT(omp_chunk_size)
-  for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+//   SU2_OMP_FOR_STAT(omp_chunk_size)
+//   for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
-    /*--- Load the volume of the dual mesh cell ---*/
+//     /*--- Load the volume of the dual mesh cell ---*/
 
-    numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+//     numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
 
-    /*--- Get control volume size. ---*/
-    su2double Volume = geometry->nodes->GetVolume(iPoint);
+//     /*--- Get control volume size. ---*/
+//     su2double Volume = geometry->nodes->GetVolume(iPoint);
 
-    /*--- Compute the residual for this control volume and subtract. ---*/
-    for (iVar = 0; iVar < nVar; iVar++) {
-      LinSysRes(iPoint, iVar) += nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
-    }
-  }
-  END_SU2_OMP_FOR
+//     /*--- Compute the residual for this control volume and subtract. ---*/
+//     for (iVar = 0; iVar < nVar; iVar++) {
+//       LinSysRes(iPoint, iVar) += nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
+//     }
+//   }
+//   END_SU2_OMP_FOR
 
-  AD::EndNoSharedReading();
+//   AD::EndNoSharedReading();
 
-}
+// }
 
 
   /*!

SU2_CFD/include/solvers/CSolver.hpp

@@ -1573,11 +1573,11 @@ class CSolver {
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
-  inline virtual void Custom_Source_Residual(CGeometry *geometry,
-                                       CSolver **solver_container,
-                                       CNumerics **numerics_container,
-                                       CConfig *config,
-                                       unsigned short iMesh) { }
+  // inline virtual void Custom_Source_Residual(CGeometry *geometry,
+  //                                      CSolver **solver_container,
+  //                                      CNumerics **numerics_container,
+  //                                      CConfig *config,
+  //                                      unsigned short iMesh) { }
 
   /*!
    * \brief A virtual member.
@@ -4346,6 +4346,35 @@ class CSolver {
     END_SU2_OMP_FOR
   }
 
+/*--- Why not do this? ---*/
+inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_container,
+                                  CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
+
+  /*--- Pick one numerics object per thread. ---*/
+  CNumerics* numerics = numerics_container[SOURCE_SECOND_TERM + omp_get_thread_num()*MAX_TERMS];
+
+  unsigned short iVar;
+  unsigned long iPoint;
+  AD::StartNoSharedReading();
+
+  SU2_OMP_FOR_STAT(omp_chunk_size)
+  for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+
+    /*--- Load the volume of the dual mesh cell ---*/
+    numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+
+      /*--- Get control volume size. ---*/
+      su2double Volume = geometry->nodes->GetVolume(iPoint);
+      /*--- Compute the residual for this control volume and subtract. ---*/
+      for (iVar = 0; iVar < nVar; iVar++) {
+        LinSysRes(iPoint,iVar) -= base_nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
+      }
+  }
+  END_SU2_OMP_FOR
+
+  AD::EndNoSharedReading();
+}
+
 protected:
   /*!
    * \brief Allocate the memory for the verification solution, if necessary.

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -170,8 +170,8 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
-  void Custom_Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
-                       unsigned short iMesh) override;
+  // void Custom_Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
+  //                      unsigned short iMesh) override;
 
 
 /*!

SU2_CFD/src/integration/CIntegration.cpp

@@ -64,7 +64,7 @@ void CIntegration::Space_Integration(CGeometry *geometry,
   solver_container[MainSolver]->Source_Residual(geometry, solver_container, numerics, config, iMesh);
 
   /*--- Compute custom (python wrapper) source term residuals ---*/
-  solver_container[MainSolver]->Custom_Source_Residual(geometry, solver_container, numerics, config, iMesh);
+  //solver_container[MainSolver]->Custom_Source_Residual(geometry, solver_container, numerics, config, iMesh);
 
   /*--- Add viscous and convective residuals, and compute the Dual Time Source term ---*/
 

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -2053,15 +2053,12 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
   /*--- Pick one numerics object per thread. ---*/
   CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS];
 
-  unsigned short iVar;
-  unsigned long iPoint;
-
   if (body_force) {
 
     /*--- Loop over all points ---*/
     AD::StartNoSharedReading();
     SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
       /*--- Load the conservative variables ---*/
       numerics->SetConservative(nodes->GetSolution(iPoint),
@@ -2091,7 +2088,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
     /*--- Loop over all points ---*/
     AD::StartNoSharedReading();
     SU2_OMP_FOR_DYN(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
       /*--- Load the conservative variables ---*/
       numerics->SetConservative(nodes->GetSolution(iPoint),
@@ -2124,7 +2121,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
     /*--- loop over points ---*/
     AD::StartNoSharedReading();
     SU2_OMP_FOR_DYN(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
       /*--- Set solution  ---*/
       numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));
@@ -2177,7 +2174,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
     /*--- loop over points ---*/
     SU2_OMP_FOR_DYN(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
       /*--- Set solution  ---*/
       numerics->SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(iPoint));
@@ -2200,13 +2197,13 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
     /*--- loop over points ---*/
     SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
       /*--- Get control volume ---*/
       su2double Volume = geometry->nodes->GetVolume(iPoint);
 
       /*--- Get stored time spectral source term and add to residual ---*/
-      for (iVar = 0; iVar < nVar; iVar++) {
+      for (auto iVar = 0; iVar < nVar; iVar++) {
         LinSysRes(iPoint,iVar) += Volume * nodes->GetHarmonicBalance_Source(iPoint,iVar);
       }
     }
@@ -2223,7 +2220,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
     /*--- set vorticity magnitude as auxilliary variable ---*/
     SU2_OMP_FOR_STAT(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPoint; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPoint; iPoint++) {
       const su2double VorticityMag = max(GeometryToolbox::Norm(3, nodes->GetVorticity(iPoint)), 1e-12);
       nodes->SetAuxVar(iPoint, 0, VorticityMag);
     }
@@ -2233,7 +2230,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
     SetAuxVar_Gradient_GG(geometry, config);
 
     SU2_OMP_FOR_DYN(omp_chunk_size)
-    for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
       second_numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr);
       second_numerics->SetVorticity(nodes->GetVorticity(iPoint), nullptr);
       second_numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr);
@@ -2260,7 +2257,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
       /*--- Loop over points ---*/
       SU2_OMP_FOR_DYN(omp_chunk_size)
-      for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
+      for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
 
         /*--- Get control volume size. ---*/
         su2double Volume = geometry->nodes->GetVolume(iPoint);
@@ -2273,14 +2270,34 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain
         VerificationSolution->GetMMSSourceTerm(coor, time, sourceMan.data());
 
         /*--- Compute the residual for this control volume and subtract. ---*/
-        for (iVar = 0; iVar < nVar; iVar++) {
+        for (auto iVar = 0; iVar < nVar; iVar++) {
           LinSysRes(iPoint,iVar) -= sourceMan[iVar]*Volume;
         }
       }
       END_SU2_OMP_FOR
     }
   }
 
+  /*--- Custom user defined source term (from the python wrapper) ---*/
+  if (config->GetPyCustom_Source() ) {
+    SU2_OMP_FOR_STAT(omp_chunk_size)
+    for (auto iPoint = 0; iPoint < nPointDomain; iPoint++) {
+
+      /*--- Load the volume of the dual mesh cell ---*/
+
+      numerics->SetVolume(geometry->nodes->GetVolume(iPoint));
+
+      /*--- Get control volume size. ---*/
+      su2double Volume = geometry->nodes->GetVolume(iPoint);
+
+      /*--- Compute the residual for this control volume and subtract. ---*/
+      for (auto iVar = 0; iVar < nVar; iVar++) {
+        LinSysRes(iPoint, iVar) += nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
+      }
+    }
+    END_SU2_OMP_FOR
+  }
+
   AD::EndNoSharedReading();
 }