create matrixview

Author: bigfooted

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -452,6 +452,15 @@ class CDriverBase {
     return CPyWrapperMatrixView(solver->GetNodes()->GetSolution(), "Solution of " + solver->GetSolverName(), false);
   }
 
+  /*!
+   * \brief Get a read/write view of the usrr 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.
    */
@@ -759,24 +768,6 @@ 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.
-   */
-
 /// \}
 
  protected:

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -170,7 +170,6 @@ 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. */
@@ -339,7 +338,7 @@ inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_contain
 
     /*--- Compute the residual for this control volume and subtract. ---*/
     for (iVar = 0; iVar < nVar; iVar++) {
-      LinSysRes(iPoint, iVar) += PointSource(iPoint, iVar) * Volume;
+      LinSysRes(iPoint, iVar) += nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
     }
   }
   END_SU2_OMP_FOR
@@ -2182,17 +2181,7 @@ inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_contain
     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.
@@ -2247,26 +2236,7 @@ inline void Custom_Source_Residual(CGeometry *geometry, CSolver **solver_contain
       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.

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -123,7 +123,6 @@ 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) = su2double(0.0);
 
   /*--- Force definition and coefficient arrays for all of the markers ---*/
 

SU2_CFD/include/solvers/CSolver.hpp

@@ -2852,15 +2852,7 @@ class CSolver {
   inline virtual su2double GetInletFlowDir(unsigned short val_marker,
                                            unsigned long val_vertex,
                                            unsigned short val_dim) const { return 0; }
-  /*!
-   * \brief A virtual member
-   * \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 virtual su2double GetCustomPointSource(unsigned long val_point,
-                                           unsigned short val_var) const { return 0; }
+
   /*!
    * \brief A virtual member
    * \param[in] val_marker - Surface marker where the total temperature is set.
@@ -2892,15 +2884,7 @@ class CSolver {
                                       unsigned long val_vertex,
                                       unsigned short val_dim,
                                       su2double val_flowdir) { }
-  /*!
-   * \brief A virtual member
-   * \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 virtual void SetCustomPointSource(unsigned long val_Point,
-                                      vector<passivedouble> val_source) { }
+
   /*!
    * \brief Updates the components of the farfield velocity vector.
    */

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -40,7 +40,6 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
  protected:
   unsigned short Inlet_Position;             /*!< \brief Column index for scalar variables in inlet files. */
   vector<su2activematrix> Inlet_SpeciesVars; /*!< \brief Species variables at inlet profiles. */
-  su2activematrix SpeciesPointSource;        /*!< \brief User defined source term. */
 
  public:
   /*!
@@ -195,38 +194,4 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
       geometry, solver_container, conv_numerics, visc_numerics, config);
   }
 
-  /*!
-   * \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++) {
-        SpeciesPointSource[val_point][iVar] = val_source[iVar];
-      }
-    }
-  }
-
-/*!
-   * \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 SpeciesPointSource[val_point][val_var];
-  }
-
 };

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). */
@@ -289,6 +291,11 @@ class CVariable {
    */
   void Set_Solution_time_n1();
 
+  /*!
+   * \brief Set the variable solution at time n-1.
+   */
+  void Set_UserDefinedSource();
+
   /*!
    * \brief Set the variable solution at time n.
    * \param[in] iPoint - Point index.
@@ -306,6 +313,15 @@ class CVariable {
       Solution_time_n1(iPoint,iVar) = val_sol[iVar];
   }
 
+  /*!
+   * \brief Set the variable solution at time n-1.
+   * \param[in] iPoint - Point index.
+   */
+  inline void Set_UserDefinedSource(unsigned long iPoint, const su2double* val_sol) {
+    for (unsigned long iVar = 0; iVar < nVar; iVar++)
+      UserDefinedSource(iPoint,iVar) = val_sol[iVar];
+  }
+
   /*!
    * \brief Set the variable solution at time n.
    * \param[in] iPoint - Point index.
@@ -322,6 +338,14 @@ class CVariable {
     Solution_time_n1(iPoint,iVar) = val_sol;
   }
 
+  /*!
+   * \brief Set the variable solution at time n-1.
+   * \param[in] iPoint - Point index.
+   */
+  inline void Set_UserDefinedSource(unsigned long iPoint, unsigned long iVar, su2double val_sol) {
+    UserDefinedSource(iPoint,iVar) = val_sol;
+  }
+
   /*!
    * \brief Virtual Member. Specify a vector to set the velocity components of the solution.
    *        Multiplied by density for compressible cases.
@@ -485,6 +509,20 @@ class CVariable {
    */
   inline su2double *GetSolution(unsigned long iPoint) { return Solution[iPoint]; }
 
+  /*!
+   * \brief Get the entire solution of the problem.
+   * \return Reference to the solution matrix.
+   */
+  inline const MatrixType& GetUserDefinedSource() const { return UserDefinedSource; }
+  inline MatrixType& GetUserDefinedSource() { return UserDefinedSource; }
+
+  /*!
+   * \brief Get the solution of the problem.
+   * \param[in] iPoint - Point index.
+   * \return Pointer to the solution vector.
+   */
+  inline su2double *GetUserDefinedSource(unsigned long iPoint) { return UserDefinedSource[iPoint]; }
+
   /*!
    * \brief Get the old solution of the problem (Runge-Kutta method)
    * \param[in] iPoint - Point index.
@@ -2168,6 +2206,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/solvers/CIncEulerSolver.cpp

@@ -1811,7 +1811,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
     }
   }
 
-  //Custom_Source_Residual(geometry, solver_container, numerics_container, config, iMesh);
+  Custom_Source_Residual(geometry, solver_container, numerics_container, config, iMesh);
 
 }
 

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -110,9 +110,6 @@ void CSpeciesSolver::Initialize(CGeometry* geometry, CConfig* config, unsigned s
 
   nDim = geometry->GetnDim();
 
-  SpeciesPointSource.resize(nPointDomain,nVar) = su2double(0.0);
-
-
   if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) {
 
     /*--- Define some auxiliary vector related with the residual ---*/
@@ -587,10 +584,9 @@ void CSpeciesSolver::Custom_Source_Residual(CGeometry *geometry, CSolver **solve
 
       /*--- 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] -= SpeciesPointSource[iPoint][iVar] * Volume;
+        LinSysRes(iPoint,iVar) -= nodes->GetUserDefinedSource(iPoint)[iVar] * Volume;
       }
   }
   END_SU2_OMP_FOR

SU2_CFD/src/variables/CVariable.cpp

@@ -64,6 +64,10 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
   if (config->GetTime_Marching() != TIME_MARCHING::STEADY)
     Solution_time_n1.resize(nPoint,nVar) = su2double(0.0);
 
+  /*--- User defined source terms ---*/
+  UserDefinedSource.resize(nPoint,nVar) = su2double(0.0);
+
+
   if (config->GetDiscrete_Adjoint()) {
     if (adjoint && config->GetMultizone_Problem())
       External.resize(nPoint,nVar) = su2double(0.0);
@@ -98,6 +102,11 @@ void CVariable::Set_Solution_time_n1() {
   parallelCopy(Solution_time_n.size(), Solution_time_n.data(), Solution_time_n1.data());
 }
 
+void CVariable::Set_UserDefinedSource() {
+  assert(UserDefinedSource.size() == UserDefinedSource.size());
+  parallelCopy(UserDefinedSource.size(), UserDefinedSource.data(), UserDefinedSource.data());
+}
+
 void CVariable::Set_BGSSolution_k() {
   assert(Solution_BGS_k.size() == Solution.size());
   parallelCopy(Solution.size(), Solution.data(), Solution_BGS_k.data());
@@ -127,3 +136,6 @@ void CVariable::RegisterSolution_time_n() {
 void CVariable::RegisterSolution_time_n1() {
   RegisterContainer(true, Solution_time_n1);
 }
+void CVariable::RegisterUserDefinedSource() {
+  RegisterContainer(true, UserDefinedSource);
+}

TestCases/py_wrapper/custom_source_buoyancy/run.py

@@ -108,7 +108,7 @@ def main():
       for iDim in range(nDim):
         custom_source_vector[iDim+1] = -(DensityInc_i - DensityInc_0) * Body_Force_Vector[iDim] / Force_Ref
 
-      driver.SetPointCustomSource(iSOLVER, i_node, custom_source_vector)
+      #driver.SetPointCustomSource(iSOLVER, i_node, custom_source_vector)
 
     print("   ***   inner iteration:",inner_iter)
     driver.Preprocess(inner_iter)