diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index ddf570ab697..59c0e6fa8c9 100644 --- a/Common/include/option_structure.hpp +++ b/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 FluidModel_Map = { MakePair("STANDARD_AIR", STANDARD_AIR) @@ -565,6 +566,7 @@ static const MapType FluidModel_Map = { MakePair("COOLPROP", COOLPROP) MakePair("DATADRIVEN_FLUID", DATADRIVEN_FLUID) MakePair("FLUID_FLAMELET", FLUID_FLAMELET) + MakePair("USER_DEFINED", USER_DEFINED) }; /*! diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 2f79be415b3..cbb646e5def 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -475,6 +475,24 @@ 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 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 +573,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; + + /// \} }; diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp index eaa6857c911..d20de93bc00 100644 --- a/SU2_CFD/include/drivers/CDriverBase.hpp +++ b/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 GetGradient(unsigned short iSolver, unsigned long iPoint, unsigned short iVar) { + const auto nDim = GetNumberDimensions(); + auto* solver = GetSolverAndCheckMarker(iSolver); + auto* nodes = solver->GetNodes(); + + vector 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,90 @@ 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 values) { + auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver]; + + //if (values[0]>1.0e-6) + // cout << "iPoint="<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 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 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 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. + * \param[in] iPoint - Point index. + * \param[out] solutionvector - Vector values of the primitive variables. + */ +inline vector 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 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 primitiveVector) { + auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver]; + auto* nodes = solver->GetNodes(); + auto nPrimvar = GetNumberPrimitiveVars(iSolver); + vector 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: @@ -750,6 +859,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::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. */ diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 2e45916edb1..309cac88e0c 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/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. */ diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index e4cb64e4b54..e6a138ec901 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/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; + +}; diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 83c4cbe18e6..ab5240bb588 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -409,6 +409,19 @@ class CEulerSolver : public CFVMFlowSolverBase > Inlet_Ptotal; /*!< \brief Value of the Total P. */ vector > Inlet_Ttotal; /*!< \brief Value of the Total T. */ vector Inlet_FlowDir; /*!< \brief Value of the Flow Direction. */ + su2activematrix PointSource; /*!< \brief Value of the Flow Direction. */ vector > HeatFlux; /*!< \brief Heat transfer coefficient for each boundary and vertex. */ vector > HeatFluxTarget; /*!< \brief Heat transfer coefficient for each boundary and vertex. */ vector CharacPrimVar; /*!< \brief Value of the characteristic variables at each boundary. */ @@ -2148,6 +2149,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 val_marker 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,28 @@ 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 val_marker 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 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_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. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index d7a634e8227..1f0863cbc69 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -123,7 +123,9 @@ void CFVMFlowSolverBase::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 ---*/ AllocVectorOfVectors(nVertex, CPressure); diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 990c9b7095d..77bc75a49e0 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -204,6 +204,19 @@ class CIncEulerSolver : public CFVMFlowSolverBase > SlidingState; // vector of matrix of pointers... inner dim alloc'd elsewhere (welcome, to the twilight zone) vector > 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]; } + }; diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 1411de2ca78..6d54bc22552 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -36,6 +36,7 @@ CScalarSolver::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]; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 2ed78a78a9b..d4ff5d6298e 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1566,6 +1566,20 @@ class CSolver { CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { } + /*! + * \brief A virtual member. + * \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. + */ + inline virtual void Custom_Source_Residual(CGeometry *geometry, + CSolver **solver_container, + CNumerics **numerics_container, + CConfig *config, + unsigned short iMesh) { } /*! * \brief A virtual member. @@ -2840,7 +2854,15 @@ 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 val_marker 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. @@ -2872,7 +2894,15 @@ 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 val_marker 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 val_source) { } /*! * \brief Updates the components of the farfield velocity vector. */ diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 7236989afd4..fca2cc75c3a 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -40,6 +40,7 @@ class CSpeciesSolver : public CScalarSolver { protected: unsigned short Inlet_Position; /*!< \brief Column index for scalar variables in inlet files. */ vector 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 { */ 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 { }, 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 val_marker 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 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 val_marker 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]; + } + }; diff --git a/SU2_CFD/src/drivers/CDriverBase.cpp b/SU2_CFD/src/drivers/CDriverBase.cpp index 53483856b5d..2a6742fe600 100644 --- a/SU2_CFD/src/drivers/CDriverBase.cpp +++ b/SU2_CFD/src/drivers/CDriverBase.cpp @@ -165,6 +165,9 @@ unsigned long CDriverBase::GetNumberDimensions() const { return main_geometry->G unsigned long CDriverBase::GetNumberElements() const { return main_geometry->GetnElem(); } +unsigned short CDriverBase::GetNumberSolverVars(const unsigned short iSol) const { return solver_container[selected_zone][INST_0][MESH_0][iSol]->GetnVar(); } +unsigned short CDriverBase::GetNumberPrimitiveVars(const unsigned short iSol) const { return solver_container[selected_zone][INST_0][MESH_0][iSol]->GetnPrimVar(); } + unsigned long CDriverBase::GetElementGlobalIndex(unsigned long iElem) const { if (iElem >= GetNumberElements()) { SU2_MPI::Error("Element index exceeds size.", CURRENT_FUNCTION); diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index a3af58d572b..59df24016c4 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -63,6 +63,9 @@ void CIntegration::Space_Integration(CGeometry *geometry, /*--- Compute source term residuals ---*/ solver_container[MainSolver]->Source_Residual(geometry, solver_container, numerics, config, iMesh); + /*--- Compute custom (python wrapper) source term residuals ---*/ + solver_container[MainSolver]->Custom_Source_Residual(geometry, solver_container, numerics, config, iMesh); + /*--- Add viscous and convective residuals, and compute the Dual Time Source term ---*/ if (dual_time) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 6f8359a5f58..ca0f40b104f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -837,3 +837,39 @@ CNumerics::ResidualType<> CSourceRadiation::ComputeResidual(const CConfig *confi return ResidualType<>(residual, jacobian, nullptr); } + +CSourceUserDefined::CSourceUserDefined(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config) : + CSourceBase_Flow(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); +} + +CNumerics::ResidualType<> CSourceUserDefined::ComputeResidual(const CConfig *config) { + + unsigned short iDim; + + /*--- Zero the continuity contribution. ---*/ + + residual[0] = 0.0; + + /*--- Zero the momentum contribution. ---*/ + + for (iDim = 0; iDim < nDim; iDim++) + residual[iDim+1] = 0.0; + + /*--- Set the energy contribution ---*/ + + residual[nDim+1] = -UserDefinedSource[0]*Volume; + + /*--- Set the energy contribution to the Jacobian. ---*/ + + if (implicit) { + + /*--- Jacobian is set to zero on initialization. ---*/ + + jacobian[nDim+1][nDim+1] = -UserDefinedSource[1]*Volume; + + } + + return ResidualType<>(residual, jacobian, nullptr); +} diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 2a148a346d9..1008af2d06a 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -58,6 +58,12 @@ void CDriver::PreprocessPythonInterface(CConfig** config, CGeometry**** geometry unsigned long CDriver::GetNumberTimeIter() const { return config_container[selected_zone]->GetnTime_Iter(); } +unsigned long CDriver::GetNumberInnerIter() const { return config_container[selected_zone]->GetnInner_Iter(); } +unsigned long CDriver::GetNumberOuterIter() const { return config_container[selected_zone]->GetnOuter_Iter(); } + +unsigned long CDriver::GetDensity_FreeStreamND() const { return config_container[selected_zone]->GetDensity_FreeStreamND(); } +unsigned long CDriver::GetForce_Ref() const { return config_container[selected_zone]->GetForce_Ref(); } + unsigned long CDriver::GetTimeIter() const { return TimeIter; } passivedouble CDriver::GetUnsteadyTimeStep() const { @@ -66,6 +72,11 @@ passivedouble CDriver::GetUnsteadyTimeStep() const { string CDriver::GetSurfaceFileName() const { return config_container[selected_zone]->GetSurfCoeff_FileName(); } +unsigned long CDriver::GetSolution(unsigned short iSOLVER, unsigned long iPoint, unsigned short iVar) { + auto solver = solver_container[iZone][INST_0][MESH_0][iSOLVER]; + return solver->GetNodes()->GetSolution(iPoint,iVar); +} + //////////////////////////////////////////////////////////////////////////////// /* Functions related to the management of markers */ //////////////////////////////////////////////////////////////////////////////// diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 6aaf09a5183..9c2bc668ad0 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/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); @@ -1813,6 +1820,44 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont } +void CIncEulerSolver::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++) { + //cout << iPoint << " " << iVar << ",S="<< PointSource[iPoint][iVar]<< endl; + LinSysRes[iPoint*nVar+iVar] += PointSource[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(); + +} + + void CIncEulerSolver::Source_Template(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, unsigned short iMesh) { diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp index 9d7c4494c51..dddf70d199d 100644 --- a/SU2_CFD/src/solvers/CSpeciesSolver.cpp +++ b/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="<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(); + +} diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index 02b91c8ab2c..3a4d45c5ebb 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -191,6 +191,7 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, unsigned sh /*--- Add the solver name. ---*/ SolverName = "SST"; + cout << "end initialize sst" << endl; } void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, diff --git a/TestCases/py_wrapper/custom_initialization/lam_buoyancy_cavity.cfg b/TestCases/py_wrapper/custom_initialization/lam_buoyancy_cavity.cfg new file mode 100644 index 00000000000..470a895b372 --- /dev/null +++ b/TestCases/py_wrapper/custom_initialization/lam_buoyancy_cavity.cfg @@ -0,0 +1,120 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Buoyancy-driven flow inside a cavity % +% Author: Thomas D. Economon % +% Date: 2018.06.10 % +% File Version 8.1.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +% +% Uncomment a density below for a desired Rayleigh number +%INC_DENSITY_INIT= 0.00018903539834 % Ra ~ 1e3 +%INC_DENSITY_INIT= 0.00059778241716 % Ra ~ 1e4 +%INC_DENSITY_INIT= 0.00189035398341 % Ra ~ 1e5 +INC_DENSITY_INIT= 0.00597782417156 % Ra ~ 1e6 +% +INC_VELOCITY_INIT= ( 1.0, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 288.15 + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +FLUID_MODEL= INC_IDEAL_GAS +SPECIFIC_HEAT_CP= 1004.703 +MOLECULAR_WEIGHT= 28.96 + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.716e-5 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY +THERMAL_CONDUCTIVITY_CONSTANT= 0.0246295028571 + +% ----------------------- BODY FORCE DEFINITION -------------------------------% +% +BODY_FORCE= NO +BODY_FORCE_VECTOR= ( 0.0, -9.81, 0.0 ) + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_AREA= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( upper, 0.0, lower, 0.0 ) +MARKER_ISOTHERMAL= ( left, 461.04, right, 115.26 ) +MARKER_PLOTTING= ( upper, left, right, lower ) +MARKER_MONITORING= ( NONE ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 50 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 15.0, 1e10) +MAX_DELTA_TIME= 1E6 +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER=1 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 1E-15 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_cavity_65x65.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FILENAME= solution_flow.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow.dat +RESTART_ADJ_FILENAME= restart_adj.dat +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 100 +SCREEN_OUTPUT= (OUTER_ITER, INNER_ITER, RMS_PRESSURE, RMS_TEMPERATURE, LIFT, DRAG) + +% ----------------------- GEOMETRY EVALUATION PARAMETERS ----------------------% +GEO_BOUNDS= ( 0.499, 0.501) diff --git a/TestCases/py_wrapper/custom_initialization/run.py b/TestCases/py_wrapper/custom_initialization/run.py new file mode 100644 index 00000000000..757b464a9c9 --- /dev/null +++ b/TestCases/py_wrapper/custom_initialization/run.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Buoyancy force using user defines source term +# \version 8.1.0 "Harrier" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import sys +import pysu2 +# from mpi4py import MPI + +def main(): + """ + Run the flow solver with a custom inlet (function of time and space). + """ + # comm = MPI.COMM_WORLD + comm = 0 + + # Initialize the primal driver of SU2, this includes solver preprocessing. + try: + driver = pysu2.CSinglezoneDriver('lam_buoyancy_cavity.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CSinglezoneDriver : ', exception) + raise + + print("\n------------------------------ Begin Solver -----------------------------") + sys.stdout.flush() + + # we need to add a source term to the energy equation. For this, we need to get the solver, and the variable first. + # we then loop over all points and for these points, we add the source term + nDim = driver.GetNumberDimensions() + + # index to the flow solver + # C.FLOW + # INC.FLOW + # HEAT + # FLAMELET + # SPECIES + # SA + # SST + iSOLVER = driver.GetSolverIndices()['INC.FLOW'] + print("index of flow solver = ",iSOLVER) + + # all the indices and the map to the names of the primitives + primindex = driver.GetPrimitiveIndices() + print("indices of primitives=",primindex) + print("number of primitives:",len(primindex)) + + nElem = driver.GetNumberElements() + print("number of elements:",nElem) + + nVars = driver.GetNumberSolverVars(iSOLVER) + print("number of solver variables:",nVars) + varindex = primindex.copy() + for prim in varindex.copy(): + if varindex[prim] >=nVars: + del varindex[prim] + varindex = dict(sorted(varindex.items(), key=lambda item: item[1])) + + +# it is possible to get the solver type by doing +# iSOLVER = driver.GetSolverIndices()['INC.FLOW'] +# with the solver type we then get the solver variables using: + +# nVars = driver.GetNumberSolverVars(iSOLVER) +# and the solver variable names: +# print("solver variable names:",varindex) + +# we can overwrite the solution using: +# driver.SetSolutionVector(iSolver, iPoint, solutionVector) +# + + print("solver variable names:",varindex) + iDENSITY = primindex.get("DENSITY") + print("index of density = ",iDENSITY) + + index_Vel = varindex.get("VELOCITY_X") + print("index of velocity = ",index_Vel) + custom_source_vector = [0.0 for i in range(nVars)] + print("custom source vector = ", custom_source_vector) + + print("max. number of inner iterations: ",driver.GetNumberInnerIter()); + print("max nr of outer iterations: ",driver.GetNumberOuterIter()); + + # is in domain: isdomain = driver.GetNodeDomain(iPoint) + #for i_vertex in range(n_vertex) + #AllSolutions = driver.GetAllSolutions(iSOLVER) + Body_Force_Vector = [0.0, -9.81, 0.0] + DensityInc_0 = driver.GetDensity_FreeStreamND() + print("rho freestream = ",DensityInc_0) + Force_Ref = driver.GetForce_Ref() + print("reference force = ",Force_Ref) + + for inner_iter in range(500): + + # set the source term, per point + for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()): + #print("node = ",i_node) + #SolutionVector = driver.GetSolutionVector(iSOLVER,i_node) + #print("solutionvector=",SolutionVector) + PrimitiveVector = driver.GetPrimitiveVector(iSOLVER,i_node) + #print("primitivevector=",PrimitiveVector) + DensityInc_i = PrimitiveVector[iDENSITY] + + for iDim in range(nDim): + custom_source_vector[iDim+1] = -(DensityInc_i - DensityInc_0) * Body_Force_Vector[iDim] / Force_Ref + + #print("density=",DensityInc_i) + driver.SetPointCustomSource(iSOLVER, i_node,custom_source_vector) + + print(" *** inner iteration:",inner_iter) + driver.Preprocess(inner_iter) + + # Run one time iteration. + driver.Run() + driver.Postprocess() + driver.Update() + + # Monitor the solver and output solution to file if required. + #driver.Monitor(inner_iter) + driver.Output(inner_iter) + + # Finalize the solver and exit cleanly. + driver.Finalize() + +if __name__ == '__main__': + main() diff --git a/TestCases/py_wrapper/custom_source/lam_buoyancy_cavity.cfg b/TestCases/py_wrapper/custom_source/lam_buoyancy_cavity.cfg new file mode 100644 index 00000000000..470a895b372 --- /dev/null +++ b/TestCases/py_wrapper/custom_source/lam_buoyancy_cavity.cfg @@ -0,0 +1,120 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Buoyancy-driven flow inside a cavity % +% Author: Thomas D. Economon % +% Date: 2018.06.10 % +% File Version 8.1.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= VARIABLE +INC_ENERGY_EQUATION = YES +% +% Uncomment a density below for a desired Rayleigh number +%INC_DENSITY_INIT= 0.00018903539834 % Ra ~ 1e3 +%INC_DENSITY_INIT= 0.00059778241716 % Ra ~ 1e4 +%INC_DENSITY_INIT= 0.00189035398341 % Ra ~ 1e5 +INC_DENSITY_INIT= 0.00597782417156 % Ra ~ 1e6 +% +INC_VELOCITY_INIT= ( 1.0, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 288.15 + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +FLUID_MODEL= INC_IDEAL_GAS +SPECIFIC_HEAT_CP= 1004.703 +MOLECULAR_WEIGHT= 28.96 + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.716e-5 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY +THERMAL_CONDUCTIVITY_CONSTANT= 0.0246295028571 + +% ----------------------- BODY FORCE DEFINITION -------------------------------% +% +BODY_FORCE= NO +BODY_FORCE_VECTOR= ( 0.0, -9.81, 0.0 ) + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_AREA= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( upper, 0.0, lower, 0.0 ) +MARKER_ISOTHERMAL= ( left, 461.04, right, 115.26 ) +MARKER_PLOTTING= ( upper, left, right, lower ) +MARKER_MONITORING= ( NONE ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 50 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 1.5, 0.5, 15.0, 1e10) +MAX_DELTA_TIME= 1E6 +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER=1 + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 1E-15 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_cavity_65x65.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FILENAME= solution_flow.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +TABULAR_FORMAT= CSV +CONV_FILENAME= history +RESTART_FILENAME= restart_flow.dat +RESTART_ADJ_FILENAME= restart_adj.dat +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +OUTPUT_WRT_FREQ= 100 +SCREEN_OUTPUT= (OUTER_ITER, INNER_ITER, RMS_PRESSURE, RMS_TEMPERATURE, LIFT, DRAG) + +% ----------------------- GEOMETRY EVALUATION PARAMETERS ----------------------% +GEO_BOUNDS= ( 0.499, 0.501) diff --git a/TestCases/py_wrapper/custom_source/run.py b/TestCases/py_wrapper/custom_source/run.py new file mode 100644 index 00000000000..ea08aaf98bc --- /dev/null +++ b/TestCases/py_wrapper/custom_source/run.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Buoyancy force using user defines source term +# \version 8.1.0 "Harrier" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import sys +import pysu2 +# from mpi4py import MPI + +def main(): + """ + Run the flow solver with a custom inlet (function of time and space). + """ + # comm = MPI.COMM_WORLD + comm = 0 + + # Initialize the primal driver of SU2, this includes solver preprocessing. + try: + driver = pysu2.CSinglezoneDriver('lam_buoyancy_cavity.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CSinglezoneDriver : ', exception) + raise + + print("\n------------------------------ Begin Solver -----------------------------") + sys.stdout.flush() + + # we need to add a source term to the energy equation. For this, we need to get the solver, and the variable first. + # we then loop over all points and for these points, we add the source term + nDim = driver.GetNumberDimensions() + + # index to the flow solver + iSOLVER = driver.GetSolverIndices()['INC.FLOW'] + print("index of flow solver = ",iSOLVER) + + # all the indices and the map to the names of the primitives + primindex = driver.GetPrimitiveIndices() + print("indices of primitives=",primindex) + print("number of primitives:",len(primindex)) + + nElem = driver.GetNumberElements() + print("number of elements:",nElem) + + nVars = driver.GetNumberSolverVars(iSOLVER) + print("number of solver variables:",nVars) + varindex = primindex.copy() + for prim in varindex.copy(): + if varindex[prim] >=nVars: + del varindex[prim] + varindex = dict(sorted(varindex.items(), key=lambda item: item[1])) + + + + print("solver variable names:",varindex) + iDENSITY = primindex.get("DENSITY") + print("index of density = ",iDENSITY) + + index_Vel = varindex.get("VELOCITY_X") + print("index of velocity = ",index_Vel) + custom_source_vector = [0.0 for i in range(nVars)] + print("custom source vector = ", custom_source_vector) + + print("max. number of inner iterations: ",driver.GetNumberInnerIter()); + print("max nr of outer iterations: ",driver.GetNumberOuterIter()); + + # is in domain: isdomain = driver.GetNodeDomain(iPoint) + #for i_vertex in range(n_vertex) + #AllSolutions = driver.GetAllSolutions(iSOLVER) + Body_Force_Vector = [0.0, -9.81, 0.0] + DensityInc_0 = driver.GetDensity_FreeStreamND() + print("rho freestream = ",DensityInc_0) + Force_Ref = driver.GetForce_Ref() + print("reference force = ",Force_Ref) + + for inner_iter in range(500): + + # set the source term, per point + for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()): + #print("node = ",i_node) + #SolutionVector = driver.GetSolutionVector(iSOLVER,i_node) + #print("solutionvector=",SolutionVector) + PrimitiveVector = driver.GetPrimitiveVector(iSOLVER,i_node) + #print("primitivevector=",PrimitiveVector) + DensityInc_i = PrimitiveVector[iDENSITY] + + for iDim in range(nDim): + custom_source_vector[iDim+1] = -(DensityInc_i - DensityInc_0) * Body_Force_Vector[iDim] / Force_Ref + + #print("density=",DensityInc_i) + driver.SetPointCustomSource(iSOLVER, i_node,custom_source_vector) + + print(" *** inner iteration:",inner_iter) + driver.Preprocess(inner_iter) + + # Run one time iteration. + driver.Run() + driver.Postprocess() + driver.Update() + + # Monitor the solver and output solution to file if required. + #driver.Monitor(inner_iter) + driver.Output(inner_iter) + + # Finalize the solver and exit cleanly. + driver.Finalize() + +if __name__ == '__main__': + main() diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg b/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg new file mode 100644 index 00000000000..ee0cfb797a7 --- /dev/null +++ b/TestCases/py_wrapper/turbulent_premixed_psi/psi.cfg @@ -0,0 +1,167 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Species mixing with 3 species, i.e. 2 transport equations % +% Author: T. Kattmann % +% Institution: Bosch Thermotechniek B.V. % +% Date: 2021/10/14 % +% File Version 7.5.1 "Blackbird" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_RANS +KIND_TURB_MODEL= SST +SST_OPTIONS= V1994m + +MGLEVEL= 0 +% + +% temperature and density +FREESTREAM_TEMPERATURE = 673 +FREESTREAM_DENSITY = 2.55 +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +%INC_DENSITY_MODEL= CONSTANT +INC_DENSITY_MODEL= VARIABLE +INC_DENSITY_INIT= 2.55 +% +INC_VELOCITY_INIT= (40.00, 0.0, 0.0 ) +% +INC_ENERGY_EQUATION= YES +INC_TEMPERATURE_INIT= 673.0 +% +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID PROPERTIES ------------------------------------- % +% +%FLUID_MODEL= CONSTANT_DENSITY +FLUID_MODEL= INC_IDEAL_GAS +% +CONDUCTIVITY_MODEL= CONSTANT_CONDUCTIVITY +THERMAL_CONDUCTIVITY_CONSTANT= 0.0357 +% +PRANDTL_LAM= 0.72 +TURBULENT_CONDUCTIVITY_MODEL= NONE +PRANDTL_TURB= 0.90 +% +VISCOSITY_MODEL= SUTHERLAND +MU_CONSTANT= 1.716E-5 +MU_REF = 1.716e-5 +MU_T_REF= 273.15 +SUTHERLAND_CONSTANT = 110.4 + +SPECIFIC_HEAT_CP = 1150 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( wall_top, 0.0,wall_side,0.0,wall_pipe,0.0, wall_out, 0.0 ) + +% note, case is axisymmetric +MARKER_SYM= ( symmetry ) +AXISYMMETRIC= YES +% +SPECIFIED_INLET_PROFILE= YES +INLET_FILENAME= inlet.dat +%INLET_INTERPOLATION_FUNCTION= LINEAR_1D +INC_INLET_TYPE= VELOCITY_INLET +INC_INLET_DAMPING= 0.01 +MARKER_INLET= ( inlet, 673, 40.0, 1.0, 0.0, 0.0) +MARKER_INLET_TURBULENT = (inlet, 0.10, 15) +MARKER_INLET_SPECIES= (inlet, 0.0) +% +INC_OUTLET_TYPE= PRESSURE_OUTLET +INC_OUTLET_DAMPING= 0.01 +MARKER_OUTLET= ( outlet, 0.0 ) +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +% +%CFL_NUMBER= 25.0 +CFL_NUMBER= 20.0 +CFL_REDUCTION_SPECIES= 1.0 +CFL_REDUCTION_TURB= 1.0 +% +% Run commented Iter for good results +ITER= 1 +% +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-8 +LINEAR_SOLVER_ITER= 5 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= NO +SLOPE_LIMITER_FLOW = NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= SPECIES_TRANSPORT +%DIFFUSIVITY_MODEL= CONSTANT_SCHMIDT +DIFFUSIVITY_MODEL= CONSTANT_DIFFUSIVITY +SCHMIDT_NUMBER_LAMINAR= 1.0 +DIFFUSIVITY_CONSTANT= 7.56e-5 + +% according to the paper +SCHMIDT_NUMBER_TURBULENT= 0.7 + +% +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +MUSCL_SPECIES= NO +SLOPE_LIMITER_SPECIES = NONE +% +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% +SPECIES_INIT= 0.0 +SPECIES_CLIPPING= YES +SPECIES_CLIPPING_MIN= 0.0 +SPECIES_CLIPPING_MAX= 1.0 +% +% -------------------- TURBULENT TRANSPORT ---------------------------------------% +% +CONV_NUM_METHOD_TURB= BOUNDED_SCALAR +MUSCL_TURB= NO + +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_PRESSURE, RMS_VELOCITY-X, RMS_VELOCITY-Y, RMS_TKE, RMS_SPECIES +CONV_RESIDUAL_MINVAL= -18 +CONV_STARTITER= 10 +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +%MESH_FILENAME= psi_coarse.su2 +MESH_FILENAME= psi_medium.su2 +%MESH_FILENAME= psi_fine.su2 +% +SCREEN_OUTPUT= INNER_ITER WALL_TIME \ + RMS_PRESSURE RMS_VELOCITY-X RMS_VELOCITY-Y RMS_TKE RMS_DISSIPATION RMS_SPECIES_0 RMS_SPECIES_1 \ + LINSOL_ITER LINSOL_RESIDUAL \ + LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB \ + LINSOL_ITER_SPECIES LINSOL_RESIDUAL_SPECIES \ + SURFACE_SPECIES_0 +SCREEN_WRT_FREQ_INNER= 1 +% +HISTORY_OUTPUT= ITER RMS_RES LINSOL SPECIES_COEFF SPECIES_COEFF_SURF +CONV_FILENAME= history +MARKER_ANALYZE= gas_inlet, air_axial_inlet, outlet +MARKER_ANALYZE_AVERAGE= AREA +% +OUTPUT_FILES= RESTART_ASCII, PARAVIEW_MULTIBLOCK +VOLUME_OUTPUT= RESIDUAL, PRIMITIVE +OUTPUT_WRT_FREQ= 100 +% +RESTART_SOL= YES +READ_BINARY_RESTART= NO +RESTART_FILENAME= restart +SOLUTION_FILENAME= solution +% +WRT_PERFORMANCE= YES diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/run.py b/TestCases/py_wrapper/turbulent_premixed_psi/run.py new file mode 100644 index 00000000000..6e22c4c152c --- /dev/null +++ b/TestCases/py_wrapper/turbulent_premixed_psi/run.py @@ -0,0 +1,304 @@ +#!/usr/bin/env python + +## \file run.py +# \brief turbulent premixed dump combustor simulation (PSI flame) +# phi=0.5, methane-air, U=40 m/s +# \version 8.1.0 "Harrier" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import sys +import pysu2 +import numpy as np +#from mpi4py import MPI + +# unburnt temperature of the propane-air mixture +# flame temperature of the methane-air mixture (phi=0.5, P=5) +Tf = 1777 + +Tu = 673.0 +Pu = 5.0 +phi = 0.5 +# unburnt density at P=5 +rho_u = 2.52 +# unburnt thermal conductivity of methane-air at phi=0.5 (P=5) +k_u = 0.0523 +# unburnt heat capacity of methane-air at phi=0.5 (P=1) +cp_u = 1311.0 + +# P = rho*R*T +# 5 = 2.55 * R * 673 +# R = 0.0029 + +# ################################################################## # +# create a function for the initial progress variable # +# ################################################################## # +def initC(coord): + x = coord[0] + y = coord[1] + #z = coord[2] + #print("x,y = ",x," ",y) + C = 0.0 + # location where the flame should be + flame_x = 0.012 + if (x < flame_x): + C = 0.0 + else: + C = 1.0 + + return C + +# ################################################################## # +# loop over all vertices and set the species progress variable # +# ################################################################## # +def SetInitialSpecies(SU2Driver): + allCoords = SU2Driver.Coordinates() + iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES'] + print("index of species solver = ",iSPECIESSOLVER) + nVarsSpecies = SU2Driver.GetNumberSolverVars(iSPECIESSOLVER) + print("number of species solver variables:",nVarsSpecies) + for iPoint in range(SU2Driver.GetNumberNodes() - SU2Driver.GetNumberHaloNodes()): + coord = allCoords.Get(iPoint) + C = initC(coord) + # now update the initial condition + SU2Driver.SetSolutionVector(iSPECIESSOLVER, iPoint, [C]) + +def SetInitialVelocity(SU2Driver): + allCoords = SU2Driver.Coordinates() + iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] + print("index of FLOW solver = ",iFLOWSOLVER) + nVarsFlow = SU2Driver.GetNumberSolverVars(iFLOWSOLVER) + print("number of flow solver variables:",nVarsFlow) + for iPoint in range(SU2Driver.GetNumberNodes() - SU2Driver.GetNumberHaloNodes()): + coord = allCoords.Get(iPoint) + C = initC(coord) + # now update the initial condition + SU2Driver.SetSolutionVector(iFLOWSOLVER, iPoint, [C]) + +def update_temperature(SU2Driver, iPoint): + # first, get the progress variable + iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES'] + # returns a list + C = SU2Driver.GetSolutionVector(iSPECIESSOLVER, iPoint) + T = Tu*(1-C[0]) + Tf*C[0] + P0 = 101325.0 + # kg/kmol.K + R = 8.314 + # kg/kmol + M = 28.5 + RHO = P0/(R*T/M) + iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] + solvar = list(SU2Driver.GetSolutionVector(iFLOWSOLVER, iPoint)) + primvar = list(SU2Driver.GetPrimitiveVector(iFLOWSOLVER, iPoint)) + # the list with names + solindex = getsolvar(SU2Driver) + primindex = SU2Driver.GetPrimitiveIndices() + #print("primindex = ",primindex) + # the actual values + #print("solindex = ",solindex) + #print("primvar = ",primvar) + #print("solvar=",solvar) + iTEMP = solindex.get("TEMPERATURE") + iDENSITY = primindex.get("DENSITY") + #iDIFFUSIVITY = primindex.get("DENSITY") + #print("itemp=",iTEMP) + #print("irho=",iDENSITY) + #print("T=",T) + solvar[iTEMP] = T + # + # + SU2Driver.SetSolutionVector(iFLOWSOLVER, iPoint, solvar) + + # also set the primitive variable list + #primvar[iDENSITY] = 1.225 #RHO + #primvar[iTEMP] = T + #SU2Driver.SetPrimitiveVector(iFLOWSOLVER, iPoint, primvar) + + # how do we get for the scalar solver the diffusion coefficient? + + +def zimont(SU2Driver, iPoint): + + iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] + iSSTSOLVER = SU2Driver.GetSolverIndices()['SST'] + iSPECIESSOLVER = SU2Driver.GetSolverIndices()['SPECIES'] + primindex = SU2Driver.GetPrimitiveIndices() + primvar = list(SU2Driver.GetPrimitiveVector(iFLOWSOLVER, iPoint)) + + iDENSITY = primindex.get("DENSITY") + iMU = primindex.get("LAMINAR_VISCOSITY") + + # laminar burning velocity of methane-air at phi=0.5, P=5 + Slu = 0.232 + + rho = primvar[iDENSITY] + mu = primvar[iMU] + nu=mu/rho + tke, dissipation = SU2Driver.GetSolutionVector(iSSTSOLVER,iPoint) + gradc = SU2Driver.GetGradient(iSPECIESSOLVER,iPoint,0) + norm_gradc = np.sqrt(gradc[0]*gradc[0] + gradc[1]*gradc[1]) + + up = np.sqrt((2.0/3.0) * tke ) + lt = (0.09**0.75) * (tke**1.5) / dissipation + Re = up*lt/nu + Le = 1.0 + Ut = Slu * (1.0 + (0.46/Le)*np.power(Re,0.25)*np.power(up/Slu,0.3)*np.power(Pu,0.2)) + + norm_gradc = np.sqrt(gradc[0]*gradc[0] + gradc[1]*gradc[1]) + + Sc = rho_u * Ut * norm_gradc + + return Sc + +def getsolvar(SU2Driver): + primindex = SU2Driver.GetPrimitiveIndices() + iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW'] + nVars = SU2Driver.GetNumberSolverVars(iFLOWSOLVER) + varindex = primindex.copy() + for prim in varindex.copy(): + if varindex[prim] >=nVars: + del varindex[prim] + varindex = dict(sorted(varindex.items(), key=lambda item: item[1])) + return varindex + + +def main(): + """ + Run the flow solver with a custom inlet (function of time and space). + """ + # comm = MPI.COMM_WORLD + comm = 0 + + # Initialize the primal driver of SU2, this includes solver preprocessing. + try: + driver = pysu2.CSinglezoneDriver('psi.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CSinglezoneDriver : ', exception) + raise + + print("\n------------------------------ Begin Solver -----------------------------") + sys.stdout.flush() + + nDim = driver.GetNumberDimensions() + + # index to the flow solver + # C.FLOW + # INC.FLOW + # HEAT + # FLAMELET + # SPECIES + # SA + # SST + iFLOWSOLVER = driver.GetSolverIndices()['INC.FLOW'] + print("index of flow solver = ",iFLOWSOLVER) + iSPECIESSOLVER = driver.GetSolverIndices()['SPECIES'] + print("index of species solver = ",iSPECIESSOLVER) + iSSTSOLVER = driver.GetSolverIndices()['SST'] + print("index of turbulence solver = ",iSSTSOLVER) + + + # all the indices and the map to the names of the primitives + primindex = driver.GetPrimitiveIndices() + print("indices of primitives=",primindex) + print("number of primitives:",len(primindex)) + + nElem = driver.GetNumberElements() + print("number of elements:",nElem) + + nVars = driver.GetNumberSolverVars(iFLOWSOLVER) + print("number of flow solver variables:",nVars) + + nVarsSpecies = driver.GetNumberSolverVars(iSPECIESSOLVER) + print("number of species solver variables:",nVarsSpecies) + nVarsTurb = driver.GetNumberSolverVars(iSSTSOLVER) + print("number of turbulence solver variables:",nVarsTurb) + + varindex = primindex.copy() + for prim in varindex.copy(): + if varindex[prim] >=nVars: + del varindex[prim] + varindex = dict(sorted(varindex.items(), key=lambda item: item[1])) + + +# it is possible to get the solver type by doing +# iFLOWSOLVER = driver.GetSolverIndices()['INC.FLOW'] +# with the solver type we then get the solver variables using: + +# nVars = driver.GetNumberSolverVars(iFLOWSOLVER) +# and the solver variable names: +# print("solver variable names:",varindex) + +# we can overwrite the solution using: +# driver.SetSolutionVector(iSolver, iPoint, solutionVector) +# + + print("solver variable names:",varindex) + iDENSITY = primindex.get("DENSITY") + print("index of density = ",iDENSITY) + + index_Vel = varindex.get("VELOCITY_X") + print("index of velocity = ",index_Vel) + custom_source_vector = [0.0 for i in range(nVars)] + print("custom source vector = ", custom_source_vector) + + print("max. number of inner iterations: ",driver.GetNumberInnerIter()); + print("max nr of outer iterations: ",driver.GetNumberOuterIter()); + + # set initial condition + print("Start calling SetInitialSpecies") + #SetInitialSpecies(driver) + print("End calling SetInitialSpecies") + + # super important to actually push the commands. + sys.stdout.flush() + + # run 500 iterations + for inner_iter in range(500): + + driver.Preprocess(inner_iter) + driver.Run() + + # set the source term, per point + for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()): + #print("inode=",i_node) + # add source term: + # default TFC of Zimont: rho*Sc = rho_u * U_t * grad(c) + S = [zimont(driver,i_node)] + driver.SetPointCustomSource(iSPECIESSOLVER, i_node,S) + # at this point we also need to update the temperature based on the progress variable: + # sete the temperature to T = c*Tf + (1-c)*Tu + update_temperature(driver, i_node) + + + driver.Postprocess() + driver.Update() + + + + driver.Monitor(inner_iter) + + driver.Output(inner_iter) + + # Finalize the solver and exit cleanly. + driver.Finalize() + +if __name__ == '__main__': + main()