[WIP] pywrapper - custom source terms for all solvers

.github/workflows/regression.yml

@@ -209,7 +209,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

@@ -9315,12 +9315,24 @@ class CConfig {
    */
   unsigned long GetnInner_Iter(void) const { return nInnerIter; }
 
+  /*!
+   * \brief Set the number of inner iterations
+   * \return
+   */
+  void SetnInner_Iter(unsigned long val_iter) { nInnerIter = val_iter; }
+
   /*!
    * \brief Get the number of outer iterations
    * \return Number of outer iterations for the multizone problem
    */
   unsigned long GetnOuter_Iter(void) const { return nOuterIter; }
 
+  /*!
+   * \brief Set the number of outer iterations
+   * \return
+   */
+  void SetnOuter_Iter(unsigned long val_iter) { nOuterIter = val_iter; }
+
   /*!
    * \brief Get the number of time iterations
    * \return Number of time steps run

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)

SU2_CFD/include/drivers/CDriver.hpp

@@ -475,6 +475,36 @@ class CDriver : public CDriverBase {
    */
   unsigned long GetNumberTimeIter() const;
 
+  /*!
+   * \brief Get the number of inner iterations.
+   * \return Number of inner iterations.
+   */
+  unsigned long GetNumberInnerIter() const;
+
+  /*!
+   * \brief Set the number of inner iterations.
+   * \return
+   */
+  void SetNumberInnerIter(unsigned long val_iter);
+
+  /*!
+   * \brief Get the number of outer iterations.
+   * \return Number of outer iterations.
+   */
+  unsigned long GetNumberOuterIter() const;
+
+  /*!
+   * \brief Set the number of outer iterations.
+   * \return
+   */
+  void SetNumberOuterIter(unsigned long val_iter);
+
+  /*!
+   * \brief Get the current solution
+   * \return Current solution
+   */
+  unsigned long GetSolution(unsigned short iSolver, unsigned long iPoint, unsigned short iVar);
+
   /*!
    * \brief Get the current time iteration.
    * \return Current time iteration.
@@ -555,6 +585,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
+   */
+  unsigned long GetDensity_FreeStreamND() const;
+
+  /*!
+   * \brief Get the reference Body force for nondimensionalization
+   * \return reference Body Force
+   */
+  unsigned long GetForce_Ref() const;
+
+
 /// \}
 };
 

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -179,6 +179,14 @@ class CDriverBase {
    */
   unsigned long GetNumberElements() const;
 
+  /*!
+   * \brief Get the number of solution variables
+   * \return Number of solution variables.
+   */
+
+  unsigned short GetNumberSolverVars(const unsigned short iSol) const;
+  unsigned short GetNumberPrimitiveVars(const unsigned short iSol) const;
+
   /*!
    * \brief Get the global index of a mesh element.
    * \param[in] iElem - Mesh element index.
@@ -667,6 +675,23 @@ class CDriverBase {
     return sens;
   }
 
+  /*!
+   * \brief Get the gradients of a solver variable in a point.
+   * \returns Vector of gradients grad(iVar).
+   */
+  inline vector<passivedouble> GetGradient(unsigned short iSolver, unsigned long iPoint, unsigned short iVar) {
+    const auto nDim = GetNumberDimensions();
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    auto* nodes = solver->GetNodes();
+
+    vector<passivedouble> grad(nDim, 0.0);
+    for (auto iDim = 0u; iDim < nDim; ++iDim) {
+          grad[iDim] = SU2_TYPE::GetValue(nodes->GetGradient(iPoint, iVar, iDim));
+    }
+    return grad;
+  }
+
+
   /*!
    * \brief Set the adjoint of the structural displacements.
    * \note This can be the input of the FEA solver in an adjoint FSI setting.
@@ -734,6 +759,85 @@ class CDriverBase {
     }
   }
 
+  /*!
+   * \brief Set the array of variables for the source in the point
+   * \param[in] iSolver - Solver index.
+   * \param[in] iPoint - Point index.
+   * \param[in] values - Vector values of the source term.
+   */
+  void SetPointCustomSource(unsigned short iSolver, unsigned long iPoint, std::vector<passivedouble> values) {
+    auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+    solver->SetCustomPointSource(iPoint, values);
+  }
+
+  /*!
+   * \brief Get the solution vector in a point for a specific solver
+   * \param[in] iSolver - Solver index.
+   * \param[in] iPoint - Point index.
+   * \param[out] solutionvector - Vector values of the solution.
+   */
+inline vector<passivedouble> GetSolutionVector(unsigned short iSolver, unsigned long iPoint) {
+  auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+  auto* nodes = solver->GetNodes();
+  auto nVar = GetNumberSolverVars(iSolver);
+  vector<passivedouble> solutionvector(nVar, 0.0);
+  for (auto iVar = 0u; iVar < nVar; ++iVar) {
+    solutionvector[iVar] = SU2_TYPE::GetValue(nodes->GetSolution(iPoint,iVar));
+  }
+  return solutionvector;
+}
+
+/* The vector with actual solution values */
+inline void SetSolutionVector(unsigned short iSolver, unsigned long iPoint, vector<passivedouble> solutionVector) {
+  auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+  auto* nodes = solver->GetNodes();
+
+  /*--- check the size of the solver variables ---*/
+  unsigned short nVar = GetNumberSolverVars(iSolver);
+  if (nVar != solutionVector.size() )
+    SU2_MPI::Error("Solution Vector size is not equal to Solver size.", CURRENT_FUNCTION);
+  for (unsigned int iVar = 0u; iVar < nVar; ++iVar) {
+    nodes->SetSolution(iPoint,iVar, solutionVector[iVar]);
+    nodes->SetSolution_Old(iPoint,iVar, solutionVector[iVar]);
+  }
+}
+
+
+  /*!
+   * \brief Get the primitive variables vector in a point for a specific solver
+   * \param[in] iSolver - Solver index.
+   * \param[in] iPoint - Point index.
+   * \param[out] solutionvector - Vector values of the primitive variables.
+   */
+inline vector<passivedouble> GetPrimitiveVector(unsigned short iSolver, unsigned long iPoint) {
+  auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+  auto* nodes = solver->GetNodes();
+  auto nPrimvar = GetNumberPrimitiveVars(iSolver);
+  vector<passivedouble> solutionvector(nPrimvar, 0.0);
+  for (auto iVar = 0u; iVar < nPrimvar; ++iVar) {
+    solutionvector[iVar] = SU2_TYPE::GetValue(nodes->GetPrimitive(iPoint,iVar));
+  }
+  return solutionvector;
+}
+
+
+  /*!
+   * \brief Get the primitive variables vector in a point for a specific solver
+   * \param[in] iSolver - Solver index.
+   * \param[in] iPoint - Point index.
+   * \param[out] solutionvector - Vector values of the primitive variables.
+   */
+inline void SetPrimitiveVector(unsigned short iSolver, unsigned long iPoint, vector<passivedouble> primitiveVector) {
+  auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+  auto* nodes = solver->GetNodes();
+  auto nPrimvar = GetNumberPrimitiveVars(iSolver);
+  vector<passivedouble> solutionvector(nPrimvar, 0.0);
+
+  for (auto iVar = 0u; iVar < nPrimvar; ++iVar) {
+    nodes->SetPrimitive(iPoint,iVar, primitiveVector[iVar]);
+  }
+}
+
 /// \}
 
  protected:
@@ -750,6 +854,18 @@ class CDriverBase {
     return solver;
   }
 
+  /*!
+   * \brief Automates some boilerplate of accessing solution fields for the python wrapper.
+   */
+  inline CSolver* GetSolverAndCheckField(unsigned short iSolver,
+                                          unsigned long iPoint = std::numeric_limits<unsigned long>::max()) const {
+    // 1. check for the solver the number of variables
+    // 2. check for the mesh the number of points
+    auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+    if (solver == nullptr) SU2_MPI::Error("The selected solver does not exist.", CURRENT_FUNCTION);
+    return solver;
+  }
+
   /*!
    * \brief Initialize containers.
    */

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -170,6 +170,7 @@ class CFVMFlowSolverBase : public CSolver {
   vector<vector<su2double> > Inlet_Ptotal;      /*!< \brief Value of the Total P. */
   vector<vector<su2double> > Inlet_Ttotal;      /*!< \brief Value of the Total T. */
   vector<su2activematrix> Inlet_FlowDir;        /*!< \brief Value of the Flow Direction. */
+  su2activematrix PointSource;        /*!< \brief Value of the Flow Direction. */
   vector<vector<su2double> > HeatFlux;          /*!< \brief Heat transfer coefficient for each boundary and vertex. */
   vector<vector<su2double> > HeatFluxTarget;    /*!< \brief Heat transfer coefficient for each boundary and vertex. */
   vector<su2activematrix> CharacPrimVar;        /*!< \brief Value of the characteristic variables at each boundary. */
@@ -315,6 +316,39 @@ class CFVMFlowSolverBase : public CSolver {
     }
   }
 
+
+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*nVar+iVar] += PointSource[iPoint][iVar] * Volume;
+    }
+  }
+  END_SU2_OMP_FOR
+
+  AD::EndNoSharedReading();
+
+}
+
+
   /*!
    * \brief Computes and sets the required auxilliary vars (and gradients) for axisymmetric flow.
    */
