diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 311e01958006..804b37d061d5 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -4299,10 +4299,10 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } } - /*--- Grid motion is not yet supported with the incompressible solver. ---*/ + /*--- Rotating frame is not yet supported with the incompressible solver. ---*/ - if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (GetGrid_Movement())) { - SU2_MPI::Error("Support for grid movement not yet implemented for incompressible flows.", CURRENT_FUNCTION); + if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (Kind_GridMovement == ROTATING_FRAME)) { + SU2_MPI::Error("Support for rotating frame simulation not yet implemented for incompressible flows.", CURRENT_FUNCTION); } /*--- Assert that there are two markers being analyzed if the diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index c25476cb17cd..ef6b44d854ee 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5263,6 +5263,42 @@ class CSourceIncBodyForce : public CNumerics { }; +/*! + * \class CSourceIncRotatingFrame_Flow + * \brief Class for a rotating frame source term. + * \ingroup SourceDiscr + */ +class CSourceIncRotatingFrame_Flow : public CNumerics { + +private: + su2double Omega[3]; /*!< \brief Angular velocity */ + bool implicit; /*!< \brief Implicit calculation. */ + +public: + + /*! + * \brief Constructor of the class. + * \param[in] val_nDim - Number of dimensions of the problem. + * \param[in] val_nVar - Number of variables of the problem. + * \param[in] config - Definition of the particular problem. + */ + CSourceIncRotatingFrame_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config); + + /*! + * \brief Destructor of the class. + */ + ~CSourceIncRotatingFrame_Flow(void); + + /*! + * \brief Residual of the rotational frame source term. + * \param[out] val_residual - Pointer to the total residual. + * \param[out] val_Jacobian_i - Jacobian of the numerical method at node i (implicit computation). + * \param[in] config - Definition of the particular problem. + */ + void ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config); + +}; + /*! * \class CSourceBoussinesq * \brief Class for the source term integration of the Boussinesq approximation for incompressible flow. diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp index bb19df436c76..790eff1baa7f 100644 --- a/SU2_CFD/include/solver_structure.hpp +++ b/SU2_CFD/include/solver_structure.hpp @@ -567,6 +567,7 @@ class CSolver { /*! * \brief Load the geometries at the previous time states n and nM1. * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ virtual void Restart_OldGeometry(CGeometry *geometry, CConfig *config); diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index f7cbb6be899a..1f3b2228484a 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -2205,7 +2205,8 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol else if (incompressible && (config->GetKind_DensityModel() == BOUSSINESQ)) numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceBoussinesq(nDim, nVar_Flow, config); else if (config->GetRotating_Frame() == YES) - numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config); + if (incompressible) numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncRotatingFrame_Flow(nDim, nVar_Flow, config); + else numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config); else if (config->GetAxisymmetric() == YES) if (incompressible) numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config); else numerics[iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config); @@ -3831,7 +3832,7 @@ void CDriver::StartSolver(){ /*--- Perform a dynamic mesh update if required. ---*/ - if (!fem_solver) { + if (!fem_solver && !(config_container[ZONE_0]->GetGrid_Movement() && config_container[ZONE_0]->GetDiscrete_Adjoint())) { DynamicMeshUpdate(ExtIter); } diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index a048a6896a2c..95df77910bcc 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -86,8 +86,15 @@ CUpwFDSInc_Flow::~CUpwFDSInc_Flow(void) { void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { + su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0}; + su2double ProjGridVel = 0.0; + AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); AD::SetPreaccIn(Normal, nDim); + if (dynamic_grid) { + AD::SetPreaccIn(GridVel_i, nDim); + AD::SetPreaccIn(GridVel_j, nDim); + } /*--- Face area (norm or the normal vector) ---*/ @@ -120,6 +127,16 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J ProjVelocity += MeanVelocity[iDim]*Normal[iDim]; } + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (dynamic_grid) { + ProjGridVel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + } + ProjVelocity -= ProjGridVel; + } + /*--- Mean variables at points iPoint and jPoint ---*/ MeanDensity = 0.5*(DensityInc_i + DensityInc_j); @@ -214,6 +231,38 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } } + /*--- Corrections due to grid motion ---*/ + if (dynamic_grid) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = DensityInc_i; U_j[0] = DensityInc_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim]; + } + U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j; + + ProjVelocity = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + + /*--- Residual contributions ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Jacobian contributions ---*/ + /*--- Implicit terms ---*/ + if (implicit) { + for (iDim = 0; iDim < nDim; iDim++){ + val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; + val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + } + val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; + val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; + } + } + } + if (!energy) { val_residual[nDim+1] = 0.0; if (implicit) { @@ -277,6 +326,9 @@ CCentJSTInc_Flow::~CCentJSTInc_Flow(void) { void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { + su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0}; + su2double ProjGridVel = 0.0; + /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -336,7 +388,39 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - + + /*--- Corrections due to grid motion ---*/ + if (dynamic_grid) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = DensityInc_i; U_j[0] = DensityInc_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim]; + } + U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j; + + su2double ProjVelocity = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + + /*--- Residual contributions ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Jacobian contributions ---*/ + /*--- Implicit terms ---*/ + if (implicit) { + for (iDim = 0; iDim < nDim; iDim++){ + val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; + val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + } + val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; + val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; + } + } + } + /*--- Computes differences between Laplacians and conservative variables ---*/ for (iVar = 0; iVar < nVar; iVar++) { @@ -351,6 +435,17 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ /*--- Compute the local spectral radius of the preconditioned system and the stretching factor. ---*/ + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (dynamic_grid) { + ProjGridVel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + } + ProjVelocity_i -= ProjGridVel; + ProjVelocity_j -= ProjGridVel; + } + SoundSpeed_i = sqrt(BetaInc2_i*Area*Area); SoundSpeed_j = sqrt(BetaInc2_j*Area*Area); @@ -396,7 +491,6 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - } CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { @@ -442,6 +536,9 @@ CCentLaxInc_Flow::~CCentLaxInc_Flow(void) { void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { + su2double U_i[5] = {0.0,0.0,0.0,0.0,0.0}, U_j[5] = {0.0,0.0,0.0,0.0,0.0}; + su2double ProjGridVel = 0.0, ProjVelocity = 0.0; + /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -503,7 +600,38 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - + + /*--- Corrections due to grid motion ---*/ + if (dynamic_grid) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = DensityInc_i; U_j[0] = DensityInc_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = DensityInc_i*Velocity_i[iDim]; U_j[iDim+1] = DensityInc_j*Velocity_j[iDim]; + } + U_i[nDim+1] = DensityInc_i*Enthalpy_i; U_j[nDim+1] = DensityInc_j*Enthalpy_j; + + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + + /*--- Residual contributions ---*/ + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Jacobian contributions ---*/ + /*--- Implicit terms ---*/ + if (implicit) { + for (iDim = 0; iDim < nDim; iDim++){ + val_Jacobian_i[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_i; + val_Jacobian_j[iDim+1][iDim+1] -= 0.5*ProjVelocity*DensityInc_j; + } + val_Jacobian_i[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_i*Cp_i; + val_Jacobian_j[nDim+1][nDim+1] -= 0.5*ProjVelocity*DensityInc_j*Cp_j; + } + } + } + /*--- Computes differences btw. conservative variables ---*/ for (iVar = 0; iVar < nVar; iVar++) @@ -519,6 +647,17 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ SoundSpeed_i = sqrt(BetaInc2_i*Area*Area); SoundSpeed_j = sqrt(BetaInc2_j*Area*Area); + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (dynamic_grid) { + ProjGridVel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + } + ProjVelocity_i -= ProjGridVel; + ProjVelocity_j -= ProjGridVel; + } + Local_Lambda_i = fabs(ProjVelocity_i)+SoundSpeed_i; Local_Lambda_j = fabs(ProjVelocity_j)+SoundSpeed_j; @@ -558,7 +697,6 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - } CAvgGradInc_Flow::CAvgGradInc_Flow(unsigned short val_nDim, @@ -852,6 +990,72 @@ void CSourceIncBodyForce::ComputeResidual(su2double *val_residual, CConfig *conf } +CSourceIncRotatingFrame_Flow::CSourceIncRotatingFrame_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { + + implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + Gamma = config->GetGamma(); + Gamma_Minus_One = Gamma - 1.0; + + /*--- Retrieve the angular velocity vector from config. ---*/ + for (unsigned short iDim = 0; iDim < 3; iDim++) + Omega[iDim] = config->GetRotation_Rate(iDim)/config->GetOmega_Ref(); + +} + +CSourceIncRotatingFrame_Flow::~CSourceIncRotatingFrame_Flow(void) { } + +void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { + + unsigned short iDim, iVar, jVar; + su2double Momentum[3] = {0,0,0}, + Velocity_i[3] = {0,0,0}; + + /*--- Primitive variables plus momentum at the node (point i) ---*/ + + DensityInc_i = V_i[nDim+2]; + + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = V_i[iDim+1]; + Momentum[iDim] = DensityInc_i*Velocity_i[iDim]; + } + + /*--- Calculate rotating frame source term residual as ( Omega X Rho-U ) ---*/ + + if (nDim == 2) { + val_residual[0] = 0.0; + val_residual[1] = (Omega[1]*Momentum[2] - Omega[2]*Momentum[1])*Volume; + val_residual[2] = (Omega[2]*Momentum[0] - Omega[0]*Momentum[2])*Volume; + val_residual[3] = 0.0; + } else { + val_residual[0] = 0.0; + val_residual[1] = (Omega[1]*Momentum[2] - Omega[2]*Momentum[1])*Volume; + val_residual[2] = (Omega[2]*Momentum[0] - Omega[0]*Momentum[2])*Volume; + val_residual[3] = (Omega[0]*Momentum[1] - Omega[1]*Momentum[0])*Volume; + val_residual[4] = 0.0; + } + + /*--- Calculate the source term Jacobian ---*/ + + if (implicit) { + for (iVar = 0; iVar < nVar; iVar++) + for (jVar = 0; jVar < nVar; jVar++) + val_Jacobian_i[iVar][jVar] = 0.0; + if (nDim == 2) { + val_Jacobian_i[1][2] = -DensityInc_i*Omega[2]*Volume; + val_Jacobian_i[2][1] = DensityInc_i*Omega[2]*Volume; + } else { + val_Jacobian_i[1][2] = -DensityInc_i*Omega[2]*Volume; + val_Jacobian_i[1][3] = DensityInc_i*Omega[1]*Volume; + val_Jacobian_i[2][1] = DensityInc_i*Omega[2]*Volume; + val_Jacobian_i[2][3] = -DensityInc_i*Omega[0]*Volume; + val_Jacobian_i[3][1] = -DensityInc_i*Omega[1]*Volume; + val_Jacobian_i[3][2] = DensityInc_i*Omega[0]*Volume; + } + } + +} + CSourceBoussinesq::CSourceBoussinesq(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { /*--- Store the pointer to the constant body force vector. ---*/ diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 52cc3ec67a8f..d45d0aaafdce 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -2098,35 +2098,38 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont } if (rotating_frame) { - + /*--- Loop over all points ---*/ - + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Load the conservative variables ---*/ - - numerics->SetConservative(node[iPoint]->GetSolution(), - node[iPoint]->GetSolution()); - + + /*--- Load the primitive variables ---*/ + + numerics->SetPrimitive(node[iPoint]->GetPrimitive(), NULL); + + /*--- Set incompressible density ---*/ + + numerics->SetDensity(node[iPoint]->GetDensity(), 0.0); + /*--- Load the volume of the dual mesh cell ---*/ - + numerics->SetVolume(geometry->node[iPoint]->GetVolume()); - + /*--- Compute the rotating frame source residual ---*/ - + numerics->ComputeResidual(Residual, Jacobian_i, config); - + /*--- Add the source residual to the total ---*/ - + LinSysRes.AddBlock(iPoint, Residual); - + /*--- Add the implicit Jacobian contribution ---*/ - + if (implicit) Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); - + } } - + if (axisymmetric) { /*--- Zero out Jacobian structure ---*/ @@ -5781,9 +5784,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver su2double Volume_nM1, Volume_nP1, TimeStep; su2double *Normal = NULL, *GridVel_i = NULL, *GridVel_j = NULL, Residual_GCL; - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - bool variable_density = (config->GetKind_DensityModel() == VARIABLE); - bool energy = config->GetEnergy_Equation(); + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + bool energy = config->GetEnergy_Equation(); /*--- Store the physical time step ---*/ @@ -5817,7 +5819,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_nP1 = node[iPoint]->GetSolution(); /*--- Access the density and Cp at this node (constant for now). ---*/ - + Density = node[iPoint]->GetDensity(); Cp = node[iPoint]->GetSpecificHeatCp(); @@ -5863,61 +5865,13 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - - unsigned short iDim, jDim; - - su2double BetaInc2, Density, dRhodT, Temperature, Cp; - su2double Velocity[3] = {0.0,0.0,0.0}; - - /*--- Access the primitive variables at this node. ---*/ - - Density = node[iPoint]->GetDensity(); - BetaInc2 = node[iPoint]->GetBetaInc2(); - Cp = node[iPoint]->GetSpecificHeatCp(); - Temperature = node[iPoint]->GetTemperature(); - - for (iDim = 0; iDim < nDim; iDim++) - Velocity[iDim] = node[iPoint]->GetVelocity(iDim); - - /*--- We need the derivative of the equation of state to build the - preconditioning matrix. For now, the only option is the ideal gas - law, but in the future, dRhodT should be in the fluid model. ---*/ - - if (variable_density) { - dRhodT = -Density/Temperature; - } else { - dRhodT = 0.0; - } - - /*--- Calculating the inverse of the preconditioning matrix - that multiplies the time derivative during time integration. ---*/ - - /*--- For implicit calculations, we multiply the preconditioner - by the cell volume over the time step and add to the Jac diagonal. ---*/ - - Jacobian_i[0][0] = 1.0/BetaInc2; - for (iDim = 0; iDim < nDim; iDim++) - Jacobian_i[iDim+1][0] = Velocity[iDim]/BetaInc2; - - if (energy) Jacobian_i[nDim+1][0] = Cp*Temperature/BetaInc2; - else Jacobian_i[nDim+1][0] = 0.0; - - for (jDim = 0; jDim < nDim; jDim++) { - Jacobian_i[0][jDim+1] = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - if (iDim == jDim) Jacobian_i[iDim+1][jDim+1] = Density; - else Jacobian_i[iDim+1][jDim+1] = 0.0; - } - Jacobian_i[nDim+1][jDim+1] = 0.0; + SetPreconditioner(config, iPoint); + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar]; } - - Jacobian_i[0][nDim+1] = dRhodT; - for (iDim = 0; iDim < nDim; iDim++) - Jacobian_i[iDim+1][nDim+1] = Velocity[iDim]*dRhodT; - - if (energy) Jacobian_i[nDim+1][nDim+1] = Cp*(dRhodT*Temperature + Density); - else Jacobian_i[nDim+1][nDim+1] = 1.0; - + } + for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) @@ -5980,10 +5934,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = node[iPoint]->GetSolution_time_n(); /*--- Access the density and Cp at this node (constant for now). ---*/ - + Density = node[iPoint]->GetDensity(); Cp = node[iPoint]->GetSpecificHeatCp(); - + /*--- Compute the conservative variable vector for all time levels. ---*/ U_time_n[0] = Density; @@ -5992,8 +5946,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver } U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - for (iVar = 1; iVar < nVar; iVar++) + for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = U_time_n[iVar]*Residual_GCL; + + if (!energy) Residual[nDim+1] = 0.0; LinSysRes.AddBlock(iPoint, Residual); /*--- Compute the GCL component of the source term for node j ---*/ @@ -6006,8 +5962,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver } U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - for (iVar = 1; iVar < nVar; iVar++) + for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = U_time_n[iVar]*Residual_GCL; + + if (!energy) Residual[nDim+1] = 0.0; LinSysRes.SubtractBlock(jPoint, Residual); } @@ -6056,6 +6014,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = U_time_n[iVar]*Residual_GCL; + + if (!energy) Residual[nDim+1] = 0.0; LinSysRes.AddBlock(iPoint, Residual); } @@ -6127,19 +6087,31 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Store the residual and compute the Jacobian contribution due to the dual time source term. ---*/ - + if (!energy) Residual[nDim+1] = 0.0; LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - for (iVar = 1; iVar < nVar; iVar++) { - if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) - Jacobian_i[iVar][iVar] = Volume_nP1/TimeStep; - if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) - Jacobian_i[iVar][iVar] = (3.0*Volume_nP1)/(2.0*TimeStep); + SetPreconditioner(config, iPoint); + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar]; + } + } + + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) + Jacobian_i[iVar][jVar] *= Volume_nP1 / TimeStep; + if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) + Jacobian_i[iVar][jVar] *= (Volume_nP1*3.0)/(2.0*TimeStep); + } + } + + if (!energy) { + for (iVar = 0; iVar < nVar; iVar++) { + Jacobian_i[iVar][nDim+1] = 0.0; + Jacobian_i[nDim+1][iVar] = 0.0; + } } - for (iDim = 0; iDim < nDim; iDim++) - Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1]; - Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1]; - Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } } @@ -6749,7 +6721,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf } /*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy ---*/ - if (dual_time && config->GetGrid_Movement()) + if (dual_time && config->GetGrid_Movement() && (config->GetKind_GridMovement() != RIGID_MOTION)) Restart_OldGeometry(geometry[MESH_0], config); delete [] Coord; diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 711fbe32beca..8610931908fb 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -804,6 +804,18 @@ def main(): ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) + # unsteady pitching NACA0015, SA + unst_inc_turb_naca0015_sa = TestCase('unst_inc_turb_naca0015_sa') + unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" + unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" + unst_inc_turb_naca0015_sa.test_iter = 1 + unst_inc_turb_naca0015_sa.test_vals = [-3.735742, -7.020535, 1.185211, 0.283184] #last 4 columns + unst_inc_turb_naca0015_sa.su2_exec = "parallel_computation.py -f" + unst_inc_turb_naca0015_sa.timeout = 1600 + unst_inc_turb_naca0015_sa.tol = 0.00001 + unst_inc_turb_naca0015_sa.unsteady = True + test_list.append(unst_inc_turb_naca0015_sa) + ###################################### ### NICFD ### ###################################### diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index b6438ea1cdad..dbb6808f52ef 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -803,6 +803,18 @@ def main(): ddes_flatplate.unsteady = True test_list.append(ddes_flatplate) + # unsteady pitching NACA0015, SA + unst_inc_turb_naca0015_sa = TestCase('unst_inc_turb_naca0015_sa') + unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" + unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" + unst_inc_turb_naca0015_sa.test_iter = 1 + unst_inc_turb_naca0015_sa.test_vals = [-3.734989, -7.016510, 1.176112, 0.282917] #last 4 columns + unst_inc_turb_naca0015_sa.su2_exec = "SU2_CFD" + unst_inc_turb_naca0015_sa.timeout = 1600 + unst_inc_turb_naca0015_sa.tol = 0.00001 + unst_inc_turb_naca0015_sa.unsteady = True + test_list.append(unst_inc_turb_naca0015_sa) + ###################################### ### NICFD ### ###################################### diff --git a/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg b/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg new file mode 100644 index 000000000000..7905138ebaef --- /dev/null +++ b/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg @@ -0,0 +1,128 @@ +% ------------------------- PHYSICAL PROBLEM ----------------------------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= SA +MATH_PROBLEM= DIRECT +RESTART_SOL= NO + +% ------------------------- UNSTEADY SIMULATION -------------------------------% +% +UNSTEADY_SIMULATION= DUAL_TIME_STEPPING-2ND_ORDER +UNST_TIMESTEP= 0.016849% 25 time steps per period +UNST_TIME= 2.528% 6 periods +% +% Old driver +EXT_ITER= 151 +UNST_INT_ITER= 201 +% +% New driver +%SINGLEZONE_DRIVER= YES +%TIME_DOMAIN= YES +%TIME_ITER= 151 +%ITER= 201 +% ----------------------- DYNAMIC MESH DEFINITION -----------------------------% +% +GRID_MOVEMENT= RIGID_MOTION +MOTION_ORIGIN= (0.25 0.0 0.0) +PITCHING_OMEGA= (0.0 0.0 14.91675) +PITCHING_AMPL= (0.0 0.0 8.0) + +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +INC_DENSITY_MODEL= CONSTANT +INC_ENERGY_EQUATION = NO +INC_DENSITY_INIT= 0.664527479 +INC_VELOCITY_INIT= ( 41.26140059, 5.79891168, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +INC_NONDIM= DIMENSIONAL +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.84592e-05 + +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( airfoil, 0.0 ) +MARKER_FAR= ( farfield ) +MARKER_PLOTTING= ( airfoil ) +MARKER_MONITORING= ( airfoil ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +CFL_ADAPT= NO + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-14 +LINEAR_SOLVER_ITER= 2 +%LINEAR_SOLVER_ERROR= 1E-8 +%LINEAR_SOLVER_ITER= 2000 + +% ----------------------- SLOPE LIMITER DEFINITION ----------------------------% +% +VENKAT_LIMITER_COEFF= 0.03 +LIMITER_ITER= 99999 +%ADJ_SHARP_LIMITER_COEFF= 3.0 +%REF_SHARP_EDGES= 3.0 +%SENS_REMOVE_SHARP= NO + +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 0 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= JST +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% Runge-Kutta alpha coefficients +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +MUSCL_TURB= YES +SLOPE_LIMITER_TURB= NONE +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_CRITERIA= RESIDUAL +RESIDUAL_REDUCTION= 3 +RESIDUAL_MINVAL= -12 +STARTCONV_ITER= 1 +CAUCHY_ELEMS= 100 +CAUCHY_FUNC_FLOW= DRAG + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= mesh_naca0015_ogrid_m151_cvc_v2_147.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= mesh_out.su2 +SOLUTION_FLOW_FILENAME= solution_flow.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +OUTPUT_FORMAT= PARAVIEW_BINARY +CONV_FILENAME= history +RESTART_FLOW_FILENAME= solution_flow.dat +RESTART_ADJ_FILENAME= solution_adj.dat +VOLUME_FLOW_FILENAME= flow +VOLUME_ADJ_FILENAME= adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FLOW_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +WRT_SOL_FREQ= 100 +WRT_SOL_FREQ_DUALTIME= 1 +WRT_CON_FREQ= 1 +WRT_CON_FREQ_DUALTIME= 10 +