pywrapper - custom source terms for all solvers (#2388)

* working version for buoyant flow
* Turbulent combustion (TFC)
* 3D matrix view

---------

Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
Co-authored-by: Pedro Gomes <[email protected]>
Co-authored-by: Pedro Gomes <[email protected]>

Author: 3 people

.github/workflows/regression.yml

@@ -211,7 +211,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c feature_custom_source -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:240320-1536
         with:

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.*/
+  PyCustomSource,                     /*!< \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 GetPyCustomSource(void) const { return PyCustomSource; }
+
   /*!
    * \brief Get the total number of moving markers.
    * \return Total number of moving markers.

Common/include/containers/CPyWrapperMatrixView.hpp

@@ -55,6 +55,9 @@
   /*! \brief Gets the value for a (row, column) pair. */                                                         \
   passivedouble operator()(unsigned long row, unsigned long col) const { return Get(row, col); }                 \
                                                                                                                  \
+  /*! \brief Gets the values for a row of the matrix. */                                                         \
+  std::vector<passivedouble> operator()(unsigned long row) const { return Get(row); }                            \
+                                                                                                                 \
   /*! \brief Gets the value for a (row, column) pair. */                                                         \
   passivedouble Get(unsigned long row, unsigned long col) const { return SU2_TYPE::GetValue(Access(row, col)); } \
                                                                                                                  \
@@ -163,3 +166,77 @@ class CPyWrapperMarkerMatrixView {
   /*--- Use the macro to generate the interface. ---*/
   PY_WRAPPER_MATRIX_INTERFACE
 };
+
+/*!
+ * \class CPyWrapper3DMatrixView
+ * \ingroup PySU2
+ * \brief This class wraps C3DDoubleMatrix for the python wrapper matrix interface.
+ * It is generaly used to wrap access to solver gradients defined for the entire volume.
+ */
+class CPyWrapper3DMatrixView {
+ protected:
+  static_assert(su2activematrix::IsRowMajor, "");
+  su2double* data_ = nullptr;
+  unsigned long rows_ = 0, cols_ = 0, dims_ = 0;
+  std::string name_;
+  bool read_only_ = false;
+
+  /*--- Define the functions required by the interface macro. ---*/
+  inline const su2double& Access(unsigned long row, unsigned long col, unsigned long dim) const {
+    if (row > rows_ || col > cols_ || dim > dims_) SU2_MPI::Error(name_ + " out of bounds", "CPyWrapper3DMatrixView");
+    return data_[row * (cols_ * dims_) + col * dims_ + dim];
+  }
+  inline su2double& Access(unsigned long row, unsigned long col, unsigned long dim) {
+    if (read_only_) SU2_MPI::Error(name_ + " is read-only", "CPyWrapper3DMatrixView");
+    const auto& const_me = *this;
+    return const_cast<su2double&>(const_me.Access(row, col, dim));
+  }
+
+ public:
+  CPyWrapper3DMatrixView() = default;
+
+  /*!
+   * \brief Construct the view of the matrix.
+   * \note "name" should be set to the variable name being returned to give better information to users.
+   * \note "read_only" can be set to true to prevent the data from being modified.
+   */
+  CPyWrapper3DMatrixView(C3DDoubleMatrix& mat, const std::string& name, bool read_only)
+      : data_(mat.data()),
+        rows_(mat.length()),
+        cols_(mat.rows()),
+        dims_(mat.cols()),
+        name_(name),
+        read_only_(read_only) {}
+
+  /*! \brief Returns the shape of the matrix. */
+  std::vector<unsigned long> Shape() const { return {rows_, cols_, dims_}; }
+
+  /*! \brief Returns whether the data is read-only [true] or if it can be modified [false]. */
+  bool IsReadOnly() const { return read_only_; }
+
+  /*! \brief Gets the value for a (row, column, dimension) triplet. */
+  passivedouble operator()(unsigned long row, unsigned long col, unsigned long dim) const { return Get(row, col, dim); }
+
+  /*! \brief Gets the values for a row and column of the matrix. */
+  std::vector<passivedouble> operator()(unsigned long row, unsigned long col) const { return Get(row, col); }
+
+  /*! \brief Gets the value for a (row, column, dimension) triplet. */
+  passivedouble Get(unsigned long row, unsigned long col, unsigned long dim) const {
+    return SU2_TYPE::GetValue(Access(row, col, dim));
+  }
+
+  /*! \brief Gets the values for a row and column of the matrix. */
+  std::vector<passivedouble> Get(unsigned long row, unsigned long col) const {
+    std::vector<passivedouble> vals(dims_);
+    for (unsigned long j = 0; j < dims_; ++j) vals[j] = Get(row, col, j);
+    return vals;
+  }
+  /*! \brief Sets the value for a (row, column, dimension) triplet. This clears derivative information. */
+  void Set(unsigned long row, unsigned long col, unsigned long dim, passivedouble val) { Access(row, col, dim) = val; }
+
+  /*! \brief Sets the values for a row and column of the matrix. */
+  void Set(unsigned long row, unsigned long col, std::vector<passivedouble> vals) {
+    unsigned long j = 0;
+    for (const auto& val : vals) Set(row, col, j++, val);
+  }
+};

Common/include/containers/container_decorators.hpp

@@ -165,6 +165,11 @@ class C3DContainerDecorator {
   FORCEINLINE StaticContainer get(Int i, Index j = 0) const noexcept {
     return m_storage.template get<StaticContainer>(i, j * m_innerSz);
   }
+
+  /*!
+   * \brief Raw data access, for Python wrapper.
+   */
+  FORCEINLINE Scalar* data() { return m_storage.data(); }
 };
 
 /*!

Common/include/geometry/dual_grid/CPoint.hpp

@@ -82,10 +82,9 @@ class CPoint {
   su2activematrix Coord; /*!< \brief vector with the coordinates of the node. */
   su2activematrix
       Coord_Old; /*!< \brief Old coordinates vector for primal solution reloading for Disc.Adj. with dynamic grid. */
-  su2activematrix Coord_Sum; /*!< \brief Sum of coordinates vector for geometry smoothing. */
-  su2activematrix Coord_n;   /*!< \brief Coordinates at time n for use with dynamic meshes. */
-  su2activematrix Coord_n1;  /*!< \brief Coordinates at time n-1 for use with dynamic meshes. */
-  su2activematrix Coord_p1;  /*!< \brief Coordinates at time n+1 for use with dynamic meshes. */
+  su2activematrix Coord_n;  /*!< \brief Coordinates at time n for use with dynamic meshes. */
+  su2activematrix Coord_n1; /*!< \brief Coordinates at time n-1 for use with dynamic meshes. */
+  su2activematrix Coord_p1; /*!< \brief Coordinates at time n+1 for use with dynamic meshes. */
 
   su2activematrix GridVel;      /*!< \brief Velocity of the grid for dynamic mesh cases. */
   CVectorOfMatrix GridVel_Grad; /*!< \brief Gradient of the grid velocity for dynamic meshes. */

Common/include/option_structure.hpp

@@ -550,7 +550,7 @@ enum ENUM_FLUIDMODEL {
   FLUID_MIXTURE = 9,      /*!< \brief Species mixture model. */
   COOLPROP = 10,          /*!< \brief Thermodynamics library. */
   FLUID_FLAMELET = 11,    /*!< \brief lookup table (LUT) method for premixed flamelets. */
-  DATADRIVEN_FLUID = 12,           /*!< \brief multi-layer perceptron driven fluid model. */
+  DATADRIVEN_FLUID = 12,  /*!< \brief multi-layer perceptron driven fluid model. */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)

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", PyCustomSource, 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,
@@ -5600,7 +5603,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
           "to be equal to the number of entries of SPECIES_INIT +1",
           CURRENT_FUNCTION);
 