@@ -2148,6 +2182,18 @@ class CFVMFlowSolverBase : public CSolver {
     return Inlet_FlowDir[val_marker][val_vertex][val_dim];
   }
 
+  /*!
+   * \brief A component of the unit vector representing the flow direction at an inlet boundary.
+   * \param[in] val_marker - Surface marker where the flow direction is evaluated
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the flow direction is evaluated
+   * \param[in] val_dim - The component of the flow direction unit vector to be evaluated
+   * \return Component of a unit vector representing the flow direction.
+   */
+  inline su2double GetCustomPointSource(unsigned long val_point,
+                                    unsigned short val_var) const final {
+    return PointSource[val_point][val_var];
+  }
+
   /*!
    * \brief Set the value of the total temperature at an inlet boundary.
    * \param[in] val_marker - Surface marker where the total temperature is set.
@@ -2201,6 +2247,27 @@ class CFVMFlowSolverBase : public CSolver {
       Inlet_FlowDir[val_marker][val_vertex][val_dim] = val_flowdir;
   }
 
+  /*!
+   * \brief Set a component of the unit vector representing the flow direction at an inlet boundary.
+   * \param[in] val_marker - Surface marker where the flow direction is set.
+   * \param[in] val_vertex - Vertex of the marker <i>val_marker</i> where the flow direction is set.
+   * \param[in] val_dim - The component of the flow direction unit vector to be set
+   * \param[in] val_flowdir - Component of a unit vector representing the flow direction.
+   */
+  inline void SetCustomPointSource(unsigned long val_point,
+                              vector<passivedouble> val_source) final {
+    /*--- Since this call can be accessed indirectly using python, do some error
+     * checking to prevent segmentation faults ---*/
+    if (val_point > nPointDomain)
+      SU2_MPI::Error("Out-of-bounds point index used on solver.", CURRENT_FUNCTION);
+    else if (val_source.size() > nVar)
+      SU2_MPI::Error("Out-of-bounds source size used on solver.", CURRENT_FUNCTION);
+    else {
+      for (size_t iVar=0; iVar < val_source.size(); iVar++)
+        PointSource[val_point][iVar] = val_source[iVar];
+    }
+  }
+
   /*!
    * \brief Update the multi-grid structure for the customized boundary conditions.
    * \param geometry_container - Geometrical definition.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -123,6 +123,8 @@ void CFVMFlowSolverBase<V, R>::Allocate(const CConfig& config) {
   /*--- Store the value of the Flow direction at the inlet BC ---*/
 
   AllocVectorOfMatrices(nVertex, nDim, Inlet_FlowDir);
+  PointSource.resize(nPointDomain,nVar);
+  PointSource.setConstant(0.0);
 
   /*--- Force definition and coefficient arrays for all of the markers ---*/