pywrapper - custom source terms for all solvers

.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.

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -163,8 +163,7 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
   void Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
                        unsigned short iMesh) override;
 
-
-/*!
+  /*!
    * \brief Impose the fluid interface boundary condition using tranfer data.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
@@ -183,4 +182,5 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
       },
       geometry, solver_container, conv_numerics, visc_numerics, config);
   }
+
 };

SU2_CFD/include/variables/CVariable.hpp

@@ -55,6 +55,8 @@ class CVariable {
   MatrixType Solution;       /*!< \brief Solution of the problem. */
   MatrixType Solution_Old;   /*!< \brief Old solution of the problem R-K. */
 
+  MatrixType UserDefinedSource; /*!< \brief User Defined Source of the problem. */
+
   MatrixType External;       /*!< \brief External (outer) contribution in discrete adjoint multizone problems. */
 
   su2vector<bool> Non_Physical;  /*!< \brief Non-physical points in the solution (force first order). */
@@ -485,6 +487,13 @@ class CVariable {
    */
   inline su2double *GetSolution(unsigned long iPoint) { return Solution[iPoint]; }
 
+  /*!
+   * \brief Get the entire User Define Source of the problem.
+   * \return Reference to the solution matrix.
+   */
+  inline const MatrixType& GetUserDefinedSource() const { return UserDefinedSource; }
+  inline MatrixType& GetUserDefinedSource() { return UserDefinedSource; }
+
   /*!
    * \brief Get the old solution of the problem (Runge-Kutta method)
    * \param[in] iPoint - Point index.
@@ -2168,6 +2177,11 @@ class CVariable {
    */
   void RegisterSolution_time_n1();
 
+  /*!
+   * \brief Register the variables in the user defined source array as input/output variable.
+   */
+  void RegisterUserDefinedSource();
+
   /*!
    * \brief Set the adjoint values of the solution.
    * \param[in] adj_sol - The adjoint values of the solution.

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -238,7 +238,7 @@ bool CFluidIteration::Monitor(COutput* output, CIntegration**** integration, CGe
   output->SetHistoryOutput(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0],
                            config[val_iZone], config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(),
                            config[val_iZone]->GetInnerIter());
-  
+
   StopCalc = output->GetConvergence();
 
   /* --- Checking convergence of Fixed CL mode to target CL, and perform finite differencing if needed  --*/
@@ -345,12 +345,12 @@ void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry****
       bladesPrimitives.push_back(bladePrimitives);
     }
     TurbomachineryPerformance->ComputeTurbomachineryPerformance(bladesPrimitives);
-    
+
     auto nSpan = config_container[ZONE_0]->GetnSpanWiseSections();
     auto InState = TurbomachineryPerformance->GetBladesPerformances().at(ZONE_0).at(nSpan)->GetInletState();
     nSpan = config_container[nZone-1]->GetnSpanWiseSections();
     auto OutState =  TurbomachineryPerformance->GetBladesPerformances().at(nZone-1).at(nSpan)->GetOutletState();
-    
+
     TurbomachineryStagePerformance->ComputePerformanceStage(InState, OutState, config_container[nZone-1]);
   }
 }
@@ -658,7 +658,6 @@ void CFluidIteration::SetDualTime_Aeroelastic(CConfig* config) const {
         Marker_Tag = config->GetMarker_All_TagBound(iMarker);
         if (Marker_Tag == Monitoring_Tag) { owner = 1; break;
         }           owner = 0;
-
 
       }
       plunge = config->GetAeroelastic_plunge(iMarker_Monitoring);

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -73,7 +73,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
       sq_vel += Velocity_i *Velocity_i;
     }
 
-    Pressure_i = (Gamma-1.0)*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel);
+    Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel);
     Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0];
 
     residual[0] = yinv*Volume*U_i[2];

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1430,6 +1430,13 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSource(const CConfig* config) {
   /*--- Only place outputs of the "SOURCE" group for scalar transport here. ---*/
 
   switch (config->GetKind_Species_Model()) {
+    case SPECIES_MODEL::SPECIES_TRANSPORT:
+      if (config->GetPyCustomSource()) {
+        for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++){
+          AddVolumeOutput("SPECIES_UDS_" + std::to_string(iVar), "Species_UDS_" + std::to_string(iVar), "SOURCE", "Species User Defined Source " + std::to_string(iVar));
+        }
+      }
+    break;
     case SPECIES_MODEL::FLAMELET: {
       const auto& flamelet_config_options = config->GetFlameletParsedOptions();
       for (auto iCV=0u; iCV < flamelet_config_options.n_control_vars; iCV++) {
@@ -1444,8 +1451,8 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSource(const CConfig* config) {
         const auto& species_name = flamelet_config_options.user_scalar_names[iReactant];
         AddVolumeOutput("SOURCE_" + species_name, "Source_" + species_name, "SOURCE", "Source " + species_name);
       }
-      }
-      break;
+    }
+    break;
     default:
       break;
   }
@@ -1581,12 +1588,15 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
 
     case SPECIES_MODEL::SPECIES_TRANSPORT: {
       const auto Node_Species = solver[SPECIES_SOL]->GetNodes();
-      for (unsigned short iVar = 0; iVar < config->GetnSpecies(); iVar++) {
+      for (unsigned long iVar = 0; iVar < config->GetnSpecies(); iVar++) {
         SetVolumeOutputValue("SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetSolution(iPoint, iVar));
         SetVolumeOutputValue("RES_SPECIES_" + std::to_string(iVar), iPoint, solver[SPECIES_SOL]->LinSysRes(iPoint, iVar));
         SetVolumeOutputValue("DIFFUSIVITY_"+ std::to_string(iVar), iPoint, Node_Species->GetDiffusivity(iPoint,iVar));
         if (config->GetKind_SlopeLimit_Species() != LIMITER::NONE)
           SetVolumeOutputValue("LIMITER_SPECIES_" + std::to_string(iVar), iPoint, Node_Species->GetLimiter(iPoint, iVar));
+        if (config->GetPyCustomSource()){
+          SetVolumeOutputValue("SPECIES_UDS_" + std::to_string(iVar), iPoint, Node_Species->GetUserDefinedSource()(iPoint, iVar));
+        }
       }
       break;
     }

SU2_CFD/src/python_wrapper_structure.cpp

@@ -58,6 +58,14 @@ void CDriver::PreprocessPythonInterface(CConfig** config, CGeometry**** geometry
 
 unsigned long CDriver::GetNumberTimeIter() const { return config_container[selected_zone]->GetnTime_Iter(); }
 
+passivedouble CDriver::GetDensityFreeStreamND() const {
+  return SU2_TYPE::GetValue(config_container[selected_zone]->GetDensity_FreeStreamND());
+  }
+
+passivedouble CDriver::GetForceRef() const {
+  return SU2_TYPE::GetValue(config_container[selected_zone]->GetForce_Ref());
+  }
+
 unsigned long CDriver::GetTimeIter() const { return TimeIter; }
 
 passivedouble CDriver::GetUnsteadyTimeStep() const {