Turbulent combustion (TFC)

Author: bigfooted

Common/include/option_structure.hpp

@@ -549,7 +549,8 @@ 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. */
+  USER_DEFINED = 13,       /*!< \brief user defined through python wrapper. */
 };
 static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("STANDARD_AIR", STANDARD_AIR)
@@ -565,6 +566,7 @@ static const MapType<std::string, ENUM_FLUIDMODEL> FluidModel_Map = {
   MakePair("COOLPROP", COOLPROP)
   MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID)
   MakePair("FLUID_FLAMELET", FLUID_FLAMELET)
+  MakePair("USER_DEFINED", USER_DEFINED)
 };
 
 /*!

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -767,6 +767,9 @@ class CDriverBase {
    */
   void SetPointCustomSource(unsigned short iSolver, unsigned long iPoint, std::vector<passivedouble> values) {
     auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver];
+
+    //if (values[0]>1.0e-6)
+    //  cout << "iPoint="<<iPoint << ", setting custom point source" << values[0]<< endl;
     solver->SetCustomPointSource(iPoint, values);
   }
 
@@ -787,6 +790,23 @@ inline vector<passivedouble> GetSolutionVector(unsigned short iSolver, unsigned
   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);
+  //cout << "setting solution vector " << nodes->GetSolution(iPoint,0) << " " << solutionVector[0] << ", "<< nVar<< endl;
+  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.
@@ -804,6 +824,25 @@ inline vector<passivedouble> GetPrimitiveVector(unsigned short iSolver, unsigned
   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);
+  //cout << "setting primitive vector " << nodes->GetPrimitive(iPoint,0) << " " << primitiveVector[0] << ", "<< nPrimvar<< endl;
+
+  for (auto iVar = 0u; iVar < nPrimvar; ++iVar) {
+    nodes->SetPrimitive(iPoint,iVar, primitiveVector[iVar]);
+  }
+}
+
 /// \}
 
  protected:

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -2225,6 +2225,7 @@ class CFVMFlowSolverBase : public CSolver {
                               vector<passivedouble> val_source) final {
     /*--- Since this call can be accessed indirectly using python, do some error
      * checking to prevent segmentation faults ---*/
+    //cout << "flow" << endl;
     if (val_point > nPointDomain)
       SU2_MPI::Error("Out-of-bounds point index used on solver.", CURRENT_FUNCTION);
     else if (val_source.size() > nVar)

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -124,7 +124,8 @@ void CFVMFlowSolverBase<V, R>::Allocate(const CConfig& config) {
 
   AllocVectorOfMatrices(nVertex, nDim, Inlet_FlowDir);
   PointSource.resize(nPointDomain,nVar);
-
+  PointSource.setConstant(0.0);
+  
   /*--- Force definition and coefficient arrays for all of the markers ---*/
 
   AllocVectorOfVectors(nVertex, CPressure);

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -64,6 +64,7 @@ class CScalarSolver : public CSolver {
   vector<su2matrix<su2double*> > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone)
   vector<vector<int> > SlidingStateNodes;
 
+
   /*--- Shallow copy of grid coloring for OpenMP parallelization. ---*/
 
 #ifdef HAVE_OMP
@@ -578,4 +579,5 @@ class CScalarSolver : public CSolver {
     return SlidingStateNodes[val_marker][val_vertex];
   }
 
+
 };

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -36,6 +36,7 @@ CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config,
                config->GetNEMOProblem(), geometry->GetnDim(), config->GetnSpecies()) {
   nMarker = config->GetnMarker_All();
 
+
   /*--- Store the number of vertices on each marker for deallocation later ---*/
   nVertex.resize(nMarker);
   for (unsigned long iMarker = 0; iMarker < nMarker; iMarker++) nVertex[iMarker] = geometry->nVertex[iMarker];

SU2_CFD/include/solvers/CSpeciesSolver.hpp

@@ -40,6 +40,7 @@ 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 Value of the Flow Direction. */
 
  public:
   /*!
@@ -162,6 +163,16 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
    */
   void Source_Residual(CGeometry* geometry, CSolver** solver_container, CNumerics** numerics_container, CConfig* config,
                        unsigned short iMesh) override;
+/*!
+   * \brief Source term computation for axisymmetric flow.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics_container - Container for 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;
 
 
 /*!
@@ -183,4 +194,42 @@ 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 ---*/
+    //cout << "***** set custom point source species *****" << endl;
+    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];
+        //if (SpeciesPointSource[val_point][iVar] > 1.0e-6)
+        //  cout << iVar<<", setcustompointsource: "<< SpeciesPointSource[val_point][iVar] << endl;
+      }
+    }
+  }
+/*!
+   * \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/src/numerics/flow/flow_sources.cpp

@@ -427,7 +427,6 @@ CNumerics::ResidualType<> CSourceIncBodyForce::ComputeResidual(const CConfig* co
 
   residual[nDim+1] = 0.0;
 
-
   return ResidualType<>(residual, jacobian, nullptr);
 }
 

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -303,6 +303,10 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
 
       config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT/(config->GetMolecular_Weight()/1000.0));
       Pressure_Thermodynamic = Density_FreeStream*Temperature_FreeStream*config->GetGas_Constant();
+      cout << "gas constant = " << config->GetGas_Constant() << endl;
+      cout << "pressure = " << Pressure_Thermodynamic << endl;
+      cout << "rho="<< Density_FreeStream<< endl;
+      cout << "T="<< Temperature_FreeStream<< endl;
       auxFluidModel = new CIncIdealGas(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), Pressure_Thermodynamic);
       auxFluidModel->SetTDState_T(Temperature_FreeStream);
       Pressure_Thermodynamic = auxFluidModel->GetPressure();
@@ -342,6 +346,9 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
       auxFluidModel->SetTDState_T(Temperature_FreeStream, config->GetSpecies_Init());
       break;
 
+    case USER_DEFINED:
+    break;
+
     default:
 
       SU2_MPI::Error("Fluid model not implemented for incompressible solver.", CURRENT_FUNCTION);
@@ -1835,6 +1842,7 @@ void CIncEulerSolver::Custom_Source_Residual(CGeometry *geometry, CSolver **solv
 
       /*--- Compute the residual for this control volume and subtract. ---*/
       for (iVar = 0; iVar < nVar; iVar++) {
+        //cout << iPoint << " " << iVar << ",S="<< PointSource[iPoint][iVar]<< endl;
         LinSysRes[iPoint*nVar+iVar] += PointSource[iPoint][iVar] * Volume;
       }
       // cout << "source = " << iPoint << " " << PointSource[iPoint][0]*Volume 

