[WIP] pywrapper - custom source terms for all solvers

Common/include/CConfig.hpp

@@ -286,6 +286,7 @@ class CConfig {
   su2double **Giles_FlowDir;                 /*!< \brief Specified flow direction vector (unit vector) for Giles BC. */
   su2double *Inlet_Ptotal;                   /*!< \brief Specified total pressures for inlet boundaries. */
   su2double **Inlet_FlowDir;                 /*!< \brief Specified flow direction vector (unit vector) for inlet boundaries. */
+  //su2double *PointSource;                 /*!< \brief Specified flow direction vector (unit vector) for inlet boundaries. */
   su2double *Inlet_Temperature;              /*!< \brief Specified temperatures for a supersonic inlet boundaries. */
   su2double *Inlet_Pressure;                 /*!< \brief Specified static pressures for supersonic inlet boundaries. */
   su2double **Inlet_Velocity;                /*!< \brief Specified flow velocity vectors for supersonic inlet boundaries. */
@@ -9454,6 +9455,18 @@ class CConfig {
    */
   unsigned long GetnOuter_Iter(void) const { return nOuterIter; }
 
+  /*!
+   * \brief Get the number of inner iterations
+   * \return Number of inner iterations on each multizone block
+   */
+  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
+   */
+  void SetnOuter_Iter(unsigned long val_iter) { cout << "inner iter=" << val_iter << endl; nOuterIter = val_iter; }
+
   /*!
    * \brief Get the number of time iterations
    * \return Number of time steps run

Common/src/CConfig.cpp

@@ -890,6 +890,7 @@ void CConfig::SetPointersNull() {
   Inlet_FlowDir   = nullptr;    Inlet_Temperature = nullptr;    Inlet_Pressure = nullptr;
   Inlet_Velocity  = nullptr;
   Outlet_Pressure = nullptr;    Inlet_SpeciesVal  = nullptr;    Inlet_TurbVal = nullptr;
+  //PointSource = nullptr;
 
   /*--- Engine Boundary Condition settings ---*/
 
@@ -913,7 +914,7 @@ void CConfig::SetPointersNull() {
   Inlet_FlowDir               = nullptr;   Inlet_Temperature       = nullptr;   Inlet_Pressure        = nullptr;
   Inlet_Velocity              = nullptr;   Inflow_Mach             = nullptr;   Inflow_Pressure       = nullptr;
   Outlet_Pressure             = nullptr;   Isothermal_Temperature  = nullptr;
-
+  //PointSource = nullptr;
   ElasticityMod             = nullptr;     PoissonRatio                = nullptr;     MaterialDensity       = nullptr;
 
   Load_Dir = nullptr;            Load_Dir_Value = nullptr;          Load_Dir_Multiplier = nullptr;

SU2_CFD/include/drivers/CDriver.hpp

@@ -475,6 +475,37 @@ 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 Get the number of outer iterations.
+   * \return Number of outer iterations.
+   */
+  unsigned long GetNumberOuterIter() const;
+
+
+  /*!
+   * \brief Set the number of inner iterations.
+   * \return Number of inner iterations.
+   */
+  void SetNumberInnerIter(unsigned long);
+
+  /*!
+   * \brief Set the number of outer iterations.
+   * \return Number of outer iterations.
+   */
+  void SetNumberOuterIter(unsigned long);
+
+  /*!
+   * \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 +586,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,15 @@
    */
   unsigned long GetNumberElements() const;
 
+  /*!
+   * \brief Get the number of solution variables 
+   * \return Number of solution variables.
+   */
+  //unsigned long GetNumberSolverVars() const;
+
+  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 +676,23 @@
     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 +760,49 @@
     }
   }
 
+  /*!
+   * \brief Set the array of variables for the source in the point  
+   * \param[in] iMarker - Marker index.
+   * \param[in] iVertex - Marker vertex index.
+   * \param[in] value - Value of the variable.
+   */
+  void SetPointCustomSource(unsigned short iSOLVER, unsigned long iPoint, std::vector<passivedouble> values) {
+    //for (auto iVar = 0ul; iVar <  nVar; ++iDim) {
+    //GetSolverAndCheckField(FLOW_SOL, iPoint)->SetCustomPointSource (iPoint, value);
+    auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSOLVER];
+    solver->SetCustomPointSource(iPoint, values);
+  }
+
+// return solution vector 
+inline vector<passivedouble> GetSolutionVector(unsigned short iSOLVER, unsigned long iPoint) {
+  auto* solver = solver_container[iSOLVER][INST_0][MESH_0][iSOLVER];
+  auto* nodes = solver->GetNodes();
+  auto nVar = GetNumberSolverVars(iSOLVER);
+  //cout << "getting solution: "<<iSOLVER << " " << iPoint <<  " " << nVar << endl;
+  //auto val = nodes->GetSolution(iPoint,2);
+  //cout << "value : " << val << endl;
+  vector<passivedouble> solutionvector(nVar, 0.0);
+  for (auto iVar = 0u; iVar < nVar; ++iVar) {
+    solutionvector[iVar] = SU2_TYPE::GetValue(nodes->GetSolution(iPoint,iVar));
+    //cout << "vector = " << solutionvector[iVar] << endl;
+  }
+  return solutionvector;
+}
+
+inline vector<passivedouble> GetPrimitiveVector(unsigned short iSOLVER, unsigned long iPoint) {
+  auto* solver = solver_container[iSOLVER][INST_0][MESH_0][iSOLVER];
+  auto* nodes = solver->GetNodes();
+  auto nPrimvar = GetNumberPrimitiveVars(iSOLVER);
+  //cout << "getting solution: "<<iSOLVER << " " << iPoint <<  " " << nVar << endl;
+  //auto val = nodes->GetSolution(iPoint,2);
+  //cout << "value : " << val << endl;
+  vector<passivedouble> solutionvector(nPrimvar, 0.0);
+  for (auto iVar = 0u; iVar < nPrimvar; ++iVar) {
+    solutionvector[iVar] = SU2_TYPE::GetValue(nodes->GetPrimitive(iPoint,iVar));
+    //cout << "vector = " << solutionvector[iVar] << endl;
+  }
+  return solutionvector;
+}
 /// \}
 
  protected:
