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..fb7f73c46a1 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,51 @@ 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]; + solver->SetCustomPointSource(iPoint, values); + } + + /*! + * \brief Get the solution vector in a point for a specific solver + * \param[in] iSolver - Solver index. + * \param[in] iPoint - Point index. + * \param[out] solutionvector - Vector values of the solution. + */ +inline vector 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; +} + + /*! + * \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; +} + /// \} protected: @@ -750,6 +820,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 535335947ca..8d23b402bab 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,27 @@ class CFVMFlowSolverBase : public CSolver { Inlet_FlowDir[val_marker][val_vertex][val_dim] = val_flowdir; } + /*! + * \brief Set a component of the unit vector representing the flow direction at an inlet boundary. + * \param[in] val_marker - Surface marker where the flow direction is set. + * \param[in] val_vertex - Vertex of the marker 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 ---*/ + 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 4797293529f..3f10ef7885e 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -123,6 +123,7 @@ 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); /*--- Force definition and coefficient arrays for all of the markers ---*/ 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 CFVMFlowSolverBaseval_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/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..93aab49d930 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -427,6 +427,7 @@ CNumerics::ResidualType<> CSourceIncBodyForce::ComputeResidual(const CConfig* co residual[nDim+1] = 0.0; + return ResidualType<>(residual, jacobian, nullptr); } @@ -837,3 +838,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..18f60dfbffa 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1813,6 +1813,43 @@ 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++) { + 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/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()