-    // Helper function that checks scalar variable bounds,
+    /*--- Helper function that checks scalar variable bounds. ---*/
     auto checkScalarBounds = [&](su2double scalar, const string& name, su2double lowerBound, su2double upperBound) {
       if (scalar < lowerBound || scalar > upperBound)
         SU2_MPI::Error(string("Variable: ") + name + string(", is out of bounds."), CURRENT_FUNCTION);

SU2_CFD/include/drivers/CDriver.hpp

@@ -555,6 +555,19 @@ class CDriver : public CDriverBase {
    */
   void SetMarkerTranslationRate(unsigned short iMarker, passivedouble vel_x, passivedouble vel_y, passivedouble vel_z);
 
+  /*!
+   * \brief Get the Freestream Density for nondimensionalization
+   * \return Freestream Density
+   */
+  passivedouble GetDensityFreeStreamND() const;
+
+  /*!
+   * \brief Get the reference Body force for nondimensionalization
+   * \return reference Body Force
+   */
+  passivedouble GetForceRef() const;
+
+
 /// \}
 };
 

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -444,6 +444,23 @@ class CDriverBase {
     return CPyWrapperMatrixView(solver->GetNodes()->GetSolution(), "Solution of " + solver->GetSolverName(), false);
   }
 
+  /*!
+   * \brief Get read/write view of the gradients of a solver variable in a point.
+   */
+  inline CPyWrapper3DMatrixView Gradient(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    return CPyWrapper3DMatrixView(solver->GetNodes()->GetGradient(), "Gradient of " + solver->GetSolverName(), false);
+  }
+
+  /*!
+   * \brief Get a read/write view of the user defined source on all mesh nodes of a solver.
+   */
+  inline CPyWrapperMatrixView UserDefinedSource(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    return CPyWrapperMatrixView(
+      solver->GetNodes()->GetUserDefinedSource(), "User Defined Source of " + solver->GetSolverName(), false);
+  }
+
   /*!
    * \brief Get a read/write view of the current solution on the mesh nodes of a marker.
    */

SU2_CFD/include/solvers/CSolver.hpp

@@ -1557,7 +1557,6 @@ class CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] numerics_container - Description of the numerical method.
-   * \param[in] second_numerics - Description of the second numerical method.
    * \param[in] config - Definition of the particular problem.
    * \param[in] iMesh - Index of the mesh in multigrid computations.
    */
@@ -4334,6 +4333,25 @@ class CSolver {
     END_SU2_OMP_FOR
   }
 
+inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container,
+                                 CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
+
+  AD::StartNoSharedReading();
+
+  SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,2*omp_get_max_threads()))
+  for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) {
+    /*--- Get control volume size. ---*/
+    su2double Volume = geometry->nodes->GetVolume(iPoint);
+    /*--- Compute the residual for this control volume and subtract. ---*/
+    for (auto iVar = 0ul; 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.