@@ -750,6 +819,21 @@
     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
+    //if (iPoint < std::numeric_limits<unsigned long>::max() && iPoint > GetNumberMarkers()) {
+     // SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION);
+    //}
+    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/numerics/CNumerics.hpp

@@ -155,6 +155,7 @@ class CNumerics {
   su2double
   LocalGridLength_i; /*!< \brief Local grid length at point i. */
   const su2double *RadVar_Source;  /*!< \brief Source term from the radiative heat transfer equation. */
+  const su2double *UserDefinedSource;  /*!< \brief User Defined Source term. */
   const su2double
   *Coord_i,      /*!< \brief Cartesians coordinates of point i. */
   *Coord_j;      /*!< \brief Cartesians coordinates of point j. */

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -467,3 +467,27 @@ class CSourceRadiation : public CSourceBase_Flow {
   ResidualType<> ComputeResidual(const CConfig* config) override;
 
 };
+
+
+/*!
+ * \class CSourceUserDefined
+ * \brief Class for a user defined source term using the python wrapper
+ * \ingroup SourceDiscr
+ * \author Nijso Beishuizen
+ */
+class CSourceUserDefined : public CSourceBase_Flow {
+private:
+  bool implicit;
+
+public:
+
+  CSourceUserDefined(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config);
+
+  /*!
+   * \brief Source term integration for a user defined source.
+   * \param[in] config - Definition of the particular problem.
+   * \return Lightweight const-view of residual and Jacobian.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override;
+
+};

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -409,6 +409,19 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
                        CNumerics **numerics_container,
                        CConfig *config,
                        unsigned short iMesh) override;
+  /*!
+   * \brief Custom Source term integration.
+   * \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] 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;
 
   /*!
    * \brief Source term integration.

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -170,6 +170,7 @@
   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. */
@@ -2148,6 +2149,18 @@
     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 +2214,29 @@
       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_marker >= nMarker)
+    //  SU2_MPI::Error("Out-of-bounds marker index used on inlet.", CURRENT_FUNCTION);
+    //else if (val_vertex >= nVertex[val_marker])
+    //  SU2_MPI::Error("Out-of-bounds vertex index used on inlet.", CURRENT_FUNCTION);
+    //else
+      //cout << "value list size=" << val_source.size() << endl;
+      //cout << "pointsource.size = " << nVar << endl;
+      //cout << "point,npoint=" << val_point << " " << nPointDomain << endl;
+      for (unsigned short 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,9 @@
   /*--- Store the value of the Flow direction at the inlet BC ---*/
 
   AllocVectorOfMatrices(nVertex, nDim, Inlet_FlowDir);
+
+  //AllocVectorOfVectors(nPointDomain, PointSource);
+  PointSource.resize(nPointDomain,nVar);
 
   /*--- Force definition and coefficient arrays for all of the markers ---*/
 

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -204,6 +204,19 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
                        CNumerics **numerics_container,
                        CConfig *config,
                        unsigned short iMesh) final;
+  /*!
+   * \brief Source term integration.
+   * \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] 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) final;
 
   /*!
    * \brief Source term integration.