SU2_CFD/src/solvers/CSpeciesSolver.cpp

@@ -94,6 +94,7 @@ void CSpeciesSolver::Initialize(CGeometry* geometry, CConfig* config, unsigned s
   /*--- Store if an implicit scheme is used, for use during periodic boundary conditions. ---*/
   SetImplicitPeriodic(config->GetKind_TimeIntScheme_Species() == EULER_IMPLICIT);
 
+
   nPrimVar = nVar;
 
   if (nVar > MAXNVAR)
@@ -111,6 +112,12 @@ void CSpeciesSolver::Initialize(CGeometry* geometry, CConfig* config, unsigned s
   nDim = geometry->GetnDim();
 
 
+
+  std::cout << "resize species pointsource, nVar="<<nVar <<" npointdomain="<<nPointDomain << endl;
+  SpeciesPointSource.resize(nPointDomain,nVar);
+  SpeciesPointSource.setConstant(0.0);
+
+
 if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) {
 
     /*--- Define some auxiliary vector related with the residual ---*/
@@ -181,7 +188,7 @@ if (iMesh == MESH_0 || config->GetMGCycle() == FULLMG_CYCLE) {
 void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter,
                                  bool val_update_geo) {
   /*--- Restart the solution from file information ---*/
-
+  cout<<"speciessolver:loadrestart" << endl;
   const string restart_filename = config->GetFilename(config->GetSolution_FileName(), "", val_iter);
 
   /*--- To make this routine safe to call in parallel most of it can only be executed by one thread. ---*/
@@ -569,3 +576,42 @@ void CSpeciesSolver::Source_Residual(CGeometry *geometry, CSolver **solver_conta
     END_SU2_OMP_FOR
   }
 }
+
+
+void CSpeciesSolver::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++) {
+        //if (SpeciesPointSource[iPoint][iVar] > 1.0e-6)
+        //cout << iPoint << " " << iVar << ",Species ="<< SpeciesPointSource[iPoint][iVar]<< endl;
+        LinSysRes[iPoint*nVar+iVar] -= SpeciesPointSource[iPoint][iVar] * Volume;
+      }
+      // cout << "source = " << iPoint << " " << PointSource[iPoint][0]*Volume 
+      //                               << " " << PointSource[iPoint][1]*Volume  
+      //                               << " " << PointSource[iPoint][2]*Volume  
+      //                               << " " << PointSource[iPoint][3]*Volume << endl;
+
+  }
+  END_SU2_OMP_FOR
+
+  AD::EndNoSharedReading();
+
+}