From cbe714304de08e2b3f30ffb118e7687b7ae6e311 Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Thu, 8 Oct 2020 13:05:01 -0700 Subject: [PATCH 01/12] implementing viscous axisymmetric terms for flow solver: runs but diverges --- Common/include/option_structure.hpp | 37 +++++++------ SU2_CFD/include/numerics/CNumerics.hpp | 14 +++++ SU2_CFD/include/solvers/CSolver.hpp | 22 +++++++- SU2_CFD/include/variables/CEulerVariable.hpp | 8 ++- SU2_CFD/include/variables/CVariable.hpp | 58 +++++++++++++++++++- SU2_CFD/src/numerics/flow/flow_sources.cpp | 40 +++++++++++++- SU2_CFD/src/solvers/CEulerSolver.cpp | 53 +++++++++++++++++- SU2_CFD/src/solvers/CSolver.cpp | 40 +++++++++++++- SU2_CFD/src/variables/CEulerVariable.cpp | 7 ++- 9 files changed, 252 insertions(+), 27 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1c9bb1d0e62f..11e892392463 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file option_structure.hpp * \brief Defines classes for referencing options for easy input in CConfig * \author J. Hicken, B. Tracey @@ -2165,23 +2165,24 @@ enum MPI_QUANTITIES { MAX_EIGENVALUE = 11, /*!< \brief Maximum eigenvalue communication. */ SENSOR = 12, /*!< \brief Dissipation sensor communication. */ AUXVAR_GRADIENT = 13, /*!< \brief Auxiliary variable gradient communication. */ - COORDINATES = 14, /*!< \brief Vertex coordinates communication. */ - COORDINATES_OLD = 15, /*!< \brief Old vertex coordinates communication. */ - MAX_LENGTH = 16, /*!< \brief Maximum length communication. */ - GRID_VELOCITY = 17, /*!< \brief Grid velocity communication. */ - CROSS_TERM = 18, /*!< \brief Cross term communication. */ - CROSS_TERM_GEOMETRY = 19, /*!< \brief Geometric cross term communication. */ - REF_GEOMETRY = 20, /*!< \brief Reference geometry communication. */ - SOLUTION_EDDY = 21, /*!< \brief Turbulent solution plus eddy viscosity communication. */ - SOLUTION_MATRIX = 22, /*!< \brief Matrix solution communication. */ - SOLUTION_MATRIXTRANS = 23, /*!< \brief Matrix transposed solution communication. */ - NEIGHBORS = 24, /*!< \brief Neighbor point count communication (for JST). */ - SOLUTION_FEA = 25, /*!< \brief FEA solution communication. */ - SOLUTION_FEA_OLD = 26, /*!< \brief FEA solution old communication. */ - MESH_DISPLACEMENTS = 27, /*!< \brief Mesh displacements at the interface. */ - SOLUTION_TIME_N = 28, /*!< \brief Solution at time n. */ - SOLUTION_TIME_N1 = 29, /*!< \brief Solution at time n-1. */ - PRIMITIVE = 30 /*!< \brief Primitive solution communication. */ + AXIAUXVAR_GRADIENT = 14, /*!< \brief Axisymmetric auxiliary variable gradient communication. */ + COORDINATES = 15, /*!< \brief Vertex coordinates communication. */ + COORDINATES_OLD = 16, /*!< \brief Old vertex coordinates communication. */ + MAX_LENGTH = 17, /*!< \brief Maximum length communication. */ + GRID_VELOCITY = 18, /*!< \brief Grid velocity communication. */ + CROSS_TERM = 19, /*!< \brief Cross term communication. */ + CROSS_TERM_GEOMETRY = 20, /*!< \brief Geometric cross term communication. */ + REF_GEOMETRY = 21, /*!< \brief Reference geometry communication. */ + SOLUTION_EDDY = 22, /*!< \brief Turbulent solution plus eddy viscosity communication. */ + SOLUTION_MATRIX = 23, /*!< \brief Matrix solution communication. */ + SOLUTION_MATRIXTRANS = 24, /*!< \brief Matrix transposed solution communication. */ + NEIGHBORS = 25, /*!< \brief Neighbor point count communication (for JST). */ + SOLUTION_FEA = 26, /*!< \brief FEA solution communication. */ + SOLUTION_FEA_OLD = 27, /*!< \brief FEA solution old communication. */ + MESH_DISPLACEMENTS = 28, /*!< \brief Mesh displacements at the interface. */ + SOLUTION_TIME_N = 29, /*!< \brief Solution at time n. */ + SOLUTION_TIME_N1 = 30, /*!< \brief Solution at time n-1. */ + PRIMITIVE = 31 /*!< \brief Primitive solution communication. */ }; /*! diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 0b7eb53a3111..92f722e8d51b 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -183,6 +183,9 @@ class CNumerics { su2double *AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ *AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ + su2double + **AxiAuxVar_Grad_i, /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */ + **AxiAuxVar_Grad_j; /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ su2double *Coord_i, /*!< \brief Cartesians coordinates of point i. */ @@ -490,6 +493,17 @@ class CNumerics { AuxVar_Grad_j = val_auxvargrad_j; } + /*! + * \brief Set the gradient of the axisymmetric variables. + * \param[in] val_primvar_grad_i - Gradient of the primitive variable at point i. + * \param[in] val_primvar_grad_j - Gradient of the primitive variable at point j. + */ + inline void SetAxiAuxVarGrad(su2double **val_auxvar_grad_i, + su2double **val_auxvar_grad_j) { + AxiAuxVar_Grad_i = val_auxvar_grad_i; + AxiAuxVar_Grad_j = val_auxvar_grad_j; + } + /*! * \brief Set the diffusion coefficient * \param[in] val_diffusioncoeff_i - Value of the diffusion coefficients at i. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index c1c2eaaea026..a6f102c3f3e1 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CSolver.hpp * \brief Headers of the CSolver class which is inherited by all of the other * solvers @@ -592,6 +592,26 @@ class CSolver { */ void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); + /*! + * \brief Compute the Green-Gauss gradient of the auxiliary variable. + * \param[in] geometry - Geometrical definition of the problem. + */ + void SetAxiAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config); + + /*! + * \brief Compute the Least Squares gradient of the auxiliary variable. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetAxiAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config); + + /*! + * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetAxiAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); + /*! * \brief Add External to Solution vector. */ diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 49b9b9b9e564..14da6c5ff2b2 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerVariable.hpp * \brief Class for defining the variables of the compressible Euler solver. * \author F. Palacios, T. Economon @@ -205,6 +205,12 @@ class CEulerVariable : public CVariable { */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } + /*! + * \brief Get the value of the primitive variables gradient. + * \return Value of the primitive variables gradient. + */ + inline su2double **GetAxiAuxVarGradient(unsigned long iPoint) final { return Grad_AxiAuxVar[iPoint]; } + /*! * \brief A virtual member. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 7f4845289b42..0a976609d7e3 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CVariable.hpp * \brief Declaration and inlines of the parent class for defining problem variables, function definitions in file CVariable.cpp. @@ -80,6 +80,9 @@ class CVariable { VectorType AuxVar; /*!< \brief Auxiliar variable for gradient computation. */ MatrixType Grad_AuxVar; /*!< \brief Gradient of the auxiliar variable. */ + MatrixType AxiAuxVar; /*!< \brief Axisymmetric auxiliar variable for gradient computation. */ + CVectorOfMatrix Grad_AxiAuxVar; /*!< \brief Gradient of the axisymmetric variables of the problem. */ + VectorType Max_Lambda_Inv; /*!< \brief Maximun inviscid eingenvalue. */ VectorType Max_Lambda_Visc; /*!< \brief Maximun viscous eingenvalue. */ VectorType Lambda; /*!< \brief Value of the eingenvalue. */ @@ -624,6 +627,53 @@ class CVariable { */ inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iDim) const { return Grad_AuxVar(iPoint,iDim); } + + + + + /*! + * \brief Set axisymmetric auxiliar variables. + * \param[in] iPoint - Point index. + * \param[in] val_auxvar - Value of the auxiliar variable. + */ + inline void SetAxiAuxVar(unsigned long iPoint, const su2double *axiauxvar) { + for (unsigned long iVar = 0; iVar < 3; iVar++) AxiAuxVar(iPoint,iVar) = axiauxvar[iVar]; + } + + /*! + * \brief Get the entire axi aux vector of the problem. + * \return Reference to the axi aux matrix. + */ + inline const MatrixType& GetAxiAuxVar(void) { return AxiAuxVar; } + + /*! + * \brief Set value of axisymmetric auxillary gradients. + * \param[in] iPoint - Point index. + * \param[in] iVar - Index of the variable. + * \param[in] iDim - Index of the dimension. + * \param[in] value - Value of the gradient. + */ + inline void SetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Grad_AxiAuxVar(iPoint,iVar,iDim) = value; } + + /*! + * \brief Get the gradient of the axi auxilary variables. + * \return Reference to gradient. + */ + inline CVectorOfMatrix& GetAxiAuxVarGradient(void) { return Grad_AxiAuxVar; } + + /*! + * \brief Get the value of the axisymmetric auxilliary gradient. + * \param[in] iPoint - Point index. + * \param[in] iVar - Index of the variable. + * \param[in] iDim - Index of the dimension. + * \return Value of the solution gradient. + */ + inline su2double GetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Grad_AxiAuxVar(iPoint,iVar,iDim); } + + + + + /*! * \brief Add a value to the truncation error. * \param[in] iPoint - Point index. @@ -1955,6 +2005,12 @@ class CVariable { */ inline virtual CVectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; } + /*! + * \brief A virtual member. + * \return Value of the primitive variables gradient. + */ + inline virtual su2double **GetAxiAuxVarGradient(unsigned long iPoint) { return nullptr; } + /*! * \brief Set the blending function for the blending of k-w and k-eps. * \param[in] val_viscosity - Value of the vicosity. diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 8d5a5232b197..d9f23c7a198c 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file flow_sources.cpp * \brief Implementation of numerics classes for integration * of source terms in fluid flow problems. @@ -56,9 +56,10 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfig* config) { su2double yinv, Pressure_i, Enthalpy_i, Velocity_i, sq_vel; - unsigned short iDim, iVar, jVar; + unsigned short iDim,jDim, iVar, jVar; bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + bool viscous = (config->GetViscous()); if (Coord_i[1] > EPS) { @@ -105,6 +106,41 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi } + if (viscous){ + + Laminar_Viscosity_i = V_i[nDim+4]; + Eddy_Viscosity_i = V_i[nDim+5]; + Thermal_Conductivity_i = V_i[nDim+6]; + + su2double u = V_i[1]; + su2double v = V_i[2]; + su2double mu = (Laminar_Viscosity_i + + Eddy_Viscosity_i); + + /*--- The full stress tensor is needed for variable density ---*/ + su2double div_vel = 0.0; + for (iDim = 0 ; iDim < nDim; iDim++) + div_vel += PrimVar_Grad_i[iDim+1][iDim]; + + for (iDim = 0 ; iDim < nDim; iDim++) + for (jDim = 0 ; jDim < nDim; jDim++) + tau[iDim][jDim] = (mu*(PrimVar_Grad_i[jDim+1][iDim] + + PrimVar_Grad_i[iDim+1][jDim] ) + -TWO3*mu*div_vel*delta[iDim][jDim]); + + su2double tau_xy = tau[0][1]; + su2double tau_yyp = tau[1][1]; + su2double tau_tt = -TWO3*mu*(div_vel-2*v*yinv); + + su2double qy = -Thermal_Conductivity_i*PrimVar_Grad_i[0][1]; + + residual[0] -= 0.0; + residual[1] -= Volume*(yinv*tau_xy - TWO3*AxiAuxVar_Grad_i[0][0]); + residual[2] -= Volume*(yinv*(tau_yyp-tau_tt-TWO3*mu*v*yinv)-TWO3*AxiAuxVar_Grad_i[0][1]); + residual[3] -= Volume*(yinv*(u*tau_xy+v*tau_yyp-qy-TWO3*mu*v*v*yinv*AxiAuxVar_Grad_i[1][1]) + -AxiAuxVar_Grad_i[2][1]); + + } } else { diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index d3ddb94943b7..40681e5a579e 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerSolver.cpp * \brief Main subrotuines for solving Finite-Volume Euler flow problems. * \author F. Palacios, T. Economon @@ -3039,7 +3039,8 @@ void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh) { - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + const bool viscous = config->GetViscous(); const bool rotating_frame = config->GetRotating_Frame(); const bool axisymmetric = config->GetAxisymmetric(); const bool gravity = (config->GetGravityForce() == YES); @@ -3107,6 +3108,39 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain if (axisymmetric) { + /*--- For viscous problems, we need an additional gradient. ---*/ + if (viscous) { + + for (iPoint = 0; iPoint < nPoint; iPoint++) { + + su2double yCoord = geometry->nodes->GetCoord(iPoint, 1); + su2double yVelocity = nodes->GetVelocity(iPoint,1); + su2double xVelocity = nodes->GetVelocity(iPoint,1); + su2double Total_Viscosity = (nodes->GetLaminarViscosity(iPoint) + + nodes->GetEddyViscosity(iPoint)); + su2double AxiAuxVar[3] = {0.0}; + + if (yCoord > EPS){ + AxiAuxVar[0] = Total_Viscosity*yVelocity/yCoord; + AxiAuxVar[1] = Total_Viscosity*yVelocity*yVelocity/yCoord; + AxiAuxVar[2] = Total_Viscosity*xVelocity*yVelocity/yCoord; + } + + /*--- Set the auxilairy variable for this node. ---*/ + nodes->SetAxiAuxVar(iPoint, AxiAuxVar); + + } + + /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetAxiAuxVar_Gradient_GG(geometry, config); + } + if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetAxiAuxVar_Gradient_LS(geometry, config); + } + + } + /*--- loop over points ---*/ SU2_OMP_FOR_DYN(omp_chunk_size) for (iPoint = 0; iPoint < nPointDomain; iPoint++) { @@ -3120,6 +3154,21 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set y coordinate ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); + /*--- If viscous, we need gradients for extra terms. ---*/ + + if (viscous) { + + /*--- Primitive variables ---*/ + numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr); + + /*--- Gradient of the primitive variables ---*/ + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr); + + /*--- Load the aux variable gradient that we already computed. ---*/ + numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); + + } + /*--- Compute Source term Residual ---*/ auto residual = numerics->ComputeResidual(config); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 3daf4651cf3a..9827c48857aa 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CSolver.cpp * \brief Main subroutines for CSolver class. * \author F. Palacios, T. Economon @@ -1650,6 +1650,10 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; break; + case AXIAUXVAR_GRADIENT: + COUNT_PER_POINT = nDim*3; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; case MESH_DISPLACEMENTS: COUNT_PER_POINT = nDim; MPI_TYPE = COMM_TYPE_DOUBLE; @@ -1781,6 +1785,13 @@ void CSolver::InitiateComms(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) bufDSend[buf_offset+iDim] = base_nodes->GetAuxVarGradient(iPoint, iDim); break; + case AXIAUXVAR_GRADIENT: + for (iVar = 0; iVar < (nDim*3); iVar++){ + for (iDim = 0; iDim < nDim; iDim++){ + bufDSend[buf_offset+iVar*nDim+iDim] = base_nodes->GetAxiAuxVarGradient(iPoint, iVar, iDim); + } + } + break; case SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); @@ -1956,6 +1967,13 @@ void CSolver::CompleteComms(CGeometry *geometry, for (iDim = 0; iDim < nDim; iDim++) base_nodes->SetAuxVarGradient(iPoint, iDim, bufDRecv[buf_offset+iDim]); break; + case AXIAUXVAR_GRADIENT: + for( iVar = 0; iVar < (nDim*3); iVar++ ){ + for (iDim = 0; iDim < nDim; iDim++){ + base_nodes->SetAxiAuxVarGradient(iPoint, iVar, iDim, bufDRecv[buf_offset+iVar*nDim+iDim]); + } + } + break; case SOLUTION_FEA: for (iVar = 0; iVar < nVar; iVar++) { base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); @@ -2514,6 +2532,26 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) weighted, solution, 0, 1, gradient, rmatrix); } +void CSolver::SetAxiAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) { + + const auto& solution = base_nodes->GetAxiAuxVar(); + auto& gradient = base_nodes->GetAxiAuxVarGradient(); + + computeGradientsGreenGauss(this, AXIAUXVAR_GRADIENT, PERIODIC_NONE, *geometry, + *config, solution, 0, 3, gradient); +} + +void CSolver::SetAxiAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) { + + bool weighted = true; + const auto& solution = base_nodes->GetAxiAuxVar(); + auto& gradient = base_nodes->GetAxiAuxVarGradient(); + auto& rmatrix = base_nodes->GetRmatrix(); + + computeGradientsLeastSquares(this, AXIAUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, + weighted, solution, 0, 3, gradient, rmatrix); +} + void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config, bool reconstruction) { const auto& solution = base_nodes->GetSolution(); diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 92361c9388f7..86f88a48ab34 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerVariable.cpp * \brief Definition of the solution fields. * \author F. Palacios, T. Economon @@ -121,6 +121,11 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 Gradient_Aux.resize(nPoint,nPrimVarGrad,nDim,0.0); } + if (config->GetAxisymmetric()){ + Grad_AxiAuxVar.resize(nPoint,3,nDim,0.0); + AxiAuxVar.resize(nPoint,3) = su2double(0.0); + } + if (config->GetLeastSquaresRequired()) { Rmatrix.resize(nPoint,nDim,nDim,0.0); } From aecffd5d555a3213ba9e1e417bc771ba4541015f Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Thu, 8 Oct 2020 15:02:30 -0700 Subject: [PATCH 02/12] fixing indices. axisymmetric is operational explicit. validation work ongoing --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 12 +++++++++--- SU2_CFD/src/solvers/CEulerSolver.cpp | 2 -- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index d9f23c7a198c..91ebc0d82ced 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -108,9 +108,9 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi if (viscous){ - Laminar_Viscosity_i = V_i[nDim+4]; - Eddy_Viscosity_i = V_i[nDim+5]; - Thermal_Conductivity_i = V_i[nDim+6]; + Laminar_Viscosity_i = V_i[nDim+5]; + Eddy_Viscosity_i = V_i[nDim+6]; + Thermal_Conductivity_i = V_i[nDim+7]; su2double u = V_i[1]; su2double v = V_i[2]; @@ -134,6 +134,12 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi su2double qy = -Thermal_Conductivity_i*PrimVar_Grad_i[0][1]; +// cout <<"lam viscosity :"<SetPrimitive(nodes->GetPrimitive(iPoint), nullptr); - /*--- Gradient of the primitive variables ---*/ numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr); From ddedd780376ec233b65f1b4ab0e54d133ea3bb5d Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Mon, 12 Oct 2020 11:15:01 -0700 Subject: [PATCH 03/12] fixing spaces, remove print statements --- SU2_CFD/include/variables/CVariable.hpp | 8 -------- SU2_CFD/src/numerics/flow/flow_sources.cpp | 8 +------- SU2_CFD/src/solvers/CEulerSolver.cpp | 4 +++- 3 files changed, 4 insertions(+), 16 deletions(-) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 0a976609d7e3..17fd64e80e49 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -627,10 +627,6 @@ class CVariable { */ inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iDim) const { return Grad_AuxVar(iPoint,iDim); } - - - - /*! * \brief Set axisymmetric auxiliar variables. * \param[in] iPoint - Point index. @@ -670,10 +666,6 @@ class CVariable { */ inline su2double GetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Grad_AxiAuxVar(iPoint,iVar,iDim); } - - - - /*! * \brief Add a value to the truncation error. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 91ebc0d82ced..d64280d3e415 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -125,7 +125,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi for (iDim = 0 ; iDim < nDim; iDim++) for (jDim = 0 ; jDim < nDim; jDim++) tau[iDim][jDim] = (mu*(PrimVar_Grad_i[jDim+1][iDim] + - PrimVar_Grad_i[iDim+1][jDim] ) + PrimVar_Grad_i[iDim+1][jDim] ) -TWO3*mu*div_vel*delta[iDim][jDim]); su2double tau_xy = tau[0][1]; @@ -134,12 +134,6 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi su2double qy = -Thermal_Conductivity_i*PrimVar_Grad_i[0][1]; -// cout <<"lam viscosity :"<nodes->GetCoord(iPoint, 1); su2double yVelocity = nodes->GetVelocity(iPoint,1); - su2double xVelocity = nodes->GetVelocity(iPoint,1); + su2double xVelocity = nodes->GetVelocity(iPoint,0); su2double Total_Viscosity = (nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint)); su2double AxiAuxVar[3] = {0.0}; @@ -3157,8 +3157,10 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- If viscous, we need gradients for extra terms. ---*/ if (viscous) { + /*--- Primitive variables ---*/ numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr); + /*--- Gradient of the primitive variables ---*/ numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr); From d91643627c3cca26c579dbdfb0067c703213c584 Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Thu, 12 Nov 2020 22:46:37 -0800 Subject: [PATCH 04/12] merge in develop and adding in Florian's work --- Common/include/option_structure.hpp | 3 +- .../include/numerics/flow/flow_sources.hpp | 32 +++- SU2_CFD/include/solvers/CSolver.hpp | 2 +- SU2_CFD/include/variables/CEulerVariable.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 2 +- SU2_CFD/src/drivers/CDriver.cpp | 4 +- SU2_CFD/src/numerics/flow/flow_sources.cpp | 138 ++++++++++++++---- SU2_CFD/src/solvers/CEulerSolver.cpp | 19 ++- SU2_CFD/src/solvers/CSolver.cpp | 2 +- SU2_CFD/src/variables/CEulerVariable.cpp | 2 +- 10 files changed, 159 insertions(+), 47 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 42e8eabd9169..8ed70db6a87d 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file option_structure.hpp * \brief Defines classes for referencing options for easy input in CConfig * \author J. Hicken, B. Tracey @@ -102,6 +102,7 @@ const su2double EPS = 1.0E-16; /*!< \brief Error scale. */ const su2double TURB_EPS = 1.0E-16; /*!< \brief Turbulent Error scale. */ const su2double ONE2 = 0.5; /*!< \brief One divided by two. */ +const su2double ONE3 = 1.0 / 3.0; /*!< \brief One divided by three. */ const su2double TWO3 = 2.0 / 3.0; /*!< \brief Two divided by three. */ const su2double FOUR3 = 4.0 / 3.0; /*!< \brief Four divided by three. */ diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index 6834d8a2a380..ee889780c694 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -63,7 +63,16 @@ class CSourceBase_Flow : public CNumerics { * \ingroup SourceDiscr * \author F. Palacios */ -class CSourceAxisymmetric_Flow final : public CSourceBase_Flow { +class CSourceAxisymmetric_Flow : public CSourceBase_Flow { +protected: + bool implicit, viscous; + su2double yinv{0.0}; + + /*! + * \brief Diffusion residual of the axisymmetric source term. + */ + void ResidualDiffusion(); + public: /*! * \brief Constructor of the class. @@ -74,12 +83,31 @@ class CSourceAxisymmetric_Flow final : public CSourceBase_Flow { CSourceAxisymmetric_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config); /*! - * \brief Residual of the rotational frame source term. + * \brief Residual of the axisymmetric source term. * \param[in] config - Definition of the particular problem. * \return Lightweight const-view of residual and Jacobian. */ ResidualType<> ComputeResidual(const CConfig* config) override; + +}; +/*! + * \class CSourceGeneralAxisymmetric_Flow + * \brief Class for source term for solving axisymmetric problems for a general (non ideal) fluid. + * \ingroup SourceDiscr + * \author F. Dittmann + */ +class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow { +public: + + using CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow; + /*! + * \brief Residual of the general axisymmetric source term. + * \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/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 55d291ce5418..5106ec6de4e1 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CSolver.hpp * \brief Headers of the CSolver class which is inherited by all of the other * solvers diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index f57a05f5e9b8..def82fb09537 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerVariable.hpp * \brief Class for defining the variables of the compressible Euler solver. * \author F. Palacios, T. Economon diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index c90a7b87c401..e41a2971d4fb 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CVariable.hpp * \brief Declaration and inlines of the parent class for defining problem variables, function definitions in file CVariable.cpp. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 6e8530faad5e..16fbd073feb2 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1865,8 +1865,10 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol else if (config->GetAxisymmetric() == YES) { if (incompressible) numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config); - else + else if (ideal_gas) numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config); + else + numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGeneralAxisymmetric_Flow(nDim, nVar_Flow, config); } else if (config->GetGravityForce() == YES) { numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGravity(nDim, nVar_Flow, config); diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 39fda4bff836..23b10a20305b 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -50,16 +50,16 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi Gamma = config->GetGamma(); Gamma_Minus_One = Gamma - 1.0; + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + viscous = config->GetViscous(); } CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfig* config) { - su2double yinv, Pressure_i, Enthalpy_i, Velocity_i, sq_vel; - unsigned short iDim,jDim, iVar, jVar; - - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - bool viscous = (config->GetViscous()); + su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel; + unsigned short iDim, iVar, jVar; if (Coord_i[1] > EPS) { @@ -78,6 +78,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi residual[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0]; residual[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]); residual[3] = yinv*Volume*Enthalpy_i*U_i[2]; + + /*--- Inviscid component of the source term. ---*/ if (implicit) { jacobian[0][0] = 0.0; @@ -105,42 +107,116 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[iVar][jVar] *= yinv*Volume; } + + /*--- Add the viscous terms if necessary. ---*/ - if (viscous){ + if (viscous) ResidualDiffusion(); - Laminar_Viscosity_i = V_i[nDim+5]; - Eddy_Viscosity_i = V_i[nDim+6]; - Thermal_Conductivity_i = V_i[nDim+7]; + } - su2double u = V_i[1]; - su2double v = V_i[2]; - su2double mu = (Laminar_Viscosity_i + - Eddy_Viscosity_i); + else { - /*--- The full stress tensor is needed for variable density ---*/ - su2double div_vel = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - div_vel += PrimVar_Grad_i[iDim+1][iDim]; + for (iVar=0; iVar < nVar; iVar++) + residual[iVar] = 0.0; - for (iDim = 0 ; iDim < nDim; iDim++) - for (jDim = 0 ; jDim < nDim; jDim++) - tau[iDim][jDim] = (mu*(PrimVar_Grad_i[jDim+1][iDim] + - PrimVar_Grad_i[iDim+1][jDim] ) - -TWO3*mu*div_vel*delta[iDim][jDim]); + if (implicit) { + for (iVar=0; iVar < nVar; iVar++) { + for (jVar=0; jVar < nVar; jVar++) + jacobian[iVar][jVar] = 0.0; + } + } - su2double tau_xy = tau[0][1]; - su2double tau_yyp = tau[1][1]; - su2double tau_tt = -TWO3*mu*(div_vel-2*v*yinv); + } + + return ResidualType<>(residual, jacobian, nullptr); +} - su2double qy = -Thermal_Conductivity_i*PrimVar_Grad_i[0][1]; +void CSourceAxisymmetric_Flow::ResidualDiffusion(){ + + su2double laminar_viscosity_i = V_i[nDim+5]; + su2double eddy_viscosity_i = V_i[nDim+6]; + su2double thermal_conductivity_i = V_i[nDim+7]; + su2double heat_capacity_cp_i = V_i[nDim+8]; + + su2double total_viscosity_i = laminar_viscosity_i + eddy_viscosity_i; + su2double total_conductivity_i = thermal_conductivity_i + heat_capacity_cp_i*eddy_viscosity_i/Prandtl_Turb; + + su2double u = U_i[1]/U_i[0]; + su2double v = U_i[2]/U_i[0]; + + residual[0] -= 0.0; + residual[1] -= yinv*Volume*total_viscosity_i*(PrimVar_Grad_i[1][1]+ONE3*PrimVar_Grad_i[2][0]) + - Volume*TWO3*AxiAuxVar_Grad_i[0][0]; // - 2/3 * y * d(v*mu/y)/dx + + residual[2] -= yinv*Volume*total_viscosity_i*FOUR3*(PrimVar_Grad_i[2][1]-v*yinv) + - Volume*TWO3*AxiAuxVar_Grad_i[0][1] ; // - 2/3 * y * d(v*mu/y)/dy + + residual[3] -= yinv*Volume*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]-TWO3*(PrimVar_Grad_i[2][1]-yinv*v)) + - TWO3*v*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[1][0])) + + total_conductivity_i*PrimVar_Grad_i[0][1]) + - Volume*total_viscosity_i*TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1]); // - 2/3 * y *[ d(mu*v*v/y)/dy * d(mu*u*v/y)/dy +} + - residual[0] -= 0.0; - residual[1] -= Volume*(yinv*tau_xy - TWO3*AxiAuxVar_Grad_i[0][0]); - residual[2] -= Volume*(yinv*(tau_yyp-tau_tt-TWO3*mu*v*yinv)-TWO3*AxiAuxVar_Grad_i[0][1]); - residual[3] -= Volume*(yinv*(u*tau_xy+v*tau_yyp-qy-TWO3*mu*v*v*yinv*AxiAuxVar_Grad_i[1][1]) - -AxiAuxVar_Grad_i[2][1]); +CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const CConfig* config) { + unsigned short iVar, jVar; + + if (Coord_i[1] > EPS) { + + yinv = 1.0/Coord_i[1]; + + su2double Density_i = U_i[0]; + su2double Velocity1_i = U_i[1]/U_i[0]; + su2double Velocity2_i = U_i[2]/U_i[0]; + su2double Energy_i = U_i[3]/U_i[0]; + + su2double Pressure_i = V_j[3]; + su2double Enthalpy_i = Energy_i + Pressure_i/Density_i; + + /*--- Inviscid component of the source term. ---*/ + + residual[0] = yinv*Volume*U_i[2]; + residual[1] = yinv*Volume*U_i[1]*Velocity2_i; + residual[2] = yinv*Volume*U_i[2]*Velocity2_i; + residual[3] = yinv*Volume*U_i[2]*Enthalpy_i; + + if (implicit) { + + su2double dPdrho_e_i = S_i[0]; + su2double dPde_rho_i = S_i[1]; + + jacobian[0][0] = 0.0; + jacobian[0][1] = 0.0; + jacobian[0][2] = 1.0; + jacobian[0][3] = 0.0; + + jacobian[1][0] = -Velocity1_i*Velocity2_i; + jacobian[1][1] = Velocity2_i; + jacobian[1][2] = Velocity1_i; + jacobian[1][3] = 0.0; + + jacobian[2][0] = -Velocity2_i*Velocity2_i; + jacobian[2][1] = 0.0; + jacobian[2][2] = 2*Velocity2_i; + jacobian[2][3] = 0.0; + + jacobian[3][0] = Velocity2_i*(dPdrho_e_i + dPde_rho_i/Density_i*(Velocity1_i*Velocity1_i + + Velocity2_i*Velocity2_i + - Energy_i) - Enthalpy_i); + jacobian[3][1] = -Velocity1_i*Velocity2_i/Density_i *dPde_rho_i; + jacobian[3][2] = Enthalpy_i - Velocity2_i*Velocity2_i/Density_i *dPde_rho_i; + jacobian[3][3] = Velocity2_i + Velocity2_i/Density_i *dPde_rho_i; + + for (iVar=0; iVar < nVar; iVar++) + for (jVar=0; jVar < nVar; jVar++) + jacobian[iVar][jVar] *= yinv*Volume; } + + /*--- Add the viscous terms if necessary. ---*/ + + if (viscous) ResidualDiffusion(); + } else { diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 0aae7ba94e49..3e892d7e1638 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerSolver.cpp * \brief Main subrotuines for solving Finite-Volume Euler flow problems. * \author F. Palacios, T. Economon @@ -3047,6 +3047,8 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain const bool harmonic_balance = (config->GetTime_Marching() == HARMONIC_BALANCE); const bool windgust = config->GetWind_Gust(); const bool body_force = config->GetBody_Force(); + const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || + (config->GetKind_FluidModel() == IDEAL_GAS); /*--- Pick one numerics object per thread. ---*/ CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; @@ -3154,15 +3156,18 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set y coordinate ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); - /*--- If viscous, we need gradients for extra terms. ---*/ + /*--- Set primitive variables ---*/ + numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint)); - if (viscous) { + /*--- Set gradient of primitive variables ---*/ + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); - /*--- Primitive variables ---*/ - numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr); + if (!ideal_gas) { + /*--- Set secondary variables ---*/ + numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint)); + } - /*--- Gradient of the primitive variables ---*/ - numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nullptr); + if (viscous){ /*--- Load the aux variable gradient that we already computed. ---*/ numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index db094dfc6491..450b8f6ef601 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CSolver.cpp * \brief Main subroutines for CSolver class. * \author F. Palacios, T. Economon diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 5bfa0b6c321b..175d1f8d9286 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -1,4 +1,4 @@ -/*! +/*! * \file CEulerVariable.cpp * \brief Definition of the solution fields. * \author F. Palacios, T. Economon From a6a19922ea631af7ff3bf45ea133e2d85c936941 Mon Sep 17 00:00:00 2001 From: FlorianDm Date: Mon, 16 Nov 2020 16:41:46 +0100 Subject: [PATCH 05/12] correct viscous terms --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 19 ++++++++--------- SU2_CFD/src/solvers/CEulerSolver.cpp | 24 +++++++++------------- 2 files changed, 19 insertions(+), 24 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 23b10a20305b..c9c967cea86e 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -145,16 +145,15 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ su2double v = U_i[2]/U_i[0]; residual[0] -= 0.0; - residual[1] -= yinv*Volume*total_viscosity_i*(PrimVar_Grad_i[1][1]+ONE3*PrimVar_Grad_i[2][0]) - - Volume*TWO3*AxiAuxVar_Grad_i[0][0]; // - 2/3 * y * d(v*mu/y)/dx - - residual[2] -= yinv*Volume*total_viscosity_i*FOUR3*(PrimVar_Grad_i[2][1]-v*yinv) - - Volume*TWO3*AxiAuxVar_Grad_i[0][1] ; // - 2/3 * y * d(v*mu/y)/dy - - residual[3] -= yinv*Volume*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]-TWO3*(PrimVar_Grad_i[2][1]-yinv*v)) - - TWO3*v*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[1][0])) - + total_conductivity_i*PrimVar_Grad_i[0][1]) - - Volume*total_viscosity_i*TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1]); // - 2/3 * y *[ d(mu*v*v/y)/dy * d(mu*u*v/y)/dy + residual[1] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][0]) + - TWO3*AxiAuxVar_Grad_i[0][0]); + residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv) + - TWO3*AxiAuxVar_Grad_i[0][1]); + residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) + + v*(FOUR3*PrimVar_Grad_i[1][1]+TWO3*PrimVar_Grad_i[1][0]) + - TWO3*v*v*yinv) + - total_conductivity_i*PrimVar_Grad_i[0][1]) + - TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1])); } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 3e892d7e1638..848648b38189 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3156,24 +3156,20 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set y coordinate ---*/ numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(iPoint)); - /*--- Set primitive variables ---*/ - numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint)); - - /*--- Set gradient of primitive variables ---*/ - numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); - - if (!ideal_gas) { - /*--- Set secondary variables ---*/ - numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint)); - } + /*--- Set primitive variables for viscous terms and/or generalised source ---*/ + if (!ideal_gas || viscous) numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint)); - if (viscous){ + /*--- Set secondary variables for generalised source ---*/ + if (!ideal_gas) numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint)); - /*--- Load the aux variable gradient that we already computed. ---*/ + if (viscous) { + + /*--- Set gradient of primitive variables ---*/ + numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); + + /*--- Set gradient of auxillary variables ---*/ numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); - } - /*--- Compute Source term Residual ---*/ auto residual = numerics->ComputeResidual(config); From ebef3821f01f3165d1f28050d8bf0c6c12799169 Mon Sep 17 00:00:00 2001 From: FlorianDm Date: Mon, 16 Nov 2020 17:21:11 +0100 Subject: [PATCH 06/12] space --- SU2_CFD/src/solvers/CEulerSolver.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 848648b38189..7906f5fb4df5 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3170,6 +3170,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set gradient of auxillary variables ---*/ numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); } + /*--- Compute Source term Residual ---*/ auto residual = numerics->ComputeResidual(config); From b6c3b02f32901d21dbcc9b2da397aac68828875a Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Thu, 19 Nov 2020 21:25:23 -0800 Subject: [PATCH 07/12] addressing PR comments --- SU2_CFD/include/variables/CEulerVariable.hpp | 6 ------ SU2_CFD/include/variables/CVariable.hpp | 6 +++--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index def82fb09537..de454a07ee0a 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -205,12 +205,6 @@ class CEulerVariable : public CVariable { */ inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; } - /*! - * \brief Get the value of the primitive variables gradient. - * \return Value of the primitive variables gradient. - */ - inline su2double **GetAxiAuxVarGradient(unsigned long iPoint) final { return Grad_AxiAuxVar[iPoint]; } - /*! * \brief A virtual member. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index e41a2971d4fb..deddb798d324 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -640,7 +640,7 @@ class CVariable { * \brief Get the entire axi aux vector of the problem. * \return Reference to the axi aux matrix. */ - inline const MatrixType& GetAxiAuxVar(void) { return AxiAuxVar; } + inline const MatrixType& GetAxiAuxVar(void) const { return AxiAuxVar; } /*! * \brief Set value of axisymmetric auxillary gradients. @@ -1998,10 +1998,10 @@ class CVariable { inline virtual CVectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; } /*! - * \brief A virtual member. + * \brief Get the gradient value of the auxillary axisymmetry variables. * \return Value of the primitive variables gradient. */ - inline virtual su2double **GetAxiAuxVarGradient(unsigned long iPoint) { return nullptr; } + inline su2double **GetAxiAuxVarGradient(unsigned long iPoint) { return Grad_AxiAuxVar[iPoint]; } /*! * \brief Set the blending function for the blending of k-w and k-eps. From b64bbd52693bf6b7b8e910b7cca71900fd2dc9ef Mon Sep 17 00:00:00 2001 From: FlorianDm Date: Fri, 20 Nov 2020 10:16:35 +0100 Subject: [PATCH 08/12] fix and simplify source term in energy eq --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c9c967cea86e..c79700e9ff81 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -150,10 +150,9 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv) - TWO3*AxiAuxVar_Grad_i[0][1]); residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) - + v*(FOUR3*PrimVar_Grad_i[1][1]+TWO3*PrimVar_Grad_i[1][0]) - - TWO3*v*v*yinv) - - total_conductivity_i*PrimVar_Grad_i[0][1]) - - TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1])); + +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0]-v*yinv)) + -total_conductivity_i*PrimVar_Grad_i[0][1]) + -TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1])); } From 0d2ded5807f0d37239c000efa7b2a698ca8b71e4 Mon Sep 17 00:00:00 2001 From: FlorianDm Date: Fri, 20 Nov 2020 12:01:41 +0100 Subject: [PATCH 09/12] add contribution of turbulence kinetic energy to axisymmetry energy source --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 7 ++++--- SU2_CFD/src/solvers/CEulerSolver.cpp | 4 ++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c79700e9ff81..590fa300ac91 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -146,11 +146,12 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ residual[0] -= 0.0; residual[1] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][0]) - - TWO3*AxiAuxVar_Grad_i[0][0]); + -TWO3*AxiAuxVar_Grad_i[0][0]); residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv) - - TWO3*AxiAuxVar_Grad_i[0][1]); + -TWO3*AxiAuxVar_Grad_i[0][1]); residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) - +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0]-v*yinv)) + +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0] + -v*yinv+U_i[0]*turb_ke_i)) -total_conductivity_i*PrimVar_Grad_i[0][1]) -TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1])); } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 7906f5fb4df5..4d9bc3b1366c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3169,6 +3169,10 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain /*--- Set gradient of auxillary variables ---*/ numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); + + /*--- Set turbulence kinetic energy ---*/ + CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); + numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); } /*--- Compute Source term Residual ---*/ From ccfbc2890ad76137c20021bcf2a75e3e51064dae Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Fri, 20 Nov 2020 13:34:31 -0800 Subject: [PATCH 10/12] attempting to generalize the axisymm formulation --- Common/include/option_structure.hpp | 35 ++++--- SU2_CFD/include/numerics/CNumerics.hpp | 29 ++---- SU2_CFD/include/solvers/CSolver.hpp | 20 ---- SU2_CFD/include/variables/CVariable.hpp | 106 ++++++++------------ SU2_CFD/src/numerics/flow/flow_sources.cpp | 18 ++-- SU2_CFD/src/solvers/CAdjEulerSolver.cpp | 6 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 16 +-- SU2_CFD/src/solvers/CFEASolver.cpp | 2 +- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 38 +------ SU2_CFD/src/variables/CAdjEulerVariable.cpp | 4 +- SU2_CFD/src/variables/CEulerVariable.cpp | 6 +- SU2_CFD/src/variables/CIncEulerVariable.cpp | 2 +- 13 files changed, 97 insertions(+), 187 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 8ed70db6a87d..ed8bfdd3e46a 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2166,24 +2166,23 @@ enum MPI_QUANTITIES { MAX_EIGENVALUE = 11, /*!< \brief Maximum eigenvalue communication. */ SENSOR = 12, /*!< \brief Dissipation sensor communication. */ AUXVAR_GRADIENT = 13, /*!< \brief Auxiliary variable gradient communication. */ - AXIAUXVAR_GRADIENT = 14, /*!< \brief Axisymmetric auxiliary variable gradient communication. */ - COORDINATES = 15, /*!< \brief Vertex coordinates communication. */ - COORDINATES_OLD = 16, /*!< \brief Old vertex coordinates communication. */ - MAX_LENGTH = 17, /*!< \brief Maximum length communication. */ - GRID_VELOCITY = 18, /*!< \brief Grid velocity communication. */ - CROSS_TERM = 19, /*!< \brief Cross term communication. */ - CROSS_TERM_GEOMETRY = 20, /*!< \brief Geometric cross term communication. */ - REF_GEOMETRY = 21, /*!< \brief Reference geometry communication. */ - SOLUTION_EDDY = 22, /*!< \brief Turbulent solution plus eddy viscosity communication. */ - SOLUTION_MATRIX = 23, /*!< \brief Matrix solution communication. */ - SOLUTION_MATRIXTRANS = 24, /*!< \brief Matrix transposed solution communication. */ - NEIGHBORS = 25, /*!< \brief Neighbor point count communication (for JST). */ - SOLUTION_FEA = 26, /*!< \brief FEA solution communication. */ - SOLUTION_FEA_OLD = 27, /*!< \brief FEA solution old communication. */ - MESH_DISPLACEMENTS = 28, /*!< \brief Mesh displacements at the interface. */ - SOLUTION_TIME_N = 29, /*!< \brief Solution at time n. */ - SOLUTION_TIME_N1 = 30, /*!< \brief Solution at time n-1. */ - PRIMITIVE = 31 /*!< \brief Primitive solution communication. */ + COORDINATES = 14, /*!< \brief Vertex coordinates communication. */ + COORDINATES_OLD = 15, /*!< \brief Old vertex coordinates communication. */ + MAX_LENGTH = 16, /*!< \brief Maximum length communication. */ + GRID_VELOCITY = 17, /*!< \brief Grid velocity communication. */ + CROSS_TERM = 18, /*!< \brief Cross term communication. */ + CROSS_TERM_GEOMETRY = 19, /*!< \brief Geometric cross term communication. */ + REF_GEOMETRY = 20, /*!< \brief Reference geometry communication. */ + SOLUTION_EDDY = 21, /*!< \brief Turbulent solution plus eddy viscosity communication. */ + SOLUTION_MATRIX = 22, /*!< \brief Matrix solution communication. */ + SOLUTION_MATRIXTRANS = 23, /*!< \brief Matrix transposed solution communication. */ + NEIGHBORS = 24, /*!< \brief Neighbor point count communication (for JST). */ + SOLUTION_FEA = 25, /*!< \brief FEA solution communication. */ + SOLUTION_FEA_OLD = 26, /*!< \brief FEA solution old communication. */ + MESH_DISPLACEMENTS = 27, /*!< \brief Mesh displacements at the interface. */ + SOLUTION_TIME_N = 28, /*!< \brief Solution at time n. */ + SOLUTION_TIME_N1 = 29, /*!< \brief Solution at time n-1. */ + PRIMITIVE = 30 /*!< \brief Primitive solution communication. */ }; /*! diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 8d0ddf44aeaf..f4ec332e2c45 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -181,11 +181,8 @@ class CNumerics { **TurbPsi_Grad_i, /*!< \brief Gradient of adjoint turbulent variables at point i. */ **TurbPsi_Grad_j; /*!< \brief Gradient of adjoint turbulent variables at point j. */ su2double - *AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ - *AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ - su2double - **AxiAuxVar_Grad_i, /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */ - **AxiAuxVar_Grad_j; /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */ + **AuxVar_Grad_i, /*!< \brief Gradient of an auxiliary variable at point i. */ + **AuxVar_Grad_j; /*!< \brief Gradient of an auxiliary variable at point i. */ const su2double *RadVar_Source; /*!< \brief Source term from the radiative heat transfer equation. */ su2double *Coord_i, /*!< \brief Cartesians coordinates of point i. */ @@ -485,23 +482,13 @@ class CNumerics { /*! * \brief Set the gradient of the auxiliary variables. - * \param[in] val_auxvargrad_i - Gradient of the auxiliary variable at point i. - * \param[in] val_auxvargrad_j - Gradient of the auxiliary variable at point j. - */ - inline void SetAuxVarGrad(su2double *val_auxvargrad_i, su2double *val_auxvargrad_j) { - AuxVar_Grad_i = val_auxvargrad_i; - AuxVar_Grad_j = val_auxvargrad_j; - } - - /*! - * \brief Set the gradient of the axisymmetric variables. - * \param[in] val_primvar_grad_i - Gradient of the primitive variable at point i. - * \param[in] val_primvar_grad_j - Gradient of the primitive variable at point j. + * \param[in] val_auxvar_grad_i - Gradient of the auxiliary variable at point i. + * \param[in] val_auxvar_grad_j - Gradient of the auxiliary variable at point j. */ - inline void SetAxiAuxVarGrad(su2double **val_auxvar_grad_i, - su2double **val_auxvar_grad_j) { - AxiAuxVar_Grad_i = val_auxvar_grad_i; - AxiAuxVar_Grad_j = val_auxvar_grad_j; + inline void SetAuxVarGrad(su2double **val_auxvar_grad_i, + su2double **val_auxvar_grad_j) { + AuxVar_Grad_i = val_auxvar_grad_i; + AuxVar_Grad_j = val_auxvar_grad_j; } /*! diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 5106ec6de4e1..6576c775a7b2 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -592,26 +592,6 @@ class CSolver { */ void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); - /*! - * \brief Compute the Green-Gauss gradient of the auxiliary variable. - * \param[in] geometry - Geometrical definition of the problem. - */ - void SetAxiAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config); - - /*! - * \brief Compute the Least Squares gradient of the auxiliary variable. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetAxiAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config); - - /*! - * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetAxiAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); - /*! * \brief Add External to Solution vector. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index deddb798d324..e3ab68bdbf76 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -77,11 +77,8 @@ class CVariable { MatrixType Solution_Max; /*!< \brief Max solution for limiter computation. */ MatrixType Solution_Min; /*!< \brief Min solution for limiter computation. */ - VectorType AuxVar; /*!< \brief Auxiliar variable for gradient computation. */ - MatrixType Grad_AuxVar; /*!< \brief Gradient of the auxiliar variable. */ - - MatrixType AxiAuxVar; /*!< \brief Axisymmetric auxiliar variable for gradient computation. */ - CVectorOfMatrix Grad_AxiAuxVar; /*!< \brief Gradient of the axisymmetric variables of the problem. */ + MatrixType AuxVar; /*!< \brief Auxiliar variable for gradient computation. */ + CVectorOfMatrix Grad_AuxVar; /*!< \brief Gradient of the auxilliary variables of the problem. */ VectorType Max_Lambda_Inv; /*!< \brief Maximun inviscid eingenvalue. */ VectorType Max_Lambda_Visc; /*!< \brief Maximun viscous eingenvalue. */ @@ -567,104 +564,83 @@ class CVariable { inline su2double GetLocalCFL(unsigned long iPoint) const { return LocalCFL(iPoint); } /*! - * \brief Set auxiliar variables, we are looking for the gradient of that variable. + * \brief Set auxiliary variables. * \param[in] iPoint - Point index. * \param[in] val_auxvar - Value of the auxiliar variable. */ - inline void SetAuxVar(unsigned long iPoint, su2double val_auxvar) { AuxVar(iPoint) = val_auxvar; } + inline void SetAuxVar(unsigned long iPoint, const su2double *auxvar) { + for (unsigned long iVar = 0; iVar < 3; iVar++) AuxVar(iPoint,iVar) = auxvar[iVar]; + } /*! - * \brief Get the value of the auxiliary variable. + * \brief Set auxiliary variables. * \param[in] iPoint - Point index. - * \return Value of the auxiliary variable. - */ - inline su2double GetAuxVar(unsigned long iPoint) const { return AuxVar(iPoint); } - - /*! - * \brief Get the auxiliary variable. - * \return 2D view of the auxiliary variable. + * \param[in] iVar - Varriable indexs + * \param[in] val_auxvar - Value of the auxiliar variable. */ - inline C2DDummyLastView GetAuxVar(void) const { - return C2DDummyLastView(AuxVar); + inline void SetAuxVar(unsigned long iPoint, unsigned long iVar, const su2double auxvar) { + AuxVar(iPoint,iVar) = auxvar; //TODO AUX, this may be unneccesary } /*! - * \brief Set the value of the auxiliary variable gradient. - * \param[in] iPoint - Point index. - * \param[in] iDim - Index of the dimension. - * \param[in] val_gradient - Value of the gradient for the index iDim. + * \brief Get the entire Aux vector of the problem. + * \return Reference to the aux matrix. */ - inline void SetAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_gradient) { Grad_AuxVar(iPoint,iDim) = val_gradient; } + inline const MatrixType& GetAuxVar(void) const { return AuxVar; } /*! - * \brief Add a value to the auxiliary variable gradient. - * \param[in] iPoint - Point index. - * \param[in] iDim - Index of the dimension. - * \param[in] val_value - Value of the gradient to be added for the index iDim. + * \brief Get the entire Aux value at Point i. + * \return Reference to the aux matrix. */ - inline void AddAuxVarGradient(unsigned long iPoint, unsigned long iDim, su2double val_value) { Grad_AuxVar(iPoint,iDim) += val_value;} + inline const su2double GetAuxVar(unsigned long iPoint) const { return AuxVar(iPoint,0); } //TODO AUX, is this unsafe? /*! - * \brief Get the gradient of the auxiliary variable. + * \brief Set value of auxillary gradients. * \param[in] iPoint - Point index. - * \return Value of the gradient of the auxiliary variable. - */ - inline su2double *GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } - - /*! - * \brief Get the gradient of the auxiliary variable. - * \return 3D view of the gradient of the auxiliary variable. + * \param[in] iVar - Index of the variable. + * \param[in] iDim - Index of the dimension. + * \param[in] value - Value of the gradient. */ - inline C3DDummyMiddleView GetAuxVarGradient() { - return C3DDummyMiddleView(Grad_AuxVar); + inline void SetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { + Grad_AuxVar(iPoint,iVar,iDim) = value; } /*! - * \brief Get the gradient of the auxiliary variable. + * \brief Add a value to the auxiliary variable gradient. * \param[in] iPoint - Point index. * \param[in] iDim - Index of the dimension. - * \return Value of the gradient of the auxiliary variable for the dimension iDim. - */ - inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iDim) const { return Grad_AuxVar(iPoint,iDim); } - - /*! - * \brief Set axisymmetric auxiliar variables. - * \param[in] iPoint - Point index. - * \param[in] val_auxvar - Value of the auxiliar variable. + * \param[in] val_value - Value of the gradient to be added for the index iDim. */ - inline void SetAxiAuxVar(unsigned long iPoint, const su2double *axiauxvar) { - for (unsigned long iVar = 0; iVar < 3; iVar++) AxiAuxVar(iPoint,iVar) = axiauxvar[iVar]; + inline void AddAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double val_value) { + Grad_AuxVar(iPoint,iVar,iDim) += val_value; } /*! - * \brief Get the entire axi aux vector of the problem. - * \return Reference to the axi aux matrix. + * \brief Get the gradient of the auxilary variables. + * \return Reference to gradient. */ - inline const MatrixType& GetAxiAuxVar(void) const { return AxiAuxVar; } + inline CVectorOfMatrix& GetAuxVarGradient(void) { return Grad_AuxVar; } /*! - * \brief Set value of axisymmetric auxillary gradients. + * \brief Get the value of the auxilliary gradient. * \param[in] iPoint - Point index. * \param[in] iVar - Index of the variable. * \param[in] iDim - Index of the dimension. - * \param[in] value - Value of the gradient. - */ - inline void SetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Grad_AxiAuxVar(iPoint,iVar,iDim) = value; } - - /*! - * \brief Get the gradient of the axi auxilary variables. - * \return Reference to gradient. + * \return Value of the solution gradient. */ - inline CVectorOfMatrix& GetAxiAuxVarGradient(void) { return Grad_AxiAuxVar; } + inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { + return Grad_AuxVar(iPoint,iVar,iDim); + } /*! - * \brief Get the value of the axisymmetric auxilliary gradient. + * \brief Get the value of the auxilliary gradient. * \param[in] iPoint - Point index. * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. * \return Value of the solution gradient. */ - inline su2double GetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Grad_AxiAuxVar(iPoint,iVar,iDim); } + inline su2double *GetAuxVarGradient(unsigned long iPoint, unsigned long iVar) { + return Grad_AuxVar(iPoint,iVar); + } /*! * \brief Add a value to the truncation error. @@ -1998,10 +1974,10 @@ class CVariable { inline virtual CVectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; } /*! - * \brief Get the gradient value of the auxillary axisymmetry variables. + * \brief Get the gradient value of the axisymmetry variables. * \return Value of the primitive variables gradient. */ - inline su2double **GetAxiAuxVarGradient(unsigned long iPoint) { return Grad_AxiAuxVar[iPoint]; } + inline su2double **GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } /*! * \brief Set the blending function for the blending of k-w and k-eps. diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 590fa300ac91..fe439ba48a9f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -146,14 +146,14 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ residual[0] -= 0.0; residual[1] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][0]) - -TWO3*AxiAuxVar_Grad_i[0][0]); + -TWO3*AuxVar_Grad_i[0][0]); residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv) - -TWO3*AxiAuxVar_Grad_i[0][1]); + -TWO3*AuxVar_Grad_i[0][1]); residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) - +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0] - -v*yinv+U_i[0]*turb_ke_i)) - -total_conductivity_i*PrimVar_Grad_i[0][1]) - -TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1])); + +v*TWO3*(2*PrimVar_Grad_i[1][1]-PrimVar_Grad_i[1][0] + -v*yinv+U_i[0]*turb_ke_i)) + -total_conductivity_i*PrimVar_Grad_i[0][1]) + -TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][1])); } @@ -327,10 +327,10 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo /*--- Viscous terms. ---*/ residual[0] -= 0.0; - residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0]); + residual[1] -= Volume*(yinv*tau[0][1] - TWO3*AuxVar_Grad_i[0][0]); residual[2] -= Volume*(yinv*2.0*total_viscosity*PrimVar_Grad_i[2][1] - - yinv*yinv*2.0*total_viscosity*Velocity_i[1] - - TWO3*AuxVar_Grad_i[1]); + yinv* yinv*2.0*total_viscosity*Velocity_i[1] - + TWO3*AuxVar_Grad_i[0][1]); residual[3] -= Volume*yinv*Thermal_Conductivity_i*PrimVar_Grad_i[nDim+1][1]; } diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index e94ae5f4c520..4c0ccddcd78a 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -2385,7 +2385,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1]; for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1]; - nodes->SetAuxVar(iPoint,conspsi); + nodes->SetAuxVar(iPoint,0,conspsi); //TODO AUX /*--- Also load the auxiliary variable for first neighbors ---*/ @@ -2396,7 +2396,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver Enthalpy = solver_container[FLOW_SOL]->GetNodes()->GetEnthalpy(Neigh); conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1]; for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1]; - nodes->SetAuxVar(Neigh,conspsi); + nodes->SetAuxVar(Neigh,0,conspsi); //TODO AUX } } } @@ -2423,7 +2423,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver Area = sqrt(Area); PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); - ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint); + ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint,0); //TODO AUX ConsPsi = nodes->GetAuxVar(iPoint); d_press = 0.0; grad_v = 0.0; v_gradconspsi = 0.0; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 4d9bc3b1366c..0940dc3f6a50 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3120,25 +3120,25 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain su2double xVelocity = nodes->GetVelocity(iPoint,0); su2double Total_Viscosity = (nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint)); - su2double AxiAuxVar[3] = {0.0}; + su2double AuxVar[3] = {0.0}; if (yCoord > EPS){ - AxiAuxVar[0] = Total_Viscosity*yVelocity/yCoord; - AxiAuxVar[1] = Total_Viscosity*yVelocity*yVelocity/yCoord; - AxiAuxVar[2] = Total_Viscosity*xVelocity*yVelocity/yCoord; + AuxVar[0] = Total_Viscosity*yVelocity/yCoord; + AuxVar[1] = Total_Viscosity*yVelocity*yVelocity/yCoord; + AuxVar[2] = Total_Viscosity*xVelocity*yVelocity/yCoord; } /*--- Set the auxilairy variable for this node. ---*/ - nodes->SetAxiAuxVar(iPoint, AxiAuxVar); + nodes->SetAuxVar(iPoint, AuxVar); } /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { - SetAxiAuxVar_Gradient_GG(geometry, config); + SetAuxVar_Gradient_GG(geometry, config); } if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { - SetAxiAuxVar_Gradient_LS(geometry, config); + SetAuxVar_Gradient_LS(geometry, config); } } @@ -3168,7 +3168,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); /*--- Set gradient of auxillary variables ---*/ - numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr); + numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); /*--- Set turbulence kinetic energy ---*/ CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 75dcee1a7309..3c9ada1ed2cc 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -3478,7 +3478,7 @@ void CFEASolver::FilterElementDensities(CGeometry *geometry, const CConfig *conf sum += w * element_properties[iElem]->GetPhysicalDensity(); vol += w; } - nodes->SetAuxVar(iPoint, sum/vol); + nodes->SetAuxVar(iPoint, 0, sum/vol); //TODO AUX } } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index e221024c8ba8..d0defb28a87e 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1572,7 +1572,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Set the auxilairy variable for this node. ---*/ - nodes->SetAuxVar(iPoint, AuxVar); + nodes->SetAuxVar(iPoint,0, AuxVar); //TODO AUX } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 450b8f6ef601..07b3e09f677e 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1646,11 +1646,7 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case AUXVAR_GRADIENT: - COUNT_PER_POINT = nDim; - MPI_TYPE = COMM_TYPE_DOUBLE; - break; - case AXIAUXVAR_GRADIENT: + case AUXVAR_GRADIENT: //TODO AUX This need to be made general COUNT_PER_POINT = nDim*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; @@ -1782,13 +1778,9 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = base_nodes->GetLimiter_Primitive(iPoint, iVar); break; case AUXVAR_GRADIENT: - for (iDim = 0; iDim < nDim; iDim++) - bufDSend[buf_offset+iDim] = base_nodes->GetAuxVarGradient(iPoint, iDim); - break; - case AXIAUXVAR_GRADIENT: for (iVar = 0; iVar < (nDim*3); iVar++){ for (iDim = 0; iDim < nDim; iDim++){ - bufDSend[buf_offset+iVar*nDim+iDim] = base_nodes->GetAxiAuxVarGradient(iPoint, iVar, iDim); + bufDSend[buf_offset+iVar*nDim+iDim] = base_nodes->GetAuxVarGradient(iPoint, iVar, iDim); } } break; @@ -1964,13 +1956,9 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetLimiter_Primitive(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; case AUXVAR_GRADIENT: - for (iDim = 0; iDim < nDim; iDim++) - base_nodes->SetAuxVarGradient(iPoint, iDim, bufDRecv[buf_offset+iDim]); - break; - case AXIAUXVAR_GRADIENT: for( iVar = 0; iVar < (nDim*3); iVar++ ){ for (iDim = 0; iDim < nDim; iDim++){ - base_nodes->SetAxiAuxVarGradient(iPoint, iVar, iDim, bufDRecv[buf_offset+iVar*nDim+iDim]); + base_nodes->SetAuxVarGradient(iPoint, iVar, iDim, bufDRecv[buf_offset+iVar*nDim+iDim]); } } break; @@ -2532,26 +2520,6 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) weighted, solution, 0, 1, gradient, rmatrix); } -void CSolver::SetAxiAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) { - - const auto& solution = base_nodes->GetAxiAuxVar(); - auto& gradient = base_nodes->GetAxiAuxVarGradient(); - - computeGradientsGreenGauss(this, AXIAUXVAR_GRADIENT, PERIODIC_NONE, *geometry, - *config, solution, 0, 3, gradient); -} - -void CSolver::SetAxiAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) { - - bool weighted = true; - const auto& solution = base_nodes->GetAxiAuxVar(); - auto& gradient = base_nodes->GetAxiAuxVarGradient(); - auto& rmatrix = base_nodes->GetRmatrix(); - - computeGradientsLeastSquares(this, AXIAUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, - weighted, solution, 0, 3, gradient, rmatrix); -} - void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config, bool reconstruction) { const auto& solution = base_nodes->GetSolution(); diff --git a/SU2_CFD/src/variables/CAdjEulerVariable.cpp b/SU2_CFD/src/variables/CAdjEulerVariable.cpp index b6ad2fc1c11f..08500014e795 100644 --- a/SU2_CFD/src/variables/CAdjEulerVariable.cpp +++ b/SU2_CFD/src/variables/CAdjEulerVariable.cpp @@ -84,8 +84,8 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2 } /*--- Allocate auxiliar vector for sensitivity computation ---*/ - AuxVar.resize(nPoint); - Grad_AuxVar.resize(nPoint,nDim); + AuxVar.resize(1,nPoint); + Grad_AuxVar.resize(nPoint,1,nDim); //TODO make general /*--- Allocate and initializate projection vector for wall boundary condition ---*/ ForceProj_Vector.resize(nPoint,nDim) = su2double(0.0); diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index 175d1f8d9286..d2bc4b40cd0d 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -121,9 +121,9 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 Gradient_Aux.resize(nPoint,nPrimVarGrad,nDim,0.0); } - if (config->GetAxisymmetric()){ - Grad_AxiAuxVar.resize(nPoint,3,nDim,0.0); - AxiAuxVar.resize(nPoint,3) = su2double(0.0); + if (config->GetAxisymmetric()){ //TODO AUX does this need to be general + Grad_AuxVar.resize(nPoint,3,nDim,0.0); + AuxVar.resize(nPoint,3) = su2double(0.0); } if (config->GetLeastSquaresRequired()) { diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index c08f587d592a..4f602212d279 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -107,7 +107,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci /*--- If axisymmetric and viscous, we need an auxiliary gradient. ---*/ - if (axisymmetric && viscous) Grad_AuxVar.resize(nPoint,nDim); + if (axisymmetric && viscous) Grad_AuxVar.resize(nPoint,1,nDim); //TODO AUX make general? if (config->GetMultizone_Problem()) Set_BGSSolution_k(); From e394f9b90acc4a93f085d2f0079705119d98fd5d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 20 Nov 2020 22:55:40 +0000 Subject: [PATCH 11/12] fix compilation issues --- SU2_CFD/include/solvers/CAdjEulerSolver.hpp | 7 ++ SU2_CFD/include/solvers/CSolver.hpp | 7 -- SU2_CFD/include/variables/CVariable.hpp | 59 +++------- SU2_CFD/src/solvers/CAdjEulerSolver.cpp | 118 ++++++++++++++++++- SU2_CFD/src/solvers/CEulerSolver.cpp | 23 ++-- SU2_CFD/src/solvers/CFEASolver.cpp | 2 +- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 +- SU2_CFD/src/solvers/CSolver.cpp | 124 +------------------- SU2_CFD/src/variables/CAdjEulerVariable.cpp | 13 +- SU2_CFD/src/variables/CEulerVariable.cpp | 7 +- SU2_CFD/src/variables/CFEAVariable.cpp | 5 +- SU2_CFD/src/variables/CIncEulerVariable.cpp | 6 - SU2_CFD/src/variables/CIncNSVariable.cpp | 7 ++ 13 files changed, 180 insertions(+), 200 deletions(-) diff --git a/SU2_CFD/include/solvers/CAdjEulerSolver.hpp b/SU2_CFD/include/solvers/CAdjEulerSolver.hpp index ff55d3cb9bad..697ef9e56281 100644 --- a/SU2_CFD/include/solvers/CAdjEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CAdjEulerSolver.hpp @@ -78,6 +78,13 @@ class CAdjEulerSolver : public CSolver { */ inline CVariable* GetBaseClassPointerToNodes() override { return nodes; } + /*! + * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); + public: /*! diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 6576c775a7b2..d2e5de673ef6 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -585,13 +585,6 @@ class CSolver { */ void SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config); - /*! - * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config); - /*! * \brief Add External to Solution vector. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index e3ab68bdbf76..7e905a2630b4 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -105,6 +105,7 @@ class CVariable { unsigned long nPrimVarGrad = 0; /*!< \brief Number of primitives for which a gradient is computed. */ unsigned long nSecondaryVar = 0; /*!< \brief Number of secondary variables. */ unsigned long nSecondaryVarGrad = 0; /*!< \brief Number of secondaries for which a gradient is computed. */ + unsigned long nAuxVar = 0; /*!< \brief Number of auxiliary variables. */ /*--- Only allow default construction by derived classes. ---*/ CVariable() = default; @@ -137,6 +138,11 @@ class CVariable { */ virtual ~CVariable() = default; + /*! + * \brief Get the number of auxiliary variables. + */ + inline unsigned long GetnAuxVar() const { return nAuxVar; } + /*! * \brief Set the value of the solution, all variables. * \param[in] iPoint - Point index. @@ -564,13 +570,15 @@ class CVariable { inline su2double GetLocalCFL(unsigned long iPoint) const { return LocalCFL(iPoint); } /*! - * \brief Set auxiliary variables. - * \param[in] iPoint - Point index. - * \param[in] val_auxvar - Value of the auxiliar variable. + * \brief Get the entire Aux matrix of the problem. + * \return Reference to the aux var matrix. */ - inline void SetAuxVar(unsigned long iPoint, const su2double *auxvar) { - for (unsigned long iVar = 0; iVar < 3; iVar++) AuxVar(iPoint,iVar) = auxvar[iVar]; - } + inline const MatrixType& GetAuxVar(void) const { return AuxVar; } + + /*! + * \brief Get the Aux var value at Point i, variable j. + */ + inline su2double GetAuxVar(unsigned long iPoint, unsigned long iVar = 0) const { return AuxVar(iPoint,iVar); } /*! * \brief Set auxiliary variables. @@ -579,21 +587,9 @@ class CVariable { * \param[in] val_auxvar - Value of the auxiliar variable. */ inline void SetAuxVar(unsigned long iPoint, unsigned long iVar, const su2double auxvar) { - AuxVar(iPoint,iVar) = auxvar; //TODO AUX, this may be unneccesary + AuxVar(iPoint,iVar) = auxvar; } - /*! - * \brief Get the entire Aux vector of the problem. - * \return Reference to the aux matrix. - */ - inline const MatrixType& GetAuxVar(void) const { return AuxVar; } - - /*! - * \brief Get the entire Aux value at Point i. - * \return Reference to the aux matrix. - */ - inline const su2double GetAuxVar(unsigned long iPoint) const { return AuxVar(iPoint,0); } //TODO AUX, is this unsafe? - /*! * \brief Set value of auxillary gradients. * \param[in] iPoint - Point index. @@ -602,17 +598,7 @@ class CVariable { * \param[in] value - Value of the gradient. */ inline void SetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { - Grad_AuxVar(iPoint,iVar,iDim) = value; - } - - /*! - * \brief Add a value to the auxiliary variable gradient. - * \param[in] iPoint - Point index. - * \param[in] iDim - Index of the dimension. - * \param[in] val_value - Value of the gradient to be added for the index iDim. - */ - inline void AddAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double val_value) { - Grad_AuxVar(iPoint,iVar,iDim) += val_value; + Grad_AuxVar(iPoint,iVar,iDim) = value; } /*! @@ -628,18 +614,17 @@ class CVariable { * \param[in] iDim - Index of the dimension. * \return Value of the solution gradient. */ - inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { + inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Grad_AuxVar(iPoint,iVar,iDim); } /*! * \brief Get the value of the auxilliary gradient. * \param[in] iPoint - Point index. - * \param[in] iVar - Index of the variable. * \return Value of the solution gradient. */ - inline su2double *GetAuxVarGradient(unsigned long iPoint, unsigned long iVar) { - return Grad_AuxVar(iPoint,iVar); + inline su2double** GetAuxVarGradient(unsigned long iPoint) { + return Grad_AuxVar[iPoint]; } /*! @@ -1973,12 +1958,6 @@ class CVariable { */ inline virtual CVectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; } - /*! - * \brief Get the gradient value of the axisymmetry variables. - * \return Value of the primitive variables gradient. - */ - inline su2double **GetAuxVarGradient(unsigned long iPoint) { return Grad_AuxVar[iPoint]; } - /*! * \brief Set the blending function for the blending of k-w and k-eps. * \param[in] val_viscosity - Value of the vicosity. diff --git a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp index 4c0ccddcd78a..b5483452b200 100644 --- a/SU2_CFD/src/solvers/CAdjEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CAdjEulerSolver.cpp @@ -2385,7 +2385,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1]; for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1]; - nodes->SetAuxVar(iPoint,0,conspsi); //TODO AUX + nodes->SetAuxVar(iPoint,0,conspsi); /*--- Also load the auxiliary variable for first neighbors ---*/ @@ -2396,7 +2396,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver Enthalpy = solver_container[FLOW_SOL]->GetNodes()->GetEnthalpy(Neigh); conspsi = U[0]*Psi[0] + U[0]*Enthalpy*Psi[nDim+1]; for (iDim = 0; iDim < nDim; iDim++) conspsi += U[iDim+1]*Psi[iDim+1]; - nodes->SetAuxVar(Neigh,0,conspsi); //TODO AUX + nodes->SetAuxVar(Neigh,0,conspsi); } } } @@ -2423,7 +2423,7 @@ void CAdjEulerSolver::Inviscid_Sensitivity(CGeometry *geometry, CSolver **solver Area = sqrt(Area); PrimVar_Grad = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); - ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint,0); //TODO AUX + ConsPsi_Grad = nodes->GetAuxVarGradient(iPoint)[0]; ConsPsi = nodes->GetAuxVar(iPoint); d_press = 0.0; grad_v = 0.0; v_gradconspsi = 0.0; @@ -4884,3 +4884,115 @@ void CAdjEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf Restart_Vars = nullptr; Restart_Data = nullptr; } + +void CAdjEulerSolver::SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config) { + + unsigned short iDim, jDim, iNeigh, iMarker; + unsigned short nDim = geometry->GetnDim(); + unsigned long iPoint, jPoint, iVertex; + su2double *Coord_i, *Coord_j, AuxVar_i, AuxVar_j; + su2double **Smatrix, *Cvector; + + Smatrix = new su2double* [nDim]; + Cvector = new su2double [nDim]; + for (iDim = 0; iDim < nDim; iDim++) + Smatrix[iDim] = new su2double [nDim]; + + + /*--- Loop over boundary markers to select those for Euler or NS walls ---*/ + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + + if (config->GetSolid_Wall(iMarker)) { + + /*--- Loop over points on the surface (Least-Squares approximation) ---*/ + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (geometry->nodes->GetDomain(iPoint)) { + Coord_i = geometry->nodes->GetCoord(iPoint); + AuxVar_i = nodes->GetAuxVar(iPoint); + + /*--- Inizialization of variables ---*/ + for (iDim = 0; iDim < nDim; iDim++) + Cvector[iDim] = 0.0; + su2double r11 = 0.0, r12 = 0.0, r13 = 0.0, r22 = 0.0, r23 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0; + + for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { + jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + Coord_j = geometry->nodes->GetCoord(jPoint); + AuxVar_j = nodes->GetAuxVar(jPoint); + + su2double weight = 0; + for (iDim = 0; iDim < nDim; iDim++) + weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); + + /*--- Sumations for entries of upper triangular matrix R ---*/ + r11 += (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight; + r12 += (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight; + r22 += (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight; + if (nDim == 3) { + r13 += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; + r23_a += (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight; + r23_b += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; + r33 += (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight; + } + + /*--- Entries of c:= transpose(A)*b ---*/ + for (iDim = 0; iDim < nDim; iDim++) + Cvector[iDim] += (Coord_j[iDim]-Coord_i[iDim])*(AuxVar_j-AuxVar_i)/weight; + } + + /*--- Entries of upper triangular matrix R ---*/ + r11 = sqrt(r11); + r12 = r12/r11; + r22 = sqrt(r22-r12*r12); + if (nDim == 3) { + r13 = r13/r11; + r23 = r23_a/r22 - r23_b*r12/(r11*r22); + r33 = sqrt(r33-r23*r23-r13*r13); + } + /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ + if (nDim == 2) { + su2double detR2 = (r11*r22)*(r11*r22); + Smatrix[0][0] = (r12*r12+r22*r22)/detR2; + Smatrix[0][1] = -r11*r12/detR2; + Smatrix[1][0] = Smatrix[0][1]; + Smatrix[1][1] = r11*r11/detR2; + } + else { + su2double detR2 = (r11*r22*r33)*(r11*r22*r33); + su2double z11, z12, z13, z22, z23, z33; // aux vars + z11 = r22*r33; + z12 = -r12*r33; + z13 = r12*r23-r13*r22; + z22 = r11*r33; + z23 = -r11*r23; + z33 = r11*r22; + Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; + Smatrix[0][1] = (z12*z22+z13*z23)/detR2; + Smatrix[0][2] = (z13*z33)/detR2; + Smatrix[1][0] = Smatrix[0][1]; + Smatrix[1][1] = (z22*z22+z23*z23)/detR2; + Smatrix[1][2] = (z23*z33)/detR2; + Smatrix[2][0] = Smatrix[0][2]; + Smatrix[2][1] = Smatrix[1][2]; + Smatrix[2][2] = (z33*z33)/detR2; + } + /*--- Computation of the gradient: S*c ---*/ + su2double product; + for (iDim = 0; iDim < nDim; iDim++) { + product = 0.0; + for (jDim = 0; jDim < nDim; jDim++) + product += Smatrix[iDim][jDim]*Cvector[jDim]; + nodes->SetAuxVarGradient(iPoint, 0, iDim, product); + } + } + } /*--- End of loop over surface points ---*/ + } + } + + /*--- Memory deallocation ---*/ + for (iDim = 0; iDim < nDim; iDim++) + delete [] Smatrix[iDim]; + delete [] Cvector; + delete [] Smatrix; +} diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 0940dc3f6a50..8d746b035f88 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3118,19 +3118,14 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain su2double yCoord = geometry->nodes->GetCoord(iPoint, 1); su2double yVelocity = nodes->GetVelocity(iPoint,1); su2double xVelocity = nodes->GetVelocity(iPoint,0); - su2double Total_Viscosity = (nodes->GetLaminarViscosity(iPoint) + - nodes->GetEddyViscosity(iPoint)); - su2double AuxVar[3] = {0.0}; + su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint); if (yCoord > EPS){ - AuxVar[0] = Total_Viscosity*yVelocity/yCoord; - AuxVar[1] = Total_Viscosity*yVelocity*yVelocity/yCoord; - AuxVar[2] = Total_Viscosity*xVelocity*yVelocity/yCoord; + su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord; + nodes->SetAuxVar(iPoint, 0, nu_v_on_y); + nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity); + nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity); } - - /*--- Set the auxilairy variable for this node. ---*/ - nodes->SetAuxVar(iPoint, AuxVar); - } /*--- Compute the auxiliary variable gradient with GG or WLS. ---*/ @@ -3163,18 +3158,18 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain if (!ideal_gas) numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint)); if (viscous) { - + /*--- Set gradient of primitive variables ---*/ numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint)); - + /*--- Set gradient of auxillary variables ---*/ numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); - + /*--- Set turbulence kinetic energy ---*/ CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); } - + /*--- Compute Source term Residual ---*/ auto residual = numerics->ComputeResidual(config); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 3c9ada1ed2cc..6a19b843c1bb 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -3478,7 +3478,7 @@ void CFEASolver::FilterElementDensities(CGeometry *geometry, const CConfig *conf sum += w * element_properties[iElem]->GetPhysicalDensity(); vol += w; } - nodes->SetAuxVar(iPoint, 0, sum/vol); //TODO AUX + nodes->SetAuxVar(iPoint, 0, sum/vol); } } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index d0defb28a87e..8162857c14d6 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1572,7 +1572,7 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Set the auxilairy variable for this node. ---*/ - nodes->SetAuxVar(iPoint,0, AuxVar); //TODO AUX + nodes->SetAuxVar(iPoint, 0, AuxVar); } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 07b3e09f677e..bf87f5c48506 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1646,8 +1646,8 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case AUXVAR_GRADIENT: //TODO AUX This need to be made general - COUNT_PER_POINT = nDim*3; + case AUXVAR_GRADIENT: + COUNT_PER_POINT = nDim*base_nodes->GetnAuxVar(); MPI_TYPE = COMM_TYPE_DOUBLE; break; case MESH_DISPLACEMENTS: @@ -1778,7 +1778,7 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = base_nodes->GetLimiter_Primitive(iPoint, iVar); break; case AUXVAR_GRADIENT: - for (iVar = 0; iVar < (nDim*3); iVar++){ + for (iVar = 0; iVar < base_nodes->GetnAuxVar(); iVar++){ for (iDim = 0; iDim < nDim; iDim++){ bufDSend[buf_offset+iVar*nDim+iDim] = base_nodes->GetAuxVarGradient(iPoint, iVar, iDim); } @@ -1956,7 +1956,7 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetLimiter_Primitive(iPoint, iVar, bufDRecv[buf_offset+iVar]); break; case AUXVAR_GRADIENT: - for( iVar = 0; iVar < (nDim*3); iVar++ ){ + for (iVar = 0; iVar < base_nodes->GetnAuxVar(); iVar++){ for (iDim = 0; iDim < nDim; iDim++){ base_nodes->SetAuxVarGradient(iPoint, iVar, iDim, bufDRecv[buf_offset+iVar*nDim+iDim]); } @@ -2506,7 +2506,7 @@ void CSolver::SetAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config) auto gradient = base_nodes->GetAuxVarGradient(); computeGradientsGreenGauss(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, - *config, solution, 0, 1, gradient); + *config, solution, 0, base_nodes->GetnAuxVar(), gradient); } void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) { @@ -2517,7 +2517,7 @@ void CSolver::SetAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config) auto& rmatrix = base_nodes->GetRmatrix(); computeGradientsLeastSquares(this, AUXVAR_GRADIENT, PERIODIC_NONE, *geometry, *config, - weighted, solution, 0, 1, gradient, rmatrix); + weighted, solution, 0, base_nodes->GetnAuxVar(), gradient, rmatrix); } void CSolver::SetSolution_Gradient_GG(CGeometry *geometry, const CConfig *config, bool reconstruction) { @@ -2598,118 +2598,6 @@ void CSolver::SetGridVel_Gradient(CGeometry *geometry, const CConfig *config) { true, gridVel, 0, nDim, gridVelGrad, rmatrix); } -void CSolver::SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config) { - - unsigned short iDim, jDim, iNeigh, iMarker; - unsigned short nDim = geometry->GetnDim(); - unsigned long iPoint, jPoint, iVertex; - su2double *Coord_i, *Coord_j, AuxVar_i, AuxVar_j; - su2double **Smatrix, *Cvector; - - Smatrix = new su2double* [nDim]; - Cvector = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Smatrix[iDim] = new su2double [nDim]; - - - /*--- Loop over boundary markers to select those for Euler or NS walls ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - - if (config->GetSolid_Wall(iMarker)) { - - /*--- Loop over points on the surface (Least-Squares approximation) ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - if (geometry->nodes->GetDomain(iPoint)) { - Coord_i = geometry->nodes->GetCoord(iPoint); - AuxVar_i = base_nodes->GetAuxVar(iPoint); - - /*--- Inizialization of variables ---*/ - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iDim] = 0.0; - su2double r11 = 0.0, r12 = 0.0, r13 = 0.0, r22 = 0.0, r23 = 0.0, r23_a = 0.0, r23_b = 0.0, r33 = 0.0; - - for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { - jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); - Coord_j = geometry->nodes->GetCoord(jPoint); - AuxVar_j = base_nodes->GetAuxVar(jPoint); - - su2double weight = 0; - for (iDim = 0; iDim < nDim; iDim++) - weight += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - - /*--- Sumations for entries of upper triangular matrix R ---*/ - r11 += (Coord_j[0]-Coord_i[0])*(Coord_j[0]-Coord_i[0])/weight; - r12 += (Coord_j[0]-Coord_i[0])*(Coord_j[1]-Coord_i[1])/weight; - r22 += (Coord_j[1]-Coord_i[1])*(Coord_j[1]-Coord_i[1])/weight; - if (nDim == 3) { - r13 += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; - r23_a += (Coord_j[1]-Coord_i[1])*(Coord_j[2]-Coord_i[2])/weight; - r23_b += (Coord_j[0]-Coord_i[0])*(Coord_j[2]-Coord_i[2])/weight; - r33 += (Coord_j[2]-Coord_i[2])*(Coord_j[2]-Coord_i[2])/weight; - } - - /*--- Entries of c:= transpose(A)*b ---*/ - for (iDim = 0; iDim < nDim; iDim++) - Cvector[iDim] += (Coord_j[iDim]-Coord_i[iDim])*(AuxVar_j-AuxVar_i)/weight; - } - - /*--- Entries of upper triangular matrix R ---*/ - r11 = sqrt(r11); - r12 = r12/r11; - r22 = sqrt(r22-r12*r12); - if (nDim == 3) { - r13 = r13/r11; - r23 = r23_a/r22 - r23_b*r12/(r11*r22); - r33 = sqrt(r33-r23*r23-r13*r13); - } - /*--- S matrix := inv(R)*traspose(inv(R)) ---*/ - if (nDim == 2) { - su2double detR2 = (r11*r22)*(r11*r22); - Smatrix[0][0] = (r12*r12+r22*r22)/detR2; - Smatrix[0][1] = -r11*r12/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = r11*r11/detR2; - } - else { - su2double detR2 = (r11*r22*r33)*(r11*r22*r33); - su2double z11, z12, z13, z22, z23, z33; // aux vars - z11 = r22*r33; - z12 = -r12*r33; - z13 = r12*r23-r13*r22; - z22 = r11*r33; - z23 = -r11*r23; - z33 = r11*r22; - Smatrix[0][0] = (z11*z11+z12*z12+z13*z13)/detR2; - Smatrix[0][1] = (z12*z22+z13*z23)/detR2; - Smatrix[0][2] = (z13*z33)/detR2; - Smatrix[1][0] = Smatrix[0][1]; - Smatrix[1][1] = (z22*z22+z23*z23)/detR2; - Smatrix[1][2] = (z23*z33)/detR2; - Smatrix[2][0] = Smatrix[0][2]; - Smatrix[2][1] = Smatrix[1][2]; - Smatrix[2][2] = (z33*z33)/detR2; - } - /*--- Computation of the gradient: S*c ---*/ - su2double product; - for (iDim = 0; iDim < nDim; iDim++) { - product = 0.0; - for (jDim = 0; jDim < nDim; jDim++) - product += Smatrix[iDim][jDim]*Cvector[jDim]; - base_nodes->SetAuxVarGradient(iPoint, iDim, product); - } - } - } /*--- End of loop over surface points ---*/ - } - } - - /*--- Memory deallocation ---*/ - for (iDim = 0; iDim < nDim; iDim++) - delete [] Smatrix[iDim]; - delete [] Cvector; - delete [] Smatrix; -} - void CSolver::SetSolution_Limiter(CGeometry *geometry, const CConfig *config) { auto kindLimiter = static_cast(config->GetKind_SlopeLimit()); diff --git a/SU2_CFD/src/variables/CAdjEulerVariable.cpp b/SU2_CFD/src/variables/CAdjEulerVariable.cpp index 08500014e795..8d533bbb16a0 100644 --- a/SU2_CFD/src/variables/CAdjEulerVariable.cpp +++ b/SU2_CFD/src/variables/CAdjEulerVariable.cpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -52,7 +52,7 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2 if (config->GetReconstructionGradientRequired()) { Gradient_Aux.resize(nPoint,nVar,nDim,0.0); } - + if (config->GetLeastSquaresRequired()) { Rmatrix.resize(nPoint,nDim,nDim,0.0); } @@ -84,8 +84,9 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2 } /*--- Allocate auxiliar vector for sensitivity computation ---*/ - AuxVar.resize(1,nPoint); - Grad_AuxVar.resize(nPoint,1,nDim); //TODO make general + nAuxVar = 1; + AuxVar.resize(nPoint, nAuxVar) = su2double(0.0); + Grad_AuxVar.resize(nPoint, nAuxVar, nDim); /*--- Allocate and initializate projection vector for wall boundary condition ---*/ ForceProj_Vector.resize(nPoint,nDim) = su2double(0.0); @@ -101,11 +102,11 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2 Set_BGSSolution_k(); Sensor.resize(nPoint); - + /* Non-physical point (first-order) initialization. */ Non_Physical.resize(nPoint) = false; Non_Physical_Counter.resize(nPoint) = 0; - + } bool CAdjEulerVariable::SetPrimVar(unsigned long iPoint, su2double SharpEdge_Distance, bool check, CConfig *config) { diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index d2bc4b40cd0d..4b0ae95cc4b5 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -121,9 +121,10 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2 Gradient_Aux.resize(nPoint,nPrimVarGrad,nDim,0.0); } - if (config->GetAxisymmetric()){ //TODO AUX does this need to be general - Grad_AuxVar.resize(nPoint,3,nDim,0.0); - AuxVar.resize(nPoint,3) = su2double(0.0); + if (config->GetAxisymmetric()){ + nAuxVar = 3; + Grad_AuxVar.resize(nPoint,nAuxVar,nDim,0.0); + AuxVar.resize(nPoint,nAuxVar) = su2double(0.0); } if (config->GetLeastSquaresRequired()) { diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp index 37741da76bb3..e62f7a31a647 100644 --- a/SU2_CFD/src/variables/CFEAVariable.cpp +++ b/SU2_CFD/src/variables/CFEAVariable.cpp @@ -88,7 +88,10 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig if (multizone) Set_BGSSolution_k(); - if (config->GetTopology_Optimization()) AuxVar.resize(nPoint); + if (config->GetTopology_Optimization()) { + nAuxVar = 1; + AuxVar.resize(nPoint); + } } void CFEAVariable::SetSolution_Vel_time_n() { Solution_Vel_time_n = Solution_Vel; } diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 4f602212d279..9873211ec635 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -34,8 +34,6 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci bool dual_time = (config->GetTime_Marching() == DT_STEPPING_1ST) || (config->GetTime_Marching() == DT_STEPPING_2ND); - bool viscous = config->GetViscous(); - bool axisymmetric = config->GetAxisymmetric(); /*--- Allocate and initialize the primitive variables and gradients ---*/ @@ -105,10 +103,6 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci Rmatrix.resize(nPoint,nDim,nDim,0.0); } - /*--- If axisymmetric and viscous, we need an auxiliary gradient. ---*/ - - if (axisymmetric && viscous) Grad_AuxVar.resize(nPoint,1,nDim); //TODO AUX make general? - if (config->GetMultizone_Problem()) Set_BGSSolution_k(); diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index a2e45e1d9cb6..7d15bef1ea2f 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -31,10 +31,17 @@ CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config) : CIncEulerVariable(pressure, velocity, temperature, npoint, ndim, nvar, config) { + Vorticity.resize(nPoint,3); StrainMag.resize(nPoint); DES_LengthScale.resize(nPoint) = su2double(0.0); Max_Lambda_Visc.resize(nPoint); + + if (config->GetAxisymmetric()) { + nAuxVar = 1; + AuxVar.resize(nPoint,nAuxVar) = su2double(0.0); + Grad_AuxVar.resize(nPoint,nAuxVar,nDim); + } } bool CIncNSVariable::SetVorticity_StrainMag() { From f61da0a22264c419fd23995dce6ec3c0879638e3 Mon Sep 17 00:00:00 2001 From: WallyMaier Date: Tue, 24 Nov 2020 17:47:18 -0800 Subject: [PATCH 12/12] allow laminar simulations --- SU2_CFD/include/numerics/flow/flow_sources.hpp | 2 +- SU2_CFD/src/numerics/flow/flow_sources.cpp | 3 +++ SU2_CFD/src/solvers/CEulerSolver.cpp | 7 +++++-- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index ee889780c694..e6d12207342b 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -65,7 +65,7 @@ class CSourceBase_Flow : public CNumerics { */ class CSourceAxisymmetric_Flow : public CSourceBase_Flow { protected: - bool implicit, viscous; + bool implicit, viscous, rans; su2double yinv{0.0}; /*! diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index fe439ba48a9f..5935423f45c7 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -53,6 +53,7 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); viscous = config->GetViscous(); + rans = (config->GetKind_Turb_Model() != NONE); } @@ -133,6 +134,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi void CSourceAxisymmetric_Flow::ResidualDiffusion(){ + if (!rans){ turb_ke_i = 0.0; } + su2double laminar_viscosity_i = V_i[nDim+5]; su2double eddy_viscosity_i = V_i[nDim+6]; su2double thermal_conductivity_i = V_i[nDim+7]; diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 8d746b035f88..30832f52ed1f 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -3049,6 +3049,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain const bool body_force = config->GetBody_Force(); const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) || (config->GetKind_FluidModel() == IDEAL_GAS); + const bool rans = (config->GetKind_Turb_Model() != NONE); /*--- Pick one numerics object per thread. ---*/ CNumerics* numerics = numerics_container[SOURCE_FIRST_TERM + omp_get_thread_num()*MAX_TERMS]; @@ -3166,8 +3167,10 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); /*--- Set turbulence kinetic energy ---*/ - CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); - numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); + if (rans){ + CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); + numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); + } } /*--- Compute Source term Residual ---*/