From d657e3f2148d81623ad8af7e049a8bd893c38599 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 26 Sep 2018 22:35:08 +0100 Subject: [PATCH 01/31] Initial implementation of incompressible ALE --- .travis.yml | 4 +- Common/src/config_structure.cpp | 4 +- SU2_CFD/include/numerics_structure.hpp | 34 +++++ SU2_CFD/src/driver_structure.cpp | 3 +- SU2_CFD/src/numerics_direct_mean_inc.cpp | 187 ++++++++++++++++++++++- SU2_CFD/src/solver_direct_mean_inc.cpp | 14 +- 6 files changed, 235 insertions(+), 11 deletions(-) diff --git a/.travis.yml b/.travis.yml index 8dce14b32153..7facb047d34e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,11 +12,11 @@ compiler: notifications: email: recipients: - - su2code-dev@lists.stanford.edu + - charanya.crome09@ic.ac.uk branches: only: - - develop + - feature_incompressible_ale python: - 2.7 diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 06ed8dde12ad..9e3a1ae44c32 100755 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -3774,9 +3774,9 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ /*--- Grid motion is not yet supported with the incompressible solver. ---*/ - if ((Kind_Regime == INCOMPRESSIBLE) && (Grid_Movement)) { + /*if ((Kind_Regime == INCOMPRESSIBLE) && (Grid_Movement)) { SU2_MPI::Error("Support for grid movement not yet implemented for incompressible flows.", CURRENT_FUNCTION); - } + }*/ /*--- Assert that there are two markers being analyzed if the pressure drop objective function is selected. ---*/ diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index 75c086314dbe..2130f3d18ad2 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5049,6 +5049,40 @@ class CSourceIncBodyForce : public CNumerics { }; +/*! + + * \class CSourceIncRotatingFrame_Flow + * \brief Class for a rotating frame source term. + * \ingroup SourceDiscr + * \author F. Palacios, T. Economon, C. Venkatesan-Crome + */ +class CSourceIncRotatingFrame_Flow : public CNumerics { +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/src/driver_structure.cpp b/SU2_CFD/src/driver_structure.cpp index 401784bf97df..fcbea340c342 100644 --- a/SU2_CFD/src/driver_structure.cpp +++ b/SU2_CFD/src/driver_structure.cpp @@ -2005,7 +2005,8 @@ void CDriver::Numerics_Preprocessing(CNumerics *****numerics_container, else if (incompressible && (config->GetKind_DensityModel() == BOUSSINESQ)) numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceBoussinesq(nDim, nVar_Flow, config); else if (config->GetRotating_Frame() == YES) - numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config); + if (incompressible) numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncRotatingFrame_Flow(nDim, nVar_Flow, config); + else numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceRotatingFrame_Flow(nDim, nVar_Flow, config); else if (config->GetAxisymmetric() == YES) if (incompressible) numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config); else numerics_container[val_iInst][iMGlevel][FLOW_SOL][SOURCE_FIRST_TERM] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config); diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index cd84c47c0d13..7fc61a5813d7 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -84,7 +84,10 @@ 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); @@ -151,6 +154,16 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J GetInviscidIncProjFlux(&DensityInc_j, Velocity_j, &Pressure_j, &BetaInc2_j, &Enthalpy_j, Normal, ProjFlux_j); + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (grid_movement) { + ProjGridVel = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*UnitNormal[iDim]; + } + ProjVelocity -= ProjGridVel; + } + /*--- Eigenvalues of the preconditioned system ---*/ if (nDim == 2) { @@ -226,6 +239,32 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } } + /*--- Jacobian contributions due to grid motion ---*/ + + if (grid_movement) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = Density_i; U_j[0] = Density_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; + } + U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; + + ProjVelocity = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Implicit terms ---*/ + /*if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + }*/ + } + } + AD::SetPreaccOut(val_residual, nVar); AD::EndPreacc(); } @@ -275,6 +314,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, ProjVelocity = 0.0; + /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -335,6 +377,30 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } + /*--- Adjustment due to grid motion ---*/ + + if (grid_movement) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = Density_i; U_j[0] = Density_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; + } + U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; + + ProjVelocity = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar] + U_j[iVar]); + /*if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + }*/ + } + } + /*--- Computes differences between Laplacians and conservative variables ---*/ for (iVar = 0; iVar < nVar; iVar++) { @@ -352,6 +418,16 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ SoundSpeed_i = sqrt(BetaInc2_i*Area*Area); SoundSpeed_j = sqrt(BetaInc2_j*Area*Area); + /*--- Adjustment due to mesh motion ---*/ + + if (grid_movement) { + 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; @@ -439,6 +515,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]; @@ -501,6 +580,30 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } + /*--- Adjustment due to grid motion ---*/ + + if (grid_movement) { + + /*--- Recompute conservative variables ---*/ + + U_i[0] = Density_i; U_j[0] = Density_j; + for (iDim = 0; iDim < nDim; iDim++) { + U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; + } + U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; + + ProjVelocity = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + /*if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + }*/ + } + } + /*--- Computes differences btw. conservative variables ---*/ for (iVar = 0; iVar < nVar; iVar++) @@ -516,6 +619,15 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ SoundSpeed_i = sqrt(BetaInc2_i*Area*Area); SoundSpeed_j = sqrt(BetaInc2_j*Area*Area); + /*--- Adjustment due to grid motion ---*/ + if (grid_movement) { + 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; @@ -889,6 +1001,79 @@ 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) { + + Gamma = config->GetGamma(); + Gamma_Minus_One = Gamma - 1.0; + +} + +CSourceIncRotatingFrame_Flow::~CSourceIncRotatingFrame_Flow(void) { } + +void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { + + + unsigned short iDim, iVar, jVar; + su2double Omega[3] = {0,0,0}, Momentum[3] = {0,0,0}, Velocity_i[3] = {0,0,0}; + + bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + + /*--- Retrieve the angular velocity vector from config. ---*/ + + Omega[0] = config->GetRotation_Rate_X(config->GetiZone())/config->GetOmega_Ref(); + Omega[1] = config->GetRotation_Rate_Y(config->GetiZone())/config->GetOmega_Ref(); + Omega[2] = config->GetRotation_Rate_Z(config->GetiZone())/config->GetOmega_Ref(); + + /*--- Primitive variables at point i and j ---*/ + + DensityInc_i = V_i[nDim+2]; + + for (iDim = 0; iDim < nDim; iDim++) { + Velocity_i[iDim] = V_i[iDim+1]; + } + + /*--- Get the momentum vector at the current node. ---*/ + + for (iDim = 0; iDim < nDim; iDim++) { + Momentum[iDim] = DensityInc_i*Velocity_i[iDim]; + } + + /*--- Calculate rotating frame source term 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 1d576aefa240..6abdd9cf44e2 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -2921,11 +2921,15 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Load the conservative variables ---*/ + /*--- Load the primitive variables ---*/ - numerics->SetConservative(node[iPoint]->GetSolution(), - node[iPoint]->GetSolution()); + numerics->SetPrimitive(node[iPoint]->GetPrimitive(), NULL); + /*--- Set incompressible density ---*/ + + numerics->SetDensity(node[iPoint]->GetDensity(), + node[iPoint]->GetDensity()); + /*--- Load the volume of the dual mesh cell ---*/ numerics->SetVolume(geometry->node[iPoint]->GetVolume()); @@ -5959,7 +5963,7 @@ 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; LinSysRes.AddBlock(iPoint, Residual); @@ -5973,7 +5977,7 @@ 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; LinSysRes.SubtractBlock(jPoint, Residual); From 6a15d717779fd264b8ea2f5368de81786abb5d55 Mon Sep 17 00:00:00 2001 From: cvencro Date: Sun, 30 Sep 2018 08:47:04 +0100 Subject: [PATCH 02/31] Add FSI bug fixes --- SU2_CFD/include/driver_structure.hpp | 3 +- SU2_CFD/src/driver_structure.cpp | 73 +++++++++++++++++------ SU2_CFD/src/solver_adjoint_discrete.cpp | 3 + SU2_CFD/src/solver_adjoint_elasticity.cpp | 7 +++ 4 files changed, 68 insertions(+), 18 deletions(-) diff --git a/SU2_CFD/include/driver_structure.hpp b/SU2_CFD/include/driver_structure.hpp index 3362c5e34fd1..004b6fb1383d 100644 --- a/SU2_CFD/include/driver_structure.hpp +++ b/SU2_CFD/include/driver_structure.hpp @@ -1188,7 +1188,8 @@ class CDiscAdjFSIDriver : public CFSIDriver { structure_criteria, structure_criteria_rel; - + bool filterCrossTerm; + enum OF_KIND{ NO_OBJECTIVE_FUNCTION = 0, /*!< \brief Indicates that there is no objective function. */ FLOW_OBJECTIVE_FUNCTION = 1, /*!< \brief Indicates that the objective function is only flow-dependent. */ diff --git a/SU2_CFD/src/driver_structure.cpp b/SU2_CFD/src/driver_structure.cpp index fcbea340c342..08b388fac1e4 100644 --- a/SU2_CFD/src/driver_structure.cpp +++ b/SU2_CFD/src/driver_structure.cpp @@ -5795,6 +5795,8 @@ CDiscAdjFSIDriver::CDiscAdjFSIDriver(char* confFile, RecordingState = 0; CurrentRecording = 0; + filterCrossTerm = false; + switch (config_container[ZONE_0]->GetKind_ObjFunc()){ case DRAG_COEFFICIENT: case LIFT_COEFFICIENT: @@ -6377,6 +6379,8 @@ void CDiscAdjFSIDriver::Iterate_Direct(unsigned short ZONE_FLOW, unsigned short void CDiscAdjFSIDriver::Fluid_Iteration_Direct(unsigned short ZONE_FLOW, unsigned short ZONE_STRUCT) { + bool turbulent = (config_container[ZONE_FLOW]->GetKind_Solver() == DISC_ADJ_RANS); + /*-----------------------------------------------------------------*/ /*------------------- Set Dependency on Geometry ------------------*/ /*-----------------------------------------------------------------*/ @@ -6387,6 +6391,11 @@ void CDiscAdjFSIDriver::Fluid_Iteration_Direct(unsigned short ZONE_FLOW, unsigne solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->Preprocessing(geometry_container[ZONE_FLOW][INST_0][MESH_0],solver_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); + if(turbulent){ + solver_container[ZONE_FLOW][INST_0][MESH_0][TURB_SOL]->Postprocessing(geometry_container[ZONE_FLOW][INST_0][MESH_0], solver_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW], MESH_0); + solver_container[ZONE_FLOW][INST_0][MESH_0][TURB_SOL]->Set_MPI_Solution(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); + } + /*-----------------------------------------------------------------*/ /*----------------- Iterate the flow solver -----------------------*/ /*---- Sets all the cross dependencies for the flow variables -----*/ @@ -6408,6 +6417,7 @@ void CDiscAdjFSIDriver::Fluid_Iteration_Direct(unsigned short ZONE_FLOW, unsigne void CDiscAdjFSIDriver::Structural_Iteration_Direct(unsigned short ZONE_FLOW, unsigned short ZONE_STRUCT) { + bool turbulent = (config_container[ZONE_FLOW]->GetKind_Solver() == DISC_ADJ_RANS); /*-----------------------------------------------------------------*/ /*---------- Set Dependencies on Geometry and Flow ----------------*/ @@ -6419,8 +6429,15 @@ void CDiscAdjFSIDriver::Structural_Iteration_Direct(unsigned short ZONE_FLOW, un solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->Set_MPI_Solution(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); + if(!filterCrossTerm) solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->Preprocessing(geometry_container[ZONE_FLOW][INST_0][MESH_0],solver_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); + if(!filterCrossTerm) + if(turbulent){ + solver_container[ZONE_FLOW][INST_0][MESH_0][TURB_SOL]->Postprocessing(geometry_container[ZONE_FLOW][INST_0][MESH_0], solver_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW], MESH_0); + solver_container[ZONE_FLOW][INST_0][MESH_0][TURB_SOL]->Set_MPI_Solution(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); + } + /*-----------------------------------------------------------------*/ /*-------------------- Transfer Tractions -------------------------*/ /*-----------------------------------------------------------------*/ @@ -6541,20 +6558,20 @@ void CDiscAdjFSIDriver::SetRecording(unsigned short ZONE_FLOW, AD::Reset(); - if (CurrentRecording != kind_recording && (CurrentRecording != NONE) ){ - - /*--- Clear indices ---*/ - - PrepareRecording(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); - - /*--- Clear indices of coupling variables ---*/ - - SetDependencies(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); - - /*--- Run one iteration while tape is passive - this clears all indices ---*/ - Iterate_Direct(ZONE_FLOW, ZONE_STRUCT, kind_recording); - - } +// if (CurrentRecording != kind_recording && (CurrentRecording != NONE) ){ +// +// /*--- Clear indices ---*/ +// +// PrepareRecording(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); +// +// /*--- Clear indices of coupling variables ---*/ +// +// SetDependencies(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); +// +// /*--- Run one iteration while tape is passive - this clears all indices ---*/ +// Iterate_Direct(ZONE_FLOW, ZONE_STRUCT, kind_recording); +// +// } /*--- Prepare for recording ---*/ @@ -6914,7 +6931,7 @@ void CDiscAdjFSIDriver::ExtractAdjoint(unsigned short ZONE_FLOW, if (kind_recording == FLOW_CROSS_TERM) { /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/ - + if(!filterCrossTerm) { solver_container[ZONE_FLOW][INST_0][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_CrossTerm(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); @@ -6922,7 +6939,23 @@ void CDiscAdjFSIDriver::ExtractAdjoint(unsigned short ZONE_FLOW, solver_container[ZONE_FLOW][INST_0][MESH_0][ADJTURB_SOL]->ExtractAdjoint_CrossTerm(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); } - + } else + { + // this is a temporary fix to remove some spurious values in the cross term before a final solution is found + // it is implemented here instead of inside ExtractAdjoint_CrossTerm to avoid touching the solver structure + unsigned long nPoint = geometry_container[ZONE_FLOW][INST_0][MESH_0]->GetnPoint(); + unsigned short nVar = solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->GetnVar(); + su2double* solution = new su2double[nVar]; + + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) + { + solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->node[iPoint]->GetAdjointSolution(solution); + for (unsigned short iVar = 0; iVar < nVar; ++iVar) + solver_container[ZONE_FLOW][INST_0][MESH_0][ADJFLOW_SOL]->node[iPoint]->SetCross_Term_Derivative(iVar, + solver_container[ZONE_FLOW][INST_0][MESH_0][ADJFLOW_SOL]->node[iPoint]->GetCross_Term_Derivative(iVar)-solution[iVar]); + } + delete[] solution; + } } if (kind_recording == FEM_CROSS_TERM_GEOMETRY) { @@ -7116,6 +7149,9 @@ void CDiscAdjFSIDriver::Iterate_Block_FlowOF(unsigned short ZONE_FLOW, /*--- Compute cross term (dS / dFv) ---*/ Iterate_Block(ZONE_FLOW, ZONE_STRUCT, FLOW_CROSS_TERM); + filterCrossTerm=true; + Iterate_Block(ZONE_FLOW, ZONE_STRUCT, FLOW_CROSS_TERM); + filterCrossTerm=false; /*--- Compute cross term (dS / dMv) ---*/ @@ -7162,7 +7198,10 @@ void CDiscAdjFSIDriver::Iterate_Block_StructuralOF(unsigned short ZONE_FLOW, /*--- Compute cross term (dS / dFv) ---*/ Iterate_Block(ZONE_FLOW, ZONE_STRUCT, FLOW_CROSS_TERM); - + filterCrossTerm=true; + Iterate_Block(ZONE_FLOW, ZONE_STRUCT, FLOW_CROSS_TERM); + filterCrossTerm=false; + /*--- Compute cross term (dS / dMv) ---*/ Iterate_Block(ZONE_FLOW, ZONE_STRUCT, GEOMETRY_CROSS_TERM); diff --git a/SU2_CFD/src/solver_adjoint_discrete.cpp b/SU2_CFD/src/solver_adjoint_discrete.cpp index f47f6b84ce68..6e5f9ba123ae 100755 --- a/SU2_CFD/src/solver_adjoint_discrete.cpp +++ b/SU2_CFD/src/solver_adjoint_discrete.cpp @@ -780,6 +780,9 @@ void CDiscAdjSolver::SetAdjoint_OutputMesh(CGeometry *geometry, CConfig *config) // Solution_Geometry[iDim] += node[iPoint]->GetDual_Time_Derivative_Geometry(iDim); // } // } + for (iDim = 0; iDim < nDim; iDim++){ + node[iPoint]->SetSensitivity(iDim, Solution_Geometry[iDim]); + } geometry->node[iPoint]->SetAdjointCoord(Solution_Geometry); } diff --git a/SU2_CFD/src/solver_adjoint_elasticity.cpp b/SU2_CFD/src/solver_adjoint_elasticity.cpp index 0260e3dbb6cb..b159a19e2e3f 100644 --- a/SU2_CFD/src/solver_adjoint_elasticity.cpp +++ b/SU2_CFD/src/solver_adjoint_elasticity.cpp @@ -1005,12 +1005,19 @@ void CDiscAdjFEASolver::ExtractAdjoint_CrossTerm_Geometry(CGeometry *geometry, C unsigned short iVar; unsigned long iPoint; + su2double relax = config->GetAitkenStatRelax(); + for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint solution ---*/ direct_solver->node[iPoint]->GetAdjointSolution(Solution); + /*--- Relax and set the solution ---*/ + + for(iVar = 0; iVar < nVar; iVar++) + Solution[iVar] = relax*Solution[iVar] + (1.0-relax)*node[iPoint]->GetGeometry_CrossTerm_Derivative(iVar); + for (iVar = 0; iVar < nVar; iVar++) node[iPoint]->SetGeometry_CrossTerm_Derivative(iVar, Solution[iVar]); } From 039ce3481e923fa146f607555579d6a82633d483 Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 21 Dec 2018 07:14:00 +0000 Subject: [PATCH 03/31] fix duplicated line --- SU2_CFD/src/solver_adjoint_elasticity.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/SU2_CFD/src/solver_adjoint_elasticity.cpp b/SU2_CFD/src/solver_adjoint_elasticity.cpp index 22a85fa9c5c7..b7d9fad0c510 100644 --- a/SU2_CFD/src/solver_adjoint_elasticity.cpp +++ b/SU2_CFD/src/solver_adjoint_elasticity.cpp @@ -1013,8 +1013,6 @@ void CDiscAdjFEASolver::ExtractAdjoint_CrossTerm_Geometry(CGeometry *geometry, C su2double relax = config->GetAitkenStatRelax(); - su2double relax = config->GetAitkenStatRelax(); - for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint solution ---*/ From c90727289c148a36bb38fab819cabe0cc41b4e8d Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 15 Feb 2019 20:11:23 +0000 Subject: [PATCH 04/31] Tidy temporary FSI fixes --- SU2_CFD/include/driver_structure.hpp | 3 +-- SU2_CFD/src/driver_structure.cpp | 31 ++++++++++------------- SU2_CFD/src/solver_adjoint_elasticity.cpp | 5 ---- 3 files changed, 15 insertions(+), 24 deletions(-) diff --git a/SU2_CFD/include/driver_structure.hpp b/SU2_CFD/include/driver_structure.hpp index bb57a5bf1a06..6b8691e24bb2 100644 --- a/SU2_CFD/include/driver_structure.hpp +++ b/SU2_CFD/include/driver_structure.hpp @@ -1218,8 +1218,7 @@ class CDiscAdjFSIDriver : public CDriver { structure_criteria, structure_criteria_rel; - bool filterCrossTerm; - + enum OF_KIND{ NO_OBJECTIVE_FUNCTION = 0, /*!< \brief Indicates that there is no objective function. */ FLOW_OBJECTIVE_FUNCTION = 1, /*!< \brief Indicates that the objective function is only flow-dependent. */ diff --git a/SU2_CFD/src/driver_structure.cpp b/SU2_CFD/src/driver_structure.cpp index 1d24994b3515..535dbc2b7822 100644 --- a/SU2_CFD/src/driver_structure.cpp +++ b/SU2_CFD/src/driver_structure.cpp @@ -5968,8 +5968,6 @@ CDiscAdjFSIDriver::CDiscAdjFSIDriver(char* confFile, RecordingState = 0; CurrentRecording = 0; - filterCrossTerm = false; - switch (config_container[ZONE_0]->GetKind_ObjFunc()){ case DRAG_COEFFICIENT: case LIFT_COEFFICIENT: @@ -6649,7 +6647,6 @@ void CDiscAdjFSIDriver::Structural_Iteration_Direct(unsigned short ZONE_FLOW, un solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->Set_MPI_Solution(geometry_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW]); - if(!filterCrossTerm) solver_container[ZONE_FLOW][INST_0][MESH_0][FLOW_SOL]->Preprocessing(geometry_container[ZONE_FLOW][INST_0][MESH_0],solver_container[ZONE_FLOW][INST_0][MESH_0], config_container[ZONE_FLOW], MESH_0, NO_RK_ITER, RUNTIME_FLOW_SYS, true); if (turbulent && !frozen_visc) { @@ -6777,20 +6774,20 @@ void CDiscAdjFSIDriver::SetRecording(unsigned short ZONE_FLOW, AD::Reset(); -// if (CurrentRecording != kind_recording && (CurrentRecording != NONE) ){ -// -// /*--- Clear indices ---*/ -// -// PrepareRecording(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); -// -// /*--- Clear indices of coupling variables ---*/ -// -// SetDependencies(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); -// -// /*--- Run one iteration while tape is passive - this clears all indices ---*/ -// Iterate_Direct(ZONE_FLOW, ZONE_STRUCT, kind_recording); -// -// } + if (CurrentRecording != kind_recording && (CurrentRecording != NONE) ){ + + /*--- Clear indices ---*/ + + PrepareRecording(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); + + /*--- Clear indices of coupling variables ---*/ + + SetDependencies(ZONE_FLOW, ZONE_STRUCT, ALL_VARIABLES); + + /*--- Run one iteration while tape is passive - this clears all indices ---*/ + Iterate_Direct(ZONE_FLOW, ZONE_STRUCT, kind_recording); + + } /*--- Prepare for recording ---*/ diff --git a/SU2_CFD/src/solver_adjoint_elasticity.cpp b/SU2_CFD/src/solver_adjoint_elasticity.cpp index 012fc11fc9aa..15f5acb240a5 100644 --- a/SU2_CFD/src/solver_adjoint_elasticity.cpp +++ b/SU2_CFD/src/solver_adjoint_elasticity.cpp @@ -1020,11 +1020,6 @@ void CDiscAdjFEASolver::ExtractAdjoint_CrossTerm_Geometry(CGeometry *geometry, C /*--- Relax and set the solution ---*/ - for(iVar = 0; iVar < nVar; iVar++) - Solution[iVar] = relax*Solution[iVar] + (1.0-relax)*node[iPoint]->GetGeometry_CrossTerm_Derivative(iVar); - - /*--- Relax and set the solution ---*/ - for(iVar = 0; iVar < nVar; iVar++) Solution[iVar] = relax*Solution[iVar] + (1.0-relax)*node[iPoint]->GetGeometry_CrossTerm_Derivative(iVar); From 6f4695580fc48d7c1524ce057da424db2675669e Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 27 Feb 2019 22:50:14 +0000 Subject: [PATCH 05/31] Corrected density term names --- SU2_CFD/include/numerics_structure.hpp | 4 +- SU2_CFD/src/numerics_direct_mean_inc.cpp | 217 ++++++++++++----------- SU2_CFD/src/solver_direct_mean_inc.cpp | 6 +- 3 files changed, 114 insertions(+), 113 deletions(-) diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index 320d1853e876..fe3588647692 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5231,8 +5231,7 @@ class CSourceIncBodyForce : public CNumerics { }; -/*! - +/*! * \class CSourceIncRotatingFrame_Flow * \brief Class for a rotating frame source term. * \ingroup SourceDiscr @@ -5264,7 +5263,6 @@ class CSourceIncRotatingFrame_Flow : public CNumerics { }; - /*! * \class CSourceBoussinesq * \brief Class for the source term integration of the Boussinesq approximation for incompressible flow. diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 805a9723813d..f20edb8a3a55 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -84,10 +84,10 @@ 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; - + + 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); @@ -122,6 +122,16 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J ProjVelocity += MeanVelocity[iDim]*Normal[iDim]; } + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (grid_movement) { + 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); @@ -154,16 +164,6 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J GetInviscidIncProjFlux(&DensityInc_j, Velocity_j, &Pressure_j, &BetaInc2_j, &Enthalpy_j, Normal, ProjFlux_j); - /*--- Projected velocity adjustment due to mesh motion ---*/ - - if (grid_movement) { - ProjGridVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*UnitNormal[iDim]; - } - ProjVelocity -= ProjGridVel; - } - /*--- Eigenvalues of the preconditioned system ---*/ if (nDim == 2) { @@ -226,6 +226,31 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } } + /*--- Jacobian contributions due to grid motion ---*/ + if (grid_movement) { + + /*--- 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]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Implicit terms ---*/ + if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + } + } + } + if (!energy) { val_residual[nDim+1] = 0.0; if (implicit) { @@ -239,32 +264,6 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } } - /*--- Jacobian contributions due to grid motion ---*/ - - if (grid_movement) { - - /*--- Recompute conservative variables ---*/ - - U_i[0] = Density_i; U_j[0] = Density_j; - for (iDim = 0; iDim < nDim; iDim++) { - U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; - } - U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; - - ProjVelocity = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - for (iVar = 0; iVar < nVar; iVar++) { - val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - - /*--- Implicit terms ---*/ - /*if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; - }*/ - } - } - AD::SetPreaccOut(val_residual, nVar); AD::EndPreacc(); } @@ -314,9 +313,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, ProjVelocity = 0.0; - + 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]; @@ -376,30 +375,31 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - - /*--- Adjustment due to grid motion ---*/ - if (grid_movement) { + /*--- Jacobian contributions due to grid motion ---*/ + if (grid_movement) { - /*--- Recompute conservative variables ---*/ + /*--- Recompute conservative variables ---*/ - U_i[0] = Density_i; U_j[0] = Density_j; - for (iDim = 0; iDim < nDim; iDim++) { - U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; - } - U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; - - ProjVelocity = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - for (iVar = 0; iVar < nVar; iVar++) { - val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar] + U_j[iVar]); - /*if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; - }*/ - } - } + 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]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Implicit terms ---*/ + if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + } + } + } /*--- Computes differences between Laplacians and conservative variables ---*/ @@ -415,19 +415,20 @@ 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 (grid_movement) { + 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); - /*--- Adjustment due to mesh motion ---*/ - - if (grid_movement) { - 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; @@ -515,9 +516,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; - + 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]; @@ -579,30 +580,30 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - - /*--- Adjustment due to grid motion ---*/ - if (grid_movement) { + /*--- Jacobian contributions due to grid motion ---*/ + if (grid_movement) { - /*--- Recompute conservative variables ---*/ + /*--- Recompute conservative variables ---*/ - U_i[0] = Density_i; U_j[0] = Density_j; - for (iDim = 0; iDim < nDim; iDim++) { - U_i[iDim+1] = Density_i*Velocity_i[iDim]; U_j[iDim+1] = Density_j*Velocity_j[iDim]; - } - U_i[nDim+1] = Density_i*Enthalpy_i; U_j[nDim+1] = Density_j*Enthalpy_j; - - ProjVelocity = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVelocity += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - for (iVar = 0; iVar < nVar; iVar++) { - val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); - /*if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; - }*/ - } - } + 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]; + for (iVar = 0; iVar < nVar; iVar++) { + val_residual[iVar] -= ProjVelocity * 0.5*(U_i[iVar]+U_j[iVar]); + + /*--- Implicit terms ---*/ + if (implicit) { + val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; + val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + } + } + } /*--- Computes differences btw. conservative variables ---*/ @@ -619,14 +620,16 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ SoundSpeed_i = sqrt(BetaInc2_i*Area*Area); SoundSpeed_j = sqrt(BetaInc2_j*Area*Area); - /*--- Adjustment due to grid motion ---*/ - if (grid_movement) { - 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; - } + /*--- Projected velocity adjustment due to mesh motion ---*/ + + if (grid_movement) { + 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; @@ -996,7 +999,7 @@ void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2d for (iDim = 0; iDim < nDim; iDim++) { Momentum[iDim] = DensityInc_i*Velocity_i[iDim]; - } + } /*--- Calculate rotating frame source term as ( Omega X Rho-U ) ---*/ diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 4398774e0b1a..c1ff5f6be65c 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -3043,12 +3043,12 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Load the primitive variables ---*/ numerics->SetPrimitive(node[iPoint]->GetPrimitive(), NULL); - + /*--- Set incompressible density ---*/ numerics->SetDensity(node[iPoint]->GetDensity(), node[iPoint]->GetDensity()); - + /*--- Load the volume of the dual mesh cell ---*/ numerics->SetVolume(geometry->node[iPoint]->GetVolume()); @@ -6509,7 +6509,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - for (iVar = 1; iVar < nVar; iVar++) { + for (iVar = 0; iVar < nVar; iVar++) { if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) Jacobian_i[iVar][iVar] = Volume_nP1/TimeStep; if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) From 74666b9e75c3bb8419163a7ed8e9c3f07ec696b8 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 13 Mar 2019 21:24:20 +0000 Subject: [PATCH 06/31] Temporarily ignore GCL --- SU2_CFD/src/solver_direct_mean_inc.cpp | 293 +++++++++++++------------ 1 file changed, 148 insertions(+), 145 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index c1ff5f6be65c..8770ab7efbf4 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6238,7 +6238,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the dual time-stepping source term for static meshes ---*/ - if (!grid_movement) { + // if (!grid_movement) { /*--- Loop over all nodes (excluding halos) ---*/ @@ -6321,208 +6321,211 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver } } - } + // } - else { + // else { - /*--- For unsteady flows on dynamic meshes (rigidly transforming or - dynamically deforming), the Geometric Conservation Law (GCL) should be - satisfied in conjunction with the ALE formulation of the governing - equations. The GCL prevents accuracy issues caused by grid motion, i.e. - a uniform free-stream should be preserved through a moving grid. First, - we will loop over the edges and boundaries to compute the GCL component - of the dual time source term that depends on grid velocities. ---*/ + // /*--- For unsteady flows on dynamic meshes (rigidly transforming or + // dynamically deforming), the Geometric Conservation Law (GCL) should be + // satisfied in conjunction with the ALE formulation of the governing + // equations. The GCL prevents accuracy issues caused by grid motion, i.e. + // a uniform free-stream should be preserved through a moving grid. First, + // we will loop over the edges and boundaries to compute the GCL component + // of the dual time source term that depends on grid velocities. ---*/ - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Initialize the Residual / Jacobian container to zero. ---*/ + // for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; + // /*--- Initialize the Residual / Jacobian container to zero. ---*/ - /*--- Get indices for nodes i & j plus the face normal ---*/ + // for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - iPoint = geometry->edge[iEdge]->GetNode(0); - jPoint = geometry->edge[iEdge]->GetNode(1); - Normal = geometry->edge[iEdge]->GetNormal(); + // /*--- Get indices for nodes i & j plus the face normal ---*/ - /*--- Grid velocities stored at nodes i & j ---*/ + // iPoint = geometry->edge[iEdge]->GetNode(0); + // jPoint = geometry->edge[iEdge]->GetNode(1); + // Normal = geometry->edge[iEdge]->GetNormal(); - GridVel_i = geometry->node[iPoint]->GetGridVel(); - GridVel_j = geometry->node[jPoint]->GetGridVel(); + // /*--- Grid velocities stored at nodes i & j ---*/ - /*--- Compute the GCL term by averaging the grid velocities at the - edge mid-point and dotting with the face normal. ---*/ + // GridVel_i = geometry->node[iPoint]->GetGridVel(); + // GridVel_j = geometry->node[jPoint]->GetGridVel(); - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + // /*--- Compute the GCL term by averaging the grid velocities at the + // edge mid-point and dotting with the face normal. ---*/ - /*--- Compute the GCL component of the source term for node i ---*/ + // Residual_GCL = 0.0; + // for (iDim = 0; iDim < nDim; iDim++) + // Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + // /*if (iEdge == geometry->GetnEdge() - 1) { + // cout << "CVC: Debug: CIncEulerSolver::SetResidual_DualTime" << endl; + // cout << "CVC: Debug: Residual_GCL = " << Residual_GCL << endl; + // }*/ + // /*--- Compute the GCL component of the source term for node i ---*/ - V_time_n = node[iPoint]->GetSolution_time_n(); + // V_time_n = node[iPoint]->GetSolution_time_n(); - /*--- Access the density and Cp at this node (constant for now). ---*/ + // /*--- Access the density and Cp at this node (constant for now). ---*/ - Density = node[iPoint]->GetDensity(); - Cp = node[iPoint]->GetSpecificHeatCp(); + // Density = node[iPoint]->GetDensity(); + // Cp = node[iPoint]->GetSpecificHeatCp(); - /*--- Compute the conservative variable vector for all time levels. ---*/ + // /*--- Compute the conservative variable vector for all time levels. ---*/ - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + // U_time_n[0] = Density; + // for (iDim = 0; iDim < nDim; iDim++) { + // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + // } + // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; - LinSysRes.AddBlock(iPoint, Residual); + // for (iVar = 0; iVar < nVar; iVar++) + // Residual[iVar] = U_time_n[iVar]*Residual_GCL; + // LinSysRes.AddBlock(iPoint, Residual); - /*--- Compute the GCL component of the source term for node j ---*/ + // /*--- Compute the GCL component of the source term for node j ---*/ - V_time_n = node[jPoint]->GetSolution_time_n(); + // V_time_n = node[jPoint]->GetSolution_time_n(); - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + // U_time_n[0] = Density; + // for (iDim = 0; iDim < nDim; iDim++) { + // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + // } + // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; - LinSysRes.SubtractBlock(jPoint, Residual); + // for (iVar = 0; iVar < nVar; iVar++) + // Residual[iVar] = U_time_n[iVar]*Residual_GCL; + // LinSysRes.SubtractBlock(jPoint, Residual); - } + // } - /*--- Loop over the boundary edges ---*/ + // /*--- Loop over the boundary edges ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + // for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + // for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - /*--- Initialize the Residual / Jacobian container to zero. ---*/ + // /*--- Initialize the Residual / Jacobian container to zero. ---*/ - for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; + // for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - /*--- Get the index for node i plus the boundary face normal ---*/ + // /*--- Get the index for node i plus the boundary face normal ---*/ - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + // iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + // Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - /*--- Grid velocities stored at boundary node i ---*/ + // /*--- Grid velocities stored at boundary node i ---*/ - GridVel_i = geometry->node[iPoint]->GetGridVel(); + // GridVel_i = geometry->node[iPoint]->GetGridVel(); - /*--- Compute the GCL term by dotting the grid velocity with the face - normal. The normal is negated to match the boundary convention. ---*/ + // /*--- Compute the GCL term by dotting the grid velocity with the face + // normal. The normal is negated to match the boundary convention. ---*/ - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; + // Residual_GCL = 0.0; + // for (iDim = 0; iDim < nDim; iDim++) + // Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; - /*--- Compute the GCL component of the source term for node i ---*/ + // /*--- Compute the GCL component of the source term for node i ---*/ - V_time_n = node[iPoint]->GetSolution_time_n(); + // V_time_n = node[iPoint]->GetSolution_time_n(); - /*--- Access the density and Cp at this node (constant for now). ---*/ + // /*--- Access the density and Cp at this node (constant for now). ---*/ - Density = node[iPoint]->GetDensity(); - Cp = node[iPoint]->GetSpecificHeatCp(); + // Density = node[iPoint]->GetDensity(); + // Cp = node[iPoint]->GetSpecificHeatCp(); - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + // U_time_n[0] = Density; + // for (iDim = 0; iDim < nDim; iDim++) { + // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + // } + // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; - LinSysRes.AddBlock(iPoint, Residual); + // for (iVar = 0; iVar < nVar; iVar++) + // Residual[iVar] = U_time_n[iVar]*Residual_GCL; + // LinSysRes.AddBlock(iPoint, Residual); - } - } + // } + // } - /*--- Loop over all nodes (excluding halos) to compute the remainder - of the dual time-stepping source term. ---*/ + // /*--- Loop over all nodes (excluding halos) to compute the remainder + // of the dual time-stepping source term. ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + // for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Initialize the Residual / Jacobian container to zero. ---*/ + // /*--- Initialize the Residual / Jacobian container to zero. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - Residual[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } + // for (iVar = 0; iVar < nVar; iVar++) { + // Residual[iVar] = 0.0; + // if (implicit) { + // for (jVar = 0; jVar < nVar; jVar++) + // Jacobian_i[iVar][jVar] = 0.0; + // } + // } - /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that - we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, - previous solutions that are stored in memory. ---*/ + // /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that + // we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, + // previous solutions that are stored in memory. ---*/ - V_time_nM1 = node[iPoint]->GetSolution_time_n1(); - V_time_n = node[iPoint]->GetSolution_time_n(); - V_time_nP1 = node[iPoint]->GetSolution(); + // V_time_nM1 = node[iPoint]->GetSolution_time_n1(); + // V_time_n = node[iPoint]->GetSolution_time_n(); + // V_time_nP1 = node[iPoint]->GetSolution(); - /*--- Access the density and Cp at this node (constant for now). ---*/ + // /*--- Access the density and Cp at this node (constant for now). ---*/ - Density = node[iPoint]->GetDensity(); - Cp = node[iPoint]->GetSpecificHeatCp(); + // Density = node[iPoint]->GetDensity(); + // Cp = node[iPoint]->GetSpecificHeatCp(); - /*--- Compute the conservative variable vector for all time levels. ---*/ + // /*--- Compute the conservative variable vector for all time levels. ---*/ - U_time_nM1[0] = Density; - U_time_n[0] = Density; - U_time_nP1[0] = Density; + // U_time_nM1[0] = Density; + // U_time_n[0] = Density; + // U_time_nP1[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; - } + // for (iDim = 0; iDim < nDim; iDim++) { + // U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; + // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + // U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; + // } - U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; + // U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; + // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + // U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; - /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming - grids, the volumes will change. On rigidly transforming grids, the - volumes will remain constant. ---*/ + // /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming + // grids, the volumes will change. On rigidly transforming grids, the + // volumes will remain constant. ---*/ - Volume_nM1 = geometry->node[iPoint]->GetVolume_nM1(); - Volume_nP1 = geometry->node[iPoint]->GetVolume(); + // Volume_nM1 = geometry->node[iPoint]->GetVolume_nM1(); + // Volume_nP1 = geometry->node[iPoint]->GetVolume(); - /*--- Compute the dual time-stepping source residual. Due to the - introduction of the GCL term above, the remainder of the source residual - due to the time discretization has a new form.---*/ + // /*--- Compute the dual time-stepping source residual. Due to the + // introduction of the GCL term above, the remainder of the source residual + // due to the time discretization has a new form.---*/ - for (iVar = 0; iVar < nVar; iVar++) { - if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) - Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); - if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) - Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) - + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); - } + // for (iVar = 0; iVar < nVar; iVar++) { + // if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) + // Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); + // if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) + // Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) + // + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); + // } - /*--- Store the residual and compute the Jacobian contribution due - to the dual time source term. ---*/ + // /*--- Store the residual and compute the Jacobian contribution due + // to the dual time source term. ---*/ - LinSysRes.AddBlock(iPoint, Residual); - if (implicit) { - for (iVar = 0; 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); - } - 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]; + // LinSysRes.AddBlock(iPoint, Residual); + // if (implicit) { + // for (iVar = 0; 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); + // } + // 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); - } - } - } + // Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); + // } + // } + // } } From bcf715cd05420054a0b6cc0ba53924c8c698f95c Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 20 Mar 2019 09:50:15 +0000 Subject: [PATCH 07/31] Reverting temporary comments --- SU2_CFD/src/solver_direct_mean_inc.cpp | 293 ++++++++++++------------- 1 file changed, 145 insertions(+), 148 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 8770ab7efbf4..c1ff5f6be65c 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6238,7 +6238,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the dual time-stepping source term for static meshes ---*/ - // if (!grid_movement) { + if (!grid_movement) { /*--- Loop over all nodes (excluding halos) ---*/ @@ -6321,211 +6321,208 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver } } - // } + } - // else { + else { - // /*--- For unsteady flows on dynamic meshes (rigidly transforming or - // dynamically deforming), the Geometric Conservation Law (GCL) should be - // satisfied in conjunction with the ALE formulation of the governing - // equations. The GCL prevents accuracy issues caused by grid motion, i.e. - // a uniform free-stream should be preserved through a moving grid. First, - // we will loop over the edges and boundaries to compute the GCL component - // of the dual time source term that depends on grid velocities. ---*/ + /*--- For unsteady flows on dynamic meshes (rigidly transforming or + dynamically deforming), the Geometric Conservation Law (GCL) should be + satisfied in conjunction with the ALE formulation of the governing + equations. The GCL prevents accuracy issues caused by grid motion, i.e. + a uniform free-stream should be preserved through a moving grid. First, + we will loop over the edges and boundaries to compute the GCL component + of the dual time source term that depends on grid velocities. ---*/ - // for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { + + /*--- Initialize the Residual / Jacobian container to zero. ---*/ - // /*--- Initialize the Residual / Jacobian container to zero. ---*/ + for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - // for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; + /*--- Get indices for nodes i & j plus the face normal ---*/ - // /*--- Get indices for nodes i & j plus the face normal ---*/ + iPoint = geometry->edge[iEdge]->GetNode(0); + jPoint = geometry->edge[iEdge]->GetNode(1); + Normal = geometry->edge[iEdge]->GetNormal(); - // iPoint = geometry->edge[iEdge]->GetNode(0); - // jPoint = geometry->edge[iEdge]->GetNode(1); - // Normal = geometry->edge[iEdge]->GetNormal(); + /*--- Grid velocities stored at nodes i & j ---*/ - // /*--- Grid velocities stored at nodes i & j ---*/ + GridVel_i = geometry->node[iPoint]->GetGridVel(); + GridVel_j = geometry->node[jPoint]->GetGridVel(); - // GridVel_i = geometry->node[iPoint]->GetGridVel(); - // GridVel_j = geometry->node[jPoint]->GetGridVel(); + /*--- Compute the GCL term by averaging the grid velocities at the + edge mid-point and dotting with the face normal. ---*/ - // /*--- Compute the GCL term by averaging the grid velocities at the - // edge mid-point and dotting with the face normal. ---*/ + Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - // Residual_GCL = 0.0; - // for (iDim = 0; iDim < nDim; iDim++) - // Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - // /*if (iEdge == geometry->GetnEdge() - 1) { - // cout << "CVC: Debug: CIncEulerSolver::SetResidual_DualTime" << endl; - // cout << "CVC: Debug: Residual_GCL = " << Residual_GCL << endl; - // }*/ - // /*--- Compute the GCL component of the source term for node i ---*/ + /*--- Compute the GCL component of the source term for node i ---*/ - // V_time_n = node[iPoint]->GetSolution_time_n(); + V_time_n = node[iPoint]->GetSolution_time_n(); - // /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Access the density and Cp at this node (constant for now). ---*/ - // Density = node[iPoint]->GetDensity(); - // Cp = node[iPoint]->GetSpecificHeatCp(); + Density = node[iPoint]->GetDensity(); + Cp = node[iPoint]->GetSpecificHeatCp(); - // /*--- Compute the conservative variable vector for all time levels. ---*/ + /*--- Compute the conservative variable vector for all time levels. ---*/ - // U_time_n[0] = Density; - // for (iDim = 0; iDim < nDim; iDim++) { - // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - // } - // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + U_time_n[0] = Density; + for (iDim = 0; iDim < nDim; iDim++) { + U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + } + U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - // for (iVar = 0; iVar < nVar; iVar++) - // Residual[iVar] = U_time_n[iVar]*Residual_GCL; - // LinSysRes.AddBlock(iPoint, Residual); + for (iVar = 0; iVar < nVar; iVar++) + Residual[iVar] = U_time_n[iVar]*Residual_GCL; + LinSysRes.AddBlock(iPoint, Residual); - // /*--- Compute the GCL component of the source term for node j ---*/ + /*--- Compute the GCL component of the source term for node j ---*/ - // V_time_n = node[jPoint]->GetSolution_time_n(); + V_time_n = node[jPoint]->GetSolution_time_n(); - // U_time_n[0] = Density; - // for (iDim = 0; iDim < nDim; iDim++) { - // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - // } - // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + U_time_n[0] = Density; + for (iDim = 0; iDim < nDim; iDim++) { + U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + } + U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - // for (iVar = 0; iVar < nVar; iVar++) - // Residual[iVar] = U_time_n[iVar]*Residual_GCL; - // LinSysRes.SubtractBlock(jPoint, Residual); + for (iVar = 0; iVar < nVar; iVar++) + Residual[iVar] = U_time_n[iVar]*Residual_GCL; + LinSysRes.SubtractBlock(jPoint, Residual); - // } + } - // /*--- Loop over the boundary edges ---*/ + /*--- Loop over the boundary edges ---*/ - // for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - // for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - // /*--- Initialize the Residual / Jacobian container to zero. ---*/ + /*--- Initialize the Residual / Jacobian container to zero. ---*/ - // for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; + for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - // /*--- Get the index for node i plus the boundary face normal ---*/ + /*--- Get the index for node i plus the boundary face normal ---*/ - // iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - // Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - // /*--- Grid velocities stored at boundary node i ---*/ + /*--- Grid velocities stored at boundary node i ---*/ - // GridVel_i = geometry->node[iPoint]->GetGridVel(); + GridVel_i = geometry->node[iPoint]->GetGridVel(); - // /*--- Compute the GCL term by dotting the grid velocity with the face - // normal. The normal is negated to match the boundary convention. ---*/ + /*--- Compute the GCL term by dotting the grid velocity with the face + normal. The normal is negated to match the boundary convention. ---*/ - // Residual_GCL = 0.0; - // for (iDim = 0; iDim < nDim; iDim++) - // Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; + Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; - // /*--- Compute the GCL component of the source term for node i ---*/ + /*--- Compute the GCL component of the source term for node i ---*/ - // V_time_n = node[iPoint]->GetSolution_time_n(); + V_time_n = node[iPoint]->GetSolution_time_n(); - // /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Access the density and Cp at this node (constant for now). ---*/ - // Density = node[iPoint]->GetDensity(); - // Cp = node[iPoint]->GetSpecificHeatCp(); + Density = node[iPoint]->GetDensity(); + Cp = node[iPoint]->GetSpecificHeatCp(); - // U_time_n[0] = Density; - // for (iDim = 0; iDim < nDim; iDim++) { - // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - // } - // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + U_time_n[0] = Density; + for (iDim = 0; iDim < nDim; iDim++) { + U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + } + U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - // for (iVar = 0; iVar < nVar; iVar++) - // Residual[iVar] = U_time_n[iVar]*Residual_GCL; - // LinSysRes.AddBlock(iPoint, Residual); + for (iVar = 0; iVar < nVar; iVar++) + Residual[iVar] = U_time_n[iVar]*Residual_GCL; + LinSysRes.AddBlock(iPoint, Residual); - // } - // } + } + } - // /*--- Loop over all nodes (excluding halos) to compute the remainder - // of the dual time-stepping source term. ---*/ + /*--- Loop over all nodes (excluding halos) to compute the remainder + of the dual time-stepping source term. ---*/ - // for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - // /*--- Initialize the Residual / Jacobian container to zero. ---*/ + /*--- Initialize the Residual / Jacobian container to zero. ---*/ - // for (iVar = 0; iVar < nVar; iVar++) { - // Residual[iVar] = 0.0; - // if (implicit) { - // for (jVar = 0; jVar < nVar; jVar++) - // Jacobian_i[iVar][jVar] = 0.0; - // } - // } + for (iVar = 0; iVar < nVar; iVar++) { + Residual[iVar] = 0.0; + if (implicit) { + for (jVar = 0; jVar < nVar; jVar++) + Jacobian_i[iVar][jVar] = 0.0; + } + } - // /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that - // we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, - // previous solutions that are stored in memory. ---*/ + /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that + we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, + previous solutions that are stored in memory. ---*/ - // V_time_nM1 = node[iPoint]->GetSolution_time_n1(); - // V_time_n = node[iPoint]->GetSolution_time_n(); - // V_time_nP1 = node[iPoint]->GetSolution(); + V_time_nM1 = node[iPoint]->GetSolution_time_n1(); + V_time_n = node[iPoint]->GetSolution_time_n(); + V_time_nP1 = node[iPoint]->GetSolution(); - // /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Access the density and Cp at this node (constant for now). ---*/ - // Density = node[iPoint]->GetDensity(); - // Cp = node[iPoint]->GetSpecificHeatCp(); + Density = node[iPoint]->GetDensity(); + Cp = node[iPoint]->GetSpecificHeatCp(); - // /*--- Compute the conservative variable vector for all time levels. ---*/ + /*--- Compute the conservative variable vector for all time levels. ---*/ - // U_time_nM1[0] = Density; - // U_time_n[0] = Density; - // U_time_nP1[0] = Density; + U_time_nM1[0] = Density; + U_time_n[0] = Density; + U_time_nP1[0] = Density; - // for (iDim = 0; iDim < nDim; iDim++) { - // U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; - // U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - // U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; - // } + for (iDim = 0; iDim < nDim; iDim++) { + U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; + U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; + } - // U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; - // U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - // U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; + U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; + U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; + U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; - // /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming - // grids, the volumes will change. On rigidly transforming grids, the - // volumes will remain constant. ---*/ + /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming + grids, the volumes will change. On rigidly transforming grids, the + volumes will remain constant. ---*/ - // Volume_nM1 = geometry->node[iPoint]->GetVolume_nM1(); - // Volume_nP1 = geometry->node[iPoint]->GetVolume(); + Volume_nM1 = geometry->node[iPoint]->GetVolume_nM1(); + Volume_nP1 = geometry->node[iPoint]->GetVolume(); - // /*--- Compute the dual time-stepping source residual. Due to the - // introduction of the GCL term above, the remainder of the source residual - // due to the time discretization has a new form.---*/ + /*--- Compute the dual time-stepping source residual. Due to the + introduction of the GCL term above, the remainder of the source residual + due to the time discretization has a new form.---*/ - // for (iVar = 0; iVar < nVar; iVar++) { - // if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) - // Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); - // if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) - // Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) - // + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); - // } + for (iVar = 0; iVar < nVar; iVar++) { + if (config->GetUnsteady_Simulation() == DT_STEPPING_1ST) + Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); + if (config->GetUnsteady_Simulation() == DT_STEPPING_2ND) + Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) + + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); + } - // /*--- Store the residual and compute the Jacobian contribution due - // to the dual time source term. ---*/ + /*--- Store the residual and compute the Jacobian contribution due + to the dual time source term. ---*/ - // LinSysRes.AddBlock(iPoint, Residual); - // if (implicit) { - // for (iVar = 0; 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); - // } - // 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]; + LinSysRes.AddBlock(iPoint, Residual); + if (implicit) { + for (iVar = 0; 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); + } + 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); - // } - // } - // } + Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); + } + } + } } From 8eb7d1d0c3a51e2508e0998c1f585a212dc461d8 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 20 Mar 2019 10:03:36 +0000 Subject: [PATCH 08/31] Preconditioning terms added to the Jacobian contribution due to the dual-time source term --- SU2_CFD/src/solver_direct_mean_inc.cpp | 164 +++++++++++++++++++++---- 1 file changed, 140 insertions(+), 24 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index c1ff5f6be65c..5ab334238bd7 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6220,10 +6220,12 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Local variables ---*/ - unsigned short iVar, jVar, iMarker, iDim; + unsigned short iVar, jVar, iMarker, iDim, jDim; unsigned long iPoint, jPoint, iEdge, iVertex; su2double Density, Cp; + su2double BetaInc2, dRhodT, Temperature, oneOverCp; + su2double Velocity[3] = {0.0,0.0,0.0}; su2double *V_time_nM1, *V_time_n, *V_time_nP1; su2double U_time_nM1[5], U_time_n[5], U_time_nP1[5]; su2double Volume_nM1, Volume_nP1, TimeStep; @@ -6231,6 +6233,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool grid_movement = config->GetGrid_Movement(); + bool variable_density = (config->GetKind_DensityModel() == VARIABLE); + bool energy = config->GetEnergy_Equation(); /*--- Store the physical time step ---*/ @@ -6263,11 +6267,27 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = node[iPoint]->GetSolution_time_n(); V_time_nP1 = node[iPoint]->GetSolution(); - /*--- Access the density and Cp at this node (constant for now). ---*/ - + /*--- Access the primitive variables at this node. ---*/ + Density = node[iPoint]->GetDensity(); + BetaInc2 = node[iPoint]->GetBetaInc2(); Cp = node[iPoint]->GetSpecificHeatCp(); - + oneOverCp = 1.0/Cp; + 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; + } + /*--- Compute the conservative variable vector for all time levels. ---*/ U_time_nM1[0] = Density; @@ -6301,21 +6321,60 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver Residual[iVar] = ( 3.0*U_time_nP1[iVar] - 4.0*U_time_n[iVar] +1.0*U_time_nM1[iVar])*Volume_nP1 / (2.0*TimeStep); } - + + if (!energy) Residual[nDim+1] = 0.0; + /*--- Store the residual and compute the Jacobian contribution due to the dual time source term. ---*/ 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] = (Volume_nP1*3.0)/(2.0*TimeStep); + + /*--- 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; } + + Jacobian_i[0][nDim+1] = dRhodT; 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_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) + 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; + } + } Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } @@ -6361,11 +6420,27 @@ 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). ---*/ - + /*--- Access the primitive variables at this node. ---*/ + Density = node[iPoint]->GetDensity(); + BetaInc2 = node[iPoint]->GetBetaInc2(); Cp = node[iPoint]->GetSpecificHeatCp(); - + oneOverCp = 1.0/Cp; + 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; + } + /*--- Compute the conservative variable vector for all time levels. ---*/ U_time_n[0] = Density; @@ -6376,6 +6451,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); /*--- Compute the GCL component of the source term for node j ---*/ @@ -6390,6 +6467,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.SubtractBlock(jPoint, Residual); } @@ -6436,6 +6515,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); } @@ -6506,19 +6587,54 @@ 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 = 0; 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); + + /*--- 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; } + + Jacobian_i[0][nDim+1] = dRhodT; 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_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) + 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; + } + } Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); } } From 914a90b899b2ccc5ba0f0580048389d8f50dda23 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Mon, 25 Mar 2019 16:01:41 +0100 Subject: [PATCH 09/31] Only commented changes. Disturbed Initialzation. Reflected state Euler_Wall. --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 1 + SU2_CFD/src/solver_direct_mean_inc.cpp | 144 ++++++++++++++++++++++- 2 files changed, 144 insertions(+), 1 deletion(-) mode change 100644 => 100755 SU2_CFD/src/numerics_direct_mean_inc.cpp diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp old mode 100644 new mode 100755 index f20edb8a3a55..e889e6030056 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -313,6 +313,7 @@ CCentJSTInc_Flow::~CCentJSTInc_Flow(void) { void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { + //Preaccumulation? 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; diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index f55763fcf955..e037cda0fb2a 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -551,6 +551,22 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned /*--- Initialize the solution to the far-field state everywhere. ---*/ + //TK:: Disturb the initial solution slightly to see whether ALE stuff is doing anything at all + //if(rank == MASTER_NODE) cout << "Disturbing initial solution slightly, CIncEulerSolver::CIncEulerSolver" << endl; + //su2double tmpPressure_Inf = Pressure_Inf + 0.2*Pressure_Inf; + //su2double *tmpVelocity_Inf = new su2double[nDim]; + //for (iDim = 0; iDim < nDim; iDim++) { + // if (Velocity_Inf[iDim] == 0.0) { + // tmpVelocity_Inf[iDim] = 12.3; + // } else { + // tmpVelocity_Inf[iDim] = Velocity_Inf[iDim] + 0.2*Velocity_Inf[iDim]; + // } + //} + //su2double tmpTemperature_Inf = Temperature_Inf + 0.2*Temperature_Inf; + //TK:: end + + //for (iPoint = 0; iPoint < nPoint; iPoint++) + // node[iPoint] = new CIncEulerVariable(tmpPressure_Inf, tmpVelocity_Inf, tmpTemperature_Inf, nDim, nVar, config); for (iPoint = 0; iPoint < nPoint; iPoint++) node[iPoint] = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nDim, nVar, config); @@ -5376,6 +5392,132 @@ void CIncEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { } +//void CIncEulerSolver::BC_Euler_Wall(CGeometry *geometry, +// CSolver **solver_container, +// CNumerics *conv_numerics, +// CConfig *config, +// unsigned short val_marker) { +// +// unsigned short iDim, iVar; +// unsigned long iVertex, iPoint; +// +// bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); +// bool grid_movement = config->GetGrid_Movement(); +// +// /*--- Allocation of variables necessary for convective fluxes. ---*/ +// su2double Area, ProjVelocity_i; +// su2double *V_reflected, *V_domain; +// su2double *Normal = new su2double[nDim]; +// su2double *UnitNormal = new su2double[nDim]; +// +// su2double *GridVel = NULL; +// su2double ProjGridVel = 0.0; +// su2double *Velocity_b; +// Velocity_b = new su2double[nDim]; +// su2double Density_b, BetaInc2_b, Enthalpy_b; +// +// /*--- Loop over all the vertices on this boundary marker. ---*/ +// for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { +// +// iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); +// +// /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ +// if (geometry->node[iPoint]->GetDomain()) { +// +// /*-------------------------------------------------------------------------------*/ +// /*--- Step 1: For the convective fluxes, create a reflected state of the ---*/ +// /*--- Primitive variables by copying all interior values to the ---*/ +// /*--- reflected. Only the velocity is mirrored along the symmetry ---*/ +// /*--- axis. Based on the Upwind_Residual routine. ---*/ +// /*-------------------------------------------------------------------------------*/ +// +// /*--- Normal vector for a random vertex (zero) on this marker (negate for outward convention). ---*/ +// geometry->vertex[val_marker][0]->GetNormal(Normal); +// for (iDim = 0; iDim < nDim; iDim++) +// Normal[iDim] = -Normal[iDim]; +// +// /*--- Compute unit normal, to be used for unit tangential, projected velocity and velocity component gradients. ---*/ +// Area = 0.0; +// for (iDim = 0; iDim < nDim; iDim++) +// Area += Normal[iDim]*Normal[iDim]; +// Area = sqrt (Area); +// +// for (iDim = 0; iDim < nDim; iDim++) +// UnitNormal[iDim] = -Normal[iDim]/Area; +// +// /*--- Allocate the reflected state at the symmetry boundary. ---*/ +// V_reflected = GetCharacPrimVar(val_marker, iVertex); +// +// /*--- Grid movement ---*/ +// //if (grid_movement) +// // conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); +// +// /*--- Normal vector for this vertex (negate for outward convention). ---*/ +// geometry->vertex[val_marker][iVertex]->GetNormal(Normal); +// for (iDim = 0; iDim < nDim; iDim++) +// Normal[iDim] = -Normal[iDim]; +// //conv_numerics->SetNormal(Normal); +// +// /*--- Get current solution at this boundary node ---*/ +// V_domain = node[iPoint]->GetPrimitive(); +// +// /*--- Set the reflected state based on the boundary node. Scalars are copied and +// the velocity is mirrored along the symmetry boundary, i.e. the velocity in +// normal direction is substracted twice. ---*/ +// for(iVar = 0; iVar < nPrimVar; iVar++) +// V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); +// +// /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from +// velocity in normal direction: v_r = v - 2 (v*n)n ---*/ +// ProjVelocity_i = 0.0; +// for (iDim = 0; iDim < nDim; iDim++) +// ProjVelocity_i += node[iPoint]->GetVelocity(iDim)*UnitNormal[iDim]; +// +// for (iDim = 0; iDim < nDim; iDim++) +// V_reflected[iDim+1] = node[iPoint]->GetVelocity(iDim) - ProjVelocity_i*UnitNormal[iDim]; +// +// +// if (grid_movement) { +// GridVel = geometry->node[iPoint]->GetGridVel(); +// ProjGridVel = 0.0; +// for (iDim = 0; iDim < nDim; iDim++) ProjGridVel += GridVel[iDim]*UnitNormal[iDim]; +// for (iDim = 0; iDim < nDim; iDim++) V_reflected[iDim+1] += GridVel[iDim] - ProjGridVel * UnitNormal[iDim]; +// } +// +// for (iDim = 0; iDim < nDim; iDim++) +// Velocity_b[iDim] = V_reflected[iDim+1]; +// +// /*--- Set Primitive and Secondary for numerics class. ---*/ +// //conv_numerics->SetPrimitive(V_domain, V_reflected); +// //conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); +// +// /*--- Compute the residual using an upwind scheme. ---*/ +// //conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); +// +// //conv_numerics->GetInviscidProjFlux(&Density_b, Velocity_b, &Pressure_b, &Enthalpy_b, NormalArea, Residual); +// +// Density_b = node[iPoint]->GetDensity(); +// BetaInc2_b = node[iPoint]->GetBetaInc2(); +// Enthalpy_b = node[iPoint]->GetEnthalpy(); +// conv_numerics->GetInviscidIncProjFlux(&Density_b, Velocity_b, &V_reflected[0], &BetaInc2_b, &Enthalpy_b, Normal, Residual); +// +// /*--- Update residual value ---*/ +// LinSysRes.AddBlock(iPoint, Residual); +// +// /*--- Jacobian contribution for implicit integration. ---*/ +// if (implicit) { +// //Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); +// } +// +// } +// } +// +// /*--- Free locally allocated memory ---*/ +// delete [] Normal; +// delete [] UnitNormal; +// delete [] Velocity_b; +//} + void CIncEulerSolver::BC_Euler_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, unsigned short val_marker) { @@ -5505,7 +5647,7 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain for (iDim = 0; iDim < nDim; iDim++) V_infty[iDim+1] = GetVelocity_Inf(iDim); - /*--- Far-field pressure set to static pressure (0.0). ---*/ + /*--- Far-field pressure set to static pressure (0.0). ---*/ //TK:: why 0.0?! V_infty[0] = GetPressure_Inf(); From 4e68c6e14abe98740544d2d002a200235dd4ecad Mon Sep 17 00:00:00 2001 From: cvencro Date: Tue, 26 Mar 2019 07:43:56 +0000 Subject: [PATCH 10/31] Primitive additions to the Jacobian in convective schemes for grid movement --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index f20edb8a3a55..f29c869f9556 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -245,8 +245,12 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J /*--- Implicit terms ---*/ if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + 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; } } } @@ -395,8 +399,12 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ /*--- Implicit terms ---*/ if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + 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; } } } @@ -599,8 +607,12 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ /*--- Implicit terms ---*/ if (implicit) { - val_Jacobian_i[iVar][iVar] -= 0.5*ProjVelocity; - val_Jacobian_j[iVar][iVar] -= 0.5*ProjVelocity; + 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; } } } From 276a7103b14d7e2e3db141cbfe70eea43182a1fe Mon Sep 17 00:00:00 2001 From: cvencro Date: Mon, 20 May 2019 08:15:29 +0100 Subject: [PATCH 11/31] Enable unsteady adjoint with dynamic mesh movement --- Common/src/config_structure.cpp | 6 +++--- Common/src/geometry_structure.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 30af16415e35..41394dd34b62 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -3999,9 +3999,9 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Restart_Flow = false; - if (Grid_Movement) { - SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION); - } + //if (Grid_Movement) { + // SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION); + //} if (Unst_AdjointIter- long(nExtIter) < 0){ SU2_MPI::Error(string("Invalid iteration number requested for unsteady adjoint.\n" ) + diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp index 46b8c506b081..1c97dfa61e9d 100644 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -17246,7 +17246,7 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config) { if (compressible) { skipVar += skipMult*(nDim+2); } if (sst && !frozen_visc) { skipVar += skipMult*2;} if (sa && !frozen_visc) { skipVar += skipMult*1;} - if (grid_movement) { skipVar += nDim;} + //CVC: Debug: if (grid_movement) { skipVar += nDim;} } else if (Kind_Solver == DISC_ADJ_HEAT) { skipVar += 1; From e6ea7eb388111521f893b867e3c0062a461e6056 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Mon, 20 May 2019 14:17:38 +0200 Subject: [PATCH 12/31] Fixed adjoint with grid movement such that correct Coords and GridVels are in thier positions. --- Common/include/dual_grid_structure.hpp | 68 ++++++++++++++++++++++++-- Common/include/dual_grid_structure.inl | 41 ++++++++++++++++ Common/src/config_structure.cpp | 6 +-- Common/src/dual_grid_structure.cpp | 27 ++++++++-- Common/src/geometry_structure.cpp | 2 +- Common/src/grid_movement_structure.cpp | 10 ++-- SU2_CFD/include/solver_structure.hpp | 2 +- SU2_CFD/src/driver_structure.cpp | 2 +- SU2_CFD/src/iteration_structure.cpp | 61 ++++++++++++++++++----- SU2_CFD/src/solver_direct_mean.cpp | 5 +- SU2_CFD/src/solver_direct_mean_inc.cpp | 5 +- SU2_CFD/src/solver_structure.cpp | 8 +-- 12 files changed, 199 insertions(+), 38 deletions(-) mode change 100644 => 100755 Common/include/dual_grid_structure.hpp mode change 100644 => 100755 Common/include/dual_grid_structure.inl mode change 100644 => 100755 Common/src/config_structure.cpp mode change 100644 => 100755 Common/src/dual_grid_structure.cpp mode change 100644 => 100755 Common/src/geometry_structure.cpp mode change 100644 => 100755 Common/src/grid_movement_structure.cpp mode change 100644 => 100755 SU2_CFD/src/driver_structure.cpp mode change 100644 => 100755 SU2_CFD/src/iteration_structure.cpp diff --git a/Common/include/dual_grid_structure.hpp b/Common/include/dual_grid_structure.hpp old mode 100644 new mode 100755 index 2f6a3348d2ab..388ff70a9890 --- a/Common/include/dual_grid_structure.hpp +++ b/Common/include/dual_grid_structure.hpp @@ -159,6 +159,9 @@ class CPoint : public CDualGrid { *Coord_n1, /*!< \brief Coordinates at time n-1 for use with dynamic meshes. */ *Coord_p1; /*!< \brief Coordinates at time n+1 for use with dynamic meshes. */ su2double *GridVel; /*!< \brief Velocity of the grid for dynamic mesh cases. */ + su2double *GridVel_n; /*!< \brief Velocity of the grid for dynamic mesh cases of previous time step. */ + su2double *GridVel_n1; /*!< \brief Velocity of the grid for dynamic mesh cases of second to current time step (2nd order timestepping only). */ + su2double *GridVel_Old; /*!< \brief Velocity of the grid for dynamic mesh cases, intermediate container. */ su2double **GridVel_Grad; /*!< \brief Gradient of the grid velocity for dynamic meshes. */ unsigned long Parent_CV; /*!< \brief Index of the parent control volume in the agglomeration process. */ unsigned short nChildren_CV; /*!< \brief Number of children in the agglomeration process. */ @@ -281,7 +284,7 @@ class CPoint : public CDualGrid { * \return pointer to the coordinate of the point. */ su2double *GetCoord(void); - + /*! * \brief Set the coordinates for the control volume. * \param[in] val_dim - Position to store the coordinate. @@ -541,12 +544,12 @@ class CPoint : public CDualGrid { su2double* GetCoord_p1(void); /*! - * \brief Set the coordinates of the control volume at time n. + * \brief Set the coordinates of the control volume at time n to the ones in Coord. */ void SetCoord_n(void); /*! - * \brief Set the coordinates of the control volume at time n-1. + * \brief Set the coordinates of the control volume at time n-1 to the ones in Coord_n. */ void SetCoord_n1(void); @@ -659,6 +662,24 @@ class CPoint : public CDualGrid { * \return Grid velocity at the point. */ su2double *GetGridVel(void); + + /*! + * \brief Get the value of the grid velocity at the point from previous timestep. + * \return Grid velocity at the point. + */ + su2double *GetGridVel_n(void); + + /*! + * \brief Get the value of the grid velocity at the point from 2nd to current timestep (2nd order timestepping only). + * \return Grid velocity at the point. + */ + su2double *GetGridVel_n1(void); + + /*! + * \brief Get the value of the grid velocity at the point from helper container. + * \return Grid velocity at the point. + */ + su2double *GetGridVel_Old(void); /*! * \brief Get the value of the grid velocity gradient at the point. @@ -682,6 +703,11 @@ class CPoint : public CDualGrid { * \param[in] val_coord_old - Value of the coordinates. */ void SetCoord_Old(su2double *val_coord_old); + + /*! + * \brief Set the value of the vector Coord_Old to Coord. + */ + void SetCoord_Old(void); /*! * \brief Set the value of the grid velocity at the point. @@ -692,11 +718,45 @@ class CPoint : public CDualGrid { /*! * \overload + * \brief Set the value of the grid velocity at the point. * \param[in] val_gridvel - Value of the grid velocity. */ void SetGridVel(su2double *val_gridvel); - /*! + /*! + * \brief Set the value of the grid velocity at the point. + * \param[in] val_gridvel - value array of the grid velocities. + */ + void SetGridVel_Old(su2double *val_gridvel); + + /*! + * \brief Set the values of the grid velocity to the current ones GridVel. + */ + void SetGridVel_Old(void); + + /*! + * \brief Set the values of the grid velocity to the current ones GridVel. + */ + void SetGridVel_n(void); + + /*! + * \brief Set the value of the grid velocity at the point. + * \param[in] val_gridvel - value array of the grid velocities. + */ + void SetGridVel_n(su2double *val_gridvel); + + /*! + * \brief Set the values of the grid velocity to the ones in GridVel_n. + */ + void SetGridVel_n1(void); + + /*! + * \brief Set the value of the grid velocity at the point. + * \param[in] val_gridvel - value array of the grid velocities. + */ + void SetGridVel_n1(su2double *val_gridvel); + + /*! * \brief Set the gradient of the grid velocity. * \param[in] val_var - Index of the variable. * \param[in] val_dim - Index of the dimension. diff --git a/Common/include/dual_grid_structure.inl b/Common/include/dual_grid_structure.inl old mode 100644 new mode 100755 index ef9cc517199a..ed82143baf27 --- a/Common/include/dual_grid_structure.inl +++ b/Common/include/dual_grid_structure.inl @@ -110,6 +110,12 @@ inline su2double *CPoint::GetCoord_Sum(void) { return Coord_Sum; } inline su2double *CPoint::GetGridVel(void) { return GridVel; } +inline su2double *CPoint::GetGridVel_n(void) { return GridVel_n; } + +inline su2double *CPoint::GetGridVel_n1(void) { return GridVel_n1; } + +inline su2double *CPoint::GetGridVel_Old(void) { return GridVel_Old; } + inline su2double **CPoint::GetGridVel_Grad(void) { return GridVel_Grad; } inline void CPoint::SetCoord_Old(su2double *val_coord_old) { @@ -168,6 +174,36 @@ inline void CPoint::SetGridVel(su2double *val_gridvel) { GridVel[iDim] = val_gridvel[iDim]; } +inline void CPoint::SetGridVel_Old(su2double *val_gridvel) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_Old[iDim] = val_gridvel[iDim]; +} + +inline void CPoint::SetGridVel_Old(void) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_Old[iDim] = GridVel[iDim]; +} + +inline void CPoint::SetGridVel_n(void) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_n[iDim] = GridVel[iDim]; +} + +inline void CPoint::SetGridVel_n(su2double *val_gridvel) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_n[iDim] = val_gridvel[iDim]; +} + +inline void CPoint::SetGridVel_n1(void) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_n1[iDim] = GridVel_n[iDim]; +} + +inline void CPoint::SetGridVel_n1(su2double *val_gridvel) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + GridVel_n1[iDim] = val_gridvel[iDim]; +} + inline void CPoint::SetVolume_n (void) { Volume[1] = Volume[0]; } inline void CPoint::SetVolume_nM1 (void) { Volume[2] = Volume[1]; } @@ -176,6 +212,11 @@ inline su2double CPoint::GetVolume_n (void) { return Volume[1]; } inline su2double CPoint::GetVolume_nM1 (void) { return Volume[2]; } +inline void CPoint::SetCoord_Old (void) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Coord_Old[iDim] = Coord[iDim]; +} + inline void CPoint::SetCoord_n (void) { for (unsigned short iDim = 0; iDim < nDim; iDim++) Coord_n[iDim] = Coord[iDim]; diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp old mode 100644 new mode 100755 index 37be2c450bd1..d9cb993c09b9 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -4001,8 +4001,8 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ if (Unsteady_Simulation) { Restart_Flow = false; - - if (Grid_Movement) { + if(false) { + //if (Grid_Movement) { SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION); } @@ -7305,7 +7305,7 @@ string CConfig::GetUnsteady_FileName(string val_filename, int val_iter) { char buffer[50]; /*--- Check that a positive value iteration is requested (for now). ---*/ - + if(rank==MASTER_NODE) {cout << "CConfig::GetUnsteady_FileName val_iter: " << val_iter << endl;} if (val_iter < 0) { SU2_MPI::Error("Requesting a negative iteration number for the restart file!!", CURRENT_FUNCTION); } diff --git a/Common/src/dual_grid_structure.cpp b/Common/src/dual_grid_structure.cpp old mode 100644 new mode 100755 index 582497a30774..1c4adfb282a6 --- a/Common/src/dual_grid_structure.cpp +++ b/Common/src/dual_grid_structure.cpp @@ -55,7 +55,8 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; + GridVel_n = NULL; GridVel_n1 = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -107,6 +108,9 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * if ( config->GetGrid_Movement() ) { GridVel = new su2double[nDim]; + GridVel_Old = new su2double[nDim]; + GridVel_n = new su2double[nDim]; + GridVel_n1 = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim++) GridVel[iDim] = 0.0; @@ -126,6 +130,7 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; + Coord_Old = new su2double[nDim]; } } @@ -146,7 +151,8 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; + GridVel_n = NULL; GridVel_n1 = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -198,6 +204,10 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g /*--- Storage of grid velocities for dynamic meshes ---*/ if ( config->GetGrid_Movement() ) { GridVel = new su2double[nDim]; + GridVel_Old = new su2double[nDim]; + GridVel_n = new su2double[nDim]; + GridVel_n1 = new su2double[nDim]; + for (iDim = 0; iDim < nDim; iDim++) GridVel[iDim] = 0.0; @@ -215,6 +225,7 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; + Coord_Old = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim ++) { Coord_p1[iDim] = Coord[iDim]; Coord_n[iDim] = Coord[iDim]; @@ -240,7 +251,9 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; + GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; + GridVel_n = NULL; GridVel_n1 = NULL; + /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ if ( config->GetUnsteady_Simulation() == NO ) { @@ -294,6 +307,10 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord if (config->GetGrid_Movement()) { GridVel = new su2double[nDim]; + GridVel_Old = new su2double[nDim]; + GridVel_n = new su2double[nDim]; + GridVel_n1 = new su2double[nDim]; + for (iDim = 0; iDim < nDim; iDim ++) GridVel[iDim] = 0.0; @@ -311,6 +328,7 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; + Coord_Old = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim ++) { Coord_p1[iDim] = Coord[iDim]; Coord_n[iDim] = Coord[iDim]; @@ -335,6 +353,9 @@ CPoint::~CPoint() { if (Coord_n1 != NULL) delete[] Coord_n1; if (Coord_p1 != NULL) delete[] Coord_p1; if (GridVel != NULL) delete[] GridVel; + if (GridVel_Old != NULL) delete[] GridVel_Old; + if (GridVel_n != NULL) delete[] GridVel_n; + if (GridVel_n1 != NULL) delete[] GridVel_n1; if (GridVel_Grad != NULL) { for (unsigned short iDim = 0; iDim < nDim; iDim++) delete [] GridVel_Grad[iDim]; diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp old mode 100644 new mode 100755 index f7a34f24f9e8..09a4eb491471 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -18522,7 +18522,7 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config) { if (compressible) { skipVar += skipMult*(nDim+2); } if (sst && !frozen_visc) { skipVar += skipMult*2;} if (sa && !frozen_visc) { skipVar += skipMult*1;} - if (grid_movement) { skipVar += nDim;} + //if (grid_movement) { skipVar += nDim;} } else if (Kind_Solver == DISC_ADJ_HEAT) { skipVar += 1; diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp old mode 100644 new mode 100755 index f0e481359ba5..5341a26df703 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -1897,7 +1897,7 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry *geometry, CConfig *config, su2double dtheta, dphi, dpsi, cosTheta, sinTheta; su2double cosPhi, sinPhi, cosPsi, sinPsi; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = config->GetContinuous_Adjoint(); + bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); /*--- Problem dimension and physical time step ---*/ @@ -2067,7 +2067,7 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry *geometry, CConfig *config, u unsigned short nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = config->GetContinuous_Adjoint(); + bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); /*--- Retrieve values from the config file ---*/ @@ -2223,7 +2223,7 @@ void CVolumetricMovement::Rigid_Plunging(CGeometry *geometry, CConfig *config, u unsigned short iDim, nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = config->GetContinuous_Adjoint(); + bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); /*--- Retrieve values from the config file ---*/ @@ -2362,7 +2362,7 @@ void CVolumetricMovement::Rigid_Translation(CGeometry *geometry, CConfig *config unsigned short iDim, nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = config->GetContinuous_Adjoint(); + bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); // This should be consistent over all Rigid_* routines /*--- Retrieve values from the config file ---*/ @@ -6463,7 +6463,7 @@ void CSurfaceMovement::SetExternal_Deformation(CGeometry *geometry, CConfig *con string DV_Filename, UnstExt, text_line; ifstream surface_positions; bool unsteady = config->GetUnsteady_Simulation(); - bool adjoint = config->GetContinuous_Adjoint(); + bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); /*--- Load stuff from config ---*/ diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp index 3be00e84ef44..67bc8bea41a5 100644 --- a/SU2_CFD/include/solver_structure.hpp +++ b/SU2_CFD/include/solver_structure.hpp @@ -603,7 +603,7 @@ class CSolver { * \brief Load the geometries at the previous time states n and nM1. * \param[in] geometry - Geometrical definition of the problem. */ - void Restart_OldGeometry(CGeometry *geometry, CConfig *config); + void Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_iter); /*! * \brief A virtual member. diff --git a/SU2_CFD/src/driver_structure.cpp b/SU2_CFD/src/driver_structure.cpp old mode 100644 new mode 100755 index 7f6f0e00f2bf..646050aa851e --- a/SU2_CFD/src/driver_structure.cpp +++ b/SU2_CFD/src/driver_structure.cpp @@ -3777,7 +3777,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/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp old mode 100644 new mode 100755 index 6f4d014bdb9e..1a63000021c2 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2086,16 +2086,16 @@ CDiscAdjFluidIteration::CDiscAdjFluidIteration(CConfig *config) : CIteration(con CDiscAdjFluidIteration::~CDiscAdjFluidIteration(void) { } void CDiscAdjFluidIteration::Preprocess(COutput *output, - CIntegration ****integration_container, - CGeometry ****geometry_container, - CSolver *****solver_container, - CNumerics ******numerics_container, - CConfig **config_container, - CSurfaceMovement **surface_movement, - CVolumetricMovement ***grid_movement, - CFreeFormDefBox*** FFDBox, - unsigned short val_iZone, - unsigned short val_iInst) { + CIntegration ****integration_container, + CGeometry ****geometry_container, + CSolver *****solver_container, + CNumerics ******numerics_container, + CConfig **config_container, + CSurfaceMovement **surface_movement, + CVolumetricMovement ***grid_movement, + CFreeFormDefBox*** FFDBox, + unsigned short val_iZone, + unsigned short val_iInst) { unsigned long IntIter = 0, iPoint; config_container[ZONE_0]->SetIntIter(IntIter); @@ -2106,6 +2106,7 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, unsigned short iMesh; int Direct_Iter; bool heat = config_container[val_iZone]->GetWeakly_Coupled_Heat(); + bool grid_movement_bool = config_container[val_iZone]->GetGrid_Movement(); /*--- For the unsteady adjoint, load direct solutions from restart files. ---*/ @@ -2124,7 +2125,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if (dual_time_2nd) { /*--- Load solution at timestep n-2 ---*/ - LoadUnsteady_Solution(geometry_container, solver_container,config_container, val_iZone, val_iInst, Direct_Iter-2); /*--- Push solution back to correct array ---*/ @@ -2133,6 +2133,13 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(); solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n1(); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n1(); + + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n1(); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(); solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n1(); @@ -2147,7 +2154,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if (dual_time) { /*--- Load solution at timestep n-1 ---*/ - LoadUnsteady_Solution(geometry_container, solver_container,config_container, val_iZone, val_iInst, Direct_Iter-1); /*--- Push solution back to correct array ---*/ @@ -2155,6 +2161,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(); } @@ -2166,7 +2176,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, } /*--- Load solution timestep n ---*/ - LoadUnsteady_Solution(geometry_container, solver_container,config_container, val_iInst, val_iZone, Direct_Iter); } @@ -2174,6 +2183,12 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if ((ExtIter > 0) && dual_time){ + /*--- + Here the primal solutions (only working variables) are loaded and put in the correct order + into containers. For ALE the mesh coordinates and the Grid-Velocities have to be put into the + correct containers as well, i.e. follow the same logic for the solution. + ---*/ + /*--- Load solution timestep n-1 | n-2 for DualTimestepping 1st | 2nd order ---*/ if (dual_time_1st){ LoadUnsteady_Solution(geometry_container, solver_container,config_container, val_iInst, val_iZone, Direct_Iter - 1); @@ -2187,6 +2202,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_OldSolution(); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_Old(); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_Old(); + } if (turbulent){ solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_OldSolution(); } @@ -2201,6 +2220,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->SetSolution(solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_time_n()); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_n()); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_n()); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->SetSolution(solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_time_n()); } @@ -2214,6 +2237,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_Old()); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_Old()); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_Old()); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_Old()); } @@ -2228,6 +2255,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_time_n1()); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_n1()); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_n1()); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_time_n1()); } @@ -2240,6 +2271,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config_container[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n1(solver_container[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_Old()); + if (grid_movement_bool) { + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n1(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_Old()); + geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n1(geometry_container[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_Old()); + } if (turbulent) { solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n1(solver_container[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_Old()); } diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index 865a49d71b14..d53f9372ceaf 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -13917,8 +13917,9 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy ---*/ - if (dual_time && grid_movement) - Restart_OldGeometry(geometry[MESH_0], config); + //if (dual_time && grid_movement) + if (dual_time && grid_movement && false) + Restart_OldGeometry(geometry[MESH_0], config, val_iter); delete [] Coord; diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index f28b8acc0ce8..181194832f6c 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -7522,8 +7522,9 @@ 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 && grid_movement) - Restart_OldGeometry(geometry[MESH_0], config); + //if (dual_time && grid_movement) + if (dual_time && grid_movement && false) + Restart_OldGeometry(geometry[MESH_0], config, val_iter); delete [] Coord; diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp index 97fc32ed4bd2..201cf3cf8fe3 100644 --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -2015,7 +2015,7 @@ void CSolver::SolveTypicalSectionWingModel(CGeometry *geometry, su2double Cl, su } -void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { +void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_iter) { /*--- This function is intended for dual time simulations ---*/ @@ -2048,7 +2048,8 @@ void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { /*-------------------------------------------------------------------------------------------*/ /*--- Modify file name for an unsteady restart ---*/ - Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-1; + //Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-1; + Unst_RestartIter = val_iter - 1; filename_n = config->GetUnsteady_FileName(filename, Unst_RestartIter); /*--- Open the restart file, throw an error if this fails. ---*/ @@ -2119,7 +2120,8 @@ void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { string filename_n1; /*--- Modify file name for an unsteady restart ---*/ - Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-2; + //Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-2; + Unst_RestartIter = val_iter - 2; filename_n1 = config->GetUnsteady_FileName(filename, Unst_RestartIter); /*--- Open the restart file, throw an error if this fails. ---*/ From 6a678f3ad0e6e178b64271b46bc21471d457becb Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Tue, 21 May 2019 10:36:18 +0200 Subject: [PATCH 13/31] Fixing failing Regression test by reverting Restart_OldGeometry to original. --- SU2_CFD/include/solver_structure.hpp | 3 ++- SU2_CFD/src/solver_direct_mean.cpp | 5 ++--- SU2_CFD/src/solver_direct_mean_inc.cpp | 5 ++--- SU2_CFD/src/solver_structure.cpp | 8 +++----- 4 files changed, 9 insertions(+), 12 deletions(-) mode change 100644 => 100755 SU2_CFD/include/solver_structure.hpp mode change 100644 => 100755 SU2_CFD/src/solver_structure.cpp diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp old mode 100644 new mode 100755 index 8cae3ba015b0..e4ced990ccee --- a/SU2_CFD/include/solver_structure.hpp +++ b/SU2_CFD/include/solver_structure.hpp @@ -550,8 +550,9 @@ 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. */ - void Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_iter); + void Restart_OldGeometry(CGeometry *geometry, CConfig *config); /*! * \brief A virtual member. diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index 3fa1d0b4de28..5758da604d28 100755 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -13085,9 +13085,8 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy ---*/ - //if (dual_time && grid_movement) - if (dual_time && grid_movement && false) - Restart_OldGeometry(geometry[MESH_0], config, val_iter); + if (dual_time && grid_movement && (config->GetKind_GridMovement() != RIGID_MOTION)) + Restart_OldGeometry(geometry[MESH_0], config); delete [] Coord; diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index c13939b78361..c30538757824 100755 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -6707,9 +6707,8 @@ 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 && grid_movement) - if (dual_time && grid_movement && false) - Restart_OldGeometry(geometry[MESH_0], config, val_iter); + if (dual_time && grid_movement && (config->GetKind_GridMovement() != RIGID_MOTION)) + Restart_OldGeometry(geometry[MESH_0], config); delete [] Coord; diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp old mode 100644 new mode 100755 index 7151743bf940..d0b6c9e3fe8f --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -3927,7 +3927,7 @@ void CSolver::SolveTypicalSectionWingModel(CGeometry *geometry, su2double Cl, su } -void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_iter) { +void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { /*--- This function is intended for dual time simulations ---*/ @@ -3960,8 +3960,7 @@ void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_ /*-------------------------------------------------------------------------------------------*/ /*--- Modify file name for an unsteady restart ---*/ - //Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-1; - Unst_RestartIter = val_iter - 1; + Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-1; filename_n = config->GetUnsteady_FileName(filename, Unst_RestartIter); /*--- Open the restart file, throw an error if this fails. ---*/ @@ -4032,8 +4031,7 @@ void CSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config, int val_ string filename_n1; /*--- Modify file name for an unsteady restart ---*/ - //Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-2; - Unst_RestartIter = val_iter - 2; + Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_RestartIter())-2; filename_n1 = config->GetUnsteady_FileName(filename, Unst_RestartIter); /*--- Open the restart file, throw an error if this fails. ---*/ From 8e1c3be29861e3c769977166f545696a6f4d6120 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Fri, 7 Jun 2019 10:45:12 +0200 Subject: [PATCH 14/31] FDS Preaccumulation with GridVels. --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 4 ++++ SU2_CFD/src/output_su2.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) mode change 100644 => 100755 SU2_CFD/src/output_su2.cpp diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 554a0c88f5fa..4b85f66c4972 100755 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -90,6 +90,10 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); AD::SetPreaccIn(Normal, nDim); + if (grid_movement) { + AD::SetPreaccIn(GridVel_i, nDim); + AD::SetPreaccIn(GridVel_j, nDim); + } /*--- Face area (norm or the normal vector) ---*/ diff --git a/SU2_CFD/src/output_su2.cpp b/SU2_CFD/src/output_su2.cpp old mode 100644 new mode 100755 index f3994e9d1900..fe67f781fc8c --- a/SU2_CFD/src/output_su2.cpp +++ b/SU2_CFD/src/output_su2.cpp @@ -224,7 +224,7 @@ void COutput::SetSU2_MeshASCII(CConfig *config, CGeometry *geometry, unsigned sh /*--- Get the total number of periodic transformations ---*/ nPeriodic = config->GetnPeriodicIndex(); - output_file << "NPERIODIC= " << nPeriodic << endl; + //output_file << "NPERIODIC= " << nPeriodic << endl; /*--- From iPeriodic obtain the iMarker ---*/ From 97577f0c858172536569e0bd96c5b80baac68629 Mon Sep 17 00:00:00 2001 From: cvencro Date: Thu, 8 Aug 2019 00:47:39 +0100 Subject: [PATCH 15/31] Fix config and rotation rate function names --- SU2_CFD/src/iteration_structure.cpp | 2 +- SU2_CFD/src/numerics_direct_mean_inc.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 32b560bf83b3..9c2d138e7d60 100755 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2127,7 +2127,7 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, unsigned short iMesh; int Direct_Iter; bool heat = config[val_iZone]->GetWeakly_Coupled_Heat(); - bool grid_movement_bool = config_container[val_iZone]->GetGrid_Movement(); + bool grid_movement_bool = config[val_iZone]->GetGrid_Movement(); /*--- Read the target pressure for inverse design. ---------------------------------------------*/ if (config[val_iZone]->GetInvDesign_Cp() == YES) diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 4b85f66c4972..9bd365a82747 100755 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -1000,9 +1000,9 @@ void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2d /*--- Retrieve the angular velocity vector from config. ---*/ - Omega[0] = config->GetRotation_Rate_X(config->GetiZone())/config->GetOmega_Ref(); - Omega[1] = config->GetRotation_Rate_Y(config->GetiZone())/config->GetOmega_Ref(); - Omega[2] = config->GetRotation_Rate_Z(config->GetiZone())/config->GetOmega_Ref(); + for (iDim = 0; iDim < 3; iDim++){ + Omega[iDim] = config->GetRotation_Rate(iDim)/config->GetOmega_Ref(); + } /*--- Primitive variables at point i and j ---*/ From 8c9c681699c5b36cb7fa5e5848bec41bdef065c6 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 21 Aug 2019 22:19:41 +0100 Subject: [PATCH 16/31] tidy comments --- SU2_CFD/include/numerics_structure.hpp | 1 - SU2_CFD/src/solver_direct_mean_inc.cpp | 144 +------------------------ 2 files changed, 1 insertion(+), 144 deletions(-) diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index 5429ca9ff696..6b62fb0157ac 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5247,7 +5247,6 @@ class CSourceIncBodyForce : public CNumerics { * \class CSourceIncRotatingFrame_Flow * \brief Class for a rotating frame source term. * \ingroup SourceDiscr - * \author F. Palacios, T. Economon, C. Venkatesan-Crome */ class CSourceIncRotatingFrame_Flow : public CNumerics { public: diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 5ae00df14a92..07d94540b941 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -564,22 +564,6 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned /*--- Initialize the solution to the far-field state everywhere. ---*/ - //TK:: Disturb the initial solution slightly to see whether ALE stuff is doing anything at all - //if(rank == MASTER_NODE) cout << "Disturbing initial solution slightly, CIncEulerSolver::CIncEulerSolver" << endl; - //su2double tmpPressure_Inf = Pressure_Inf + 0.2*Pressure_Inf; - //su2double *tmpVelocity_Inf = new su2double[nDim]; - //for (iDim = 0; iDim < nDim; iDim++) { - // if (Velocity_Inf[iDim] == 0.0) { - // tmpVelocity_Inf[iDim] = 12.3; - // } else { - // tmpVelocity_Inf[iDim] = Velocity_Inf[iDim] + 0.2*Velocity_Inf[iDim]; - // } - //} - //su2double tmpTemperature_Inf = Temperature_Inf + 0.2*Temperature_Inf; - //TK:: end - - //for (iPoint = 0; iPoint < nPoint; iPoint++) - // node[iPoint] = new CIncEulerVariable(tmpPressure_Inf, tmpVelocity_Inf, tmpTemperature_Inf, nDim, nVar, config); for (iPoint = 0; iPoint < nPoint; iPoint++) node[iPoint] = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nDim, nVar, config); @@ -4619,132 +4603,6 @@ void CIncEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { } -//void CIncEulerSolver::BC_Euler_Wall(CGeometry *geometry, -// CSolver **solver_container, -// CNumerics *conv_numerics, -// CConfig *config, -// unsigned short val_marker) { -// -// unsigned short iDim, iVar; -// unsigned long iVertex, iPoint; -// -// bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); -// bool grid_movement = config->GetGrid_Movement(); -// -// /*--- Allocation of variables necessary for convective fluxes. ---*/ -// su2double Area, ProjVelocity_i; -// su2double *V_reflected, *V_domain; -// su2double *Normal = new su2double[nDim]; -// su2double *UnitNormal = new su2double[nDim]; -// -// su2double *GridVel = NULL; -// su2double ProjGridVel = 0.0; -// su2double *Velocity_b; -// Velocity_b = new su2double[nDim]; -// su2double Density_b, BetaInc2_b, Enthalpy_b; -// -// /*--- Loop over all the vertices on this boundary marker. ---*/ -// for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { -// -// iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); -// -// /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ -// if (geometry->node[iPoint]->GetDomain()) { -// -// /*-------------------------------------------------------------------------------*/ -// /*--- Step 1: For the convective fluxes, create a reflected state of the ---*/ -// /*--- Primitive variables by copying all interior values to the ---*/ -// /*--- reflected. Only the velocity is mirrored along the symmetry ---*/ -// /*--- axis. Based on the Upwind_Residual routine. ---*/ -// /*-------------------------------------------------------------------------------*/ -// -// /*--- Normal vector for a random vertex (zero) on this marker (negate for outward convention). ---*/ -// geometry->vertex[val_marker][0]->GetNormal(Normal); -// for (iDim = 0; iDim < nDim; iDim++) -// Normal[iDim] = -Normal[iDim]; -// -// /*--- Compute unit normal, to be used for unit tangential, projected velocity and velocity component gradients. ---*/ -// Area = 0.0; -// for (iDim = 0; iDim < nDim; iDim++) -// Area += Normal[iDim]*Normal[iDim]; -// Area = sqrt (Area); -// -// for (iDim = 0; iDim < nDim; iDim++) -// UnitNormal[iDim] = -Normal[iDim]/Area; -// -// /*--- Allocate the reflected state at the symmetry boundary. ---*/ -// V_reflected = GetCharacPrimVar(val_marker, iVertex); -// -// /*--- Grid movement ---*/ -// //if (grid_movement) -// // conv_numerics->SetGridVel(geometry->node[iPoint]->GetGridVel(), geometry->node[iPoint]->GetGridVel()); -// -// /*--- Normal vector for this vertex (negate for outward convention). ---*/ -// geometry->vertex[val_marker][iVertex]->GetNormal(Normal); -// for (iDim = 0; iDim < nDim; iDim++) -// Normal[iDim] = -Normal[iDim]; -// //conv_numerics->SetNormal(Normal); -// -// /*--- Get current solution at this boundary node ---*/ -// V_domain = node[iPoint]->GetPrimitive(); -// -// /*--- Set the reflected state based on the boundary node. Scalars are copied and -// the velocity is mirrored along the symmetry boundary, i.e. the velocity in -// normal direction is substracted twice. ---*/ -// for(iVar = 0; iVar < nPrimVar; iVar++) -// V_reflected[iVar] = node[iPoint]->GetPrimitive(iVar); -// -// /*--- Compute velocity in normal direction (ProjVelcity_i=(v*n)) und substract twice from -// velocity in normal direction: v_r = v - 2 (v*n)n ---*/ -// ProjVelocity_i = 0.0; -// for (iDim = 0; iDim < nDim; iDim++) -// ProjVelocity_i += node[iPoint]->GetVelocity(iDim)*UnitNormal[iDim]; -// -// for (iDim = 0; iDim < nDim; iDim++) -// V_reflected[iDim+1] = node[iPoint]->GetVelocity(iDim) - ProjVelocity_i*UnitNormal[iDim]; -// -// -// if (grid_movement) { -// GridVel = geometry->node[iPoint]->GetGridVel(); -// ProjGridVel = 0.0; -// for (iDim = 0; iDim < nDim; iDim++) ProjGridVel += GridVel[iDim]*UnitNormal[iDim]; -// for (iDim = 0; iDim < nDim; iDim++) V_reflected[iDim+1] += GridVel[iDim] - ProjGridVel * UnitNormal[iDim]; -// } -// -// for (iDim = 0; iDim < nDim; iDim++) -// Velocity_b[iDim] = V_reflected[iDim+1]; -// -// /*--- Set Primitive and Secondary for numerics class. ---*/ -// //conv_numerics->SetPrimitive(V_domain, V_reflected); -// //conv_numerics->SetSecondary(node[iPoint]->GetSecondary(), node[iPoint]->GetSecondary()); -// -// /*--- Compute the residual using an upwind scheme. ---*/ -// //conv_numerics->ComputeResidual(Residual, Jacobian_i, Jacobian_j, config); -// -// //conv_numerics->GetInviscidProjFlux(&Density_b, Velocity_b, &Pressure_b, &Enthalpy_b, NormalArea, Residual); -// -// Density_b = node[iPoint]->GetDensity(); -// BetaInc2_b = node[iPoint]->GetBetaInc2(); -// Enthalpy_b = node[iPoint]->GetEnthalpy(); -// conv_numerics->GetInviscidIncProjFlux(&Density_b, Velocity_b, &V_reflected[0], &BetaInc2_b, &Enthalpy_b, Normal, Residual); -// -// /*--- Update residual value ---*/ -// LinSysRes.AddBlock(iPoint, Residual); -// -// /*--- Jacobian contribution for implicit integration. ---*/ -// if (implicit) { -// //Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); -// } -// -// } -// } -// -// /*--- Free locally allocated memory ---*/ -// delete [] Normal; -// delete [] UnitNormal; -// delete [] Velocity_b; -//} - void CIncEulerSolver::BC_Euler_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *numerics, CConfig *config, unsigned short val_marker) { @@ -4874,7 +4732,7 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain for (iDim = 0; iDim < nDim; iDim++) V_infty[iDim+1] = GetVelocity_Inf(iDim); - /*--- Far-field pressure set to static pressure (0.0). ---*/ //TK:: why 0.0?! + /*--- Far-field pressure set to static pressure (0.0). ---*/ V_infty[0] = GetPressure_Inf(); From 4c4f16302426224e99e70cc9f8ae2c3a06db6583 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 21 Aug 2019 22:21:07 +0100 Subject: [PATCH 17/31] revert travis.yml --- .travis.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 341ca8a29623..428f6d61f7fb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,11 +18,11 @@ compiler: notifications: email: recipients: - - charanya.crome09@ic.ac.uk + - su2code-dev@lists.stanford.edu branches: only: - - feature_incompressible_ale + - develop virtualenv: system_site_packages: true From c12d70b848123f73bc36ca3d61cc46c580d61966 Mon Sep 17 00:00:00 2001 From: cvencro Date: Wed, 21 Aug 2019 23:32:15 +0100 Subject: [PATCH 18/31] revert changes for dynamic discrete adjoint to prepare for PR --- Common/include/dual_grid_structure.hpp | 68 ++------------------------ Common/include/dual_grid_structure.inl | 41 ---------------- Common/src/config_structure.cpp | 12 ++--- Common/src/dual_grid_structure.cpp | 27 ++-------- Common/src/geometry_structure.cpp | 2 +- Common/src/grid_movement_structure.cpp | 10 ++-- SU2_CFD/src/iteration_structure.cpp | 40 +-------------- SU2_CFD/src/solver_direct_mean.cpp | 2 +- 8 files changed, 19 insertions(+), 183 deletions(-) diff --git a/Common/include/dual_grid_structure.hpp b/Common/include/dual_grid_structure.hpp index eea54c9306f4..4c984f313714 100755 --- a/Common/include/dual_grid_structure.hpp +++ b/Common/include/dual_grid_structure.hpp @@ -161,9 +161,6 @@ class CPoint : public CDualGrid { *Coord_n1, /*!< \brief Coordinates at time n-1 for use with dynamic meshes. */ *Coord_p1; /*!< \brief Coordinates at time n+1 for use with dynamic meshes. */ su2double *GridVel; /*!< \brief Velocity of the grid for dynamic mesh cases. */ - su2double *GridVel_n; /*!< \brief Velocity of the grid for dynamic mesh cases of previous time step. */ - su2double *GridVel_n1; /*!< \brief Velocity of the grid for dynamic mesh cases of second to current time step (2nd order timestepping only). */ - su2double *GridVel_Old; /*!< \brief Velocity of the grid for dynamic mesh cases, intermediate container. */ su2double **GridVel_Grad; /*!< \brief Gradient of the grid velocity for dynamic meshes. */ unsigned long Parent_CV; /*!< \brief Index of the parent control volume in the agglomeration process. */ unsigned short nChildren_CV; /*!< \brief Number of children in the agglomeration process. */ @@ -286,7 +283,7 @@ class CPoint : public CDualGrid { * \return pointer to the coordinate of the point. */ su2double *GetCoord(void); - + /*! * \brief Set the coordinates for the control volume. * \param[in] val_dim - Position to store the coordinate. @@ -570,12 +567,12 @@ class CPoint : public CDualGrid { su2double* GetCoord_p1(void); /*! - * \brief Set the coordinates of the control volume at time n to the ones in Coord. + * \brief Set the coordinates of the control volume at time n. */ void SetCoord_n(void); /*! - * \brief Set the coordinates of the control volume at time n-1 to the ones in Coord_n. + * \brief Set the coordinates of the control volume at time n-1. */ void SetCoord_n1(void); @@ -688,24 +685,6 @@ class CPoint : public CDualGrid { * \return Grid velocity at the point. */ su2double *GetGridVel(void); - - /*! - * \brief Get the value of the grid velocity at the point from previous timestep. - * \return Grid velocity at the point. - */ - su2double *GetGridVel_n(void); - - /*! - * \brief Get the value of the grid velocity at the point from 2nd to current timestep (2nd order timestepping only). - * \return Grid velocity at the point. - */ - su2double *GetGridVel_n1(void); - - /*! - * \brief Get the value of the grid velocity at the point from helper container. - * \return Grid velocity at the point. - */ - su2double *GetGridVel_Old(void); /*! * \brief Get the value of the grid velocity gradient at the point. @@ -729,11 +708,6 @@ class CPoint : public CDualGrid { * \param[in] val_coord_old - Value of the coordinates. */ void SetCoord_Old(su2double *val_coord_old); - - /*! - * \brief Set the value of the vector Coord_Old to Coord. - */ - void SetCoord_Old(void); /*! * \brief Set the value of the grid velocity at the point. @@ -744,45 +718,11 @@ class CPoint : public CDualGrid { /*! * \overload - * \brief Set the value of the grid velocity at the point. * \param[in] val_gridvel - Value of the grid velocity. */ void SetGridVel(su2double *val_gridvel); - /*! - * \brief Set the value of the grid velocity at the point. - * \param[in] val_gridvel - value array of the grid velocities. - */ - void SetGridVel_Old(su2double *val_gridvel); - - /*! - * \brief Set the values of the grid velocity to the current ones GridVel. - */ - void SetGridVel_Old(void); - - /*! - * \brief Set the values of the grid velocity to the current ones GridVel. - */ - void SetGridVel_n(void); - - /*! - * \brief Set the value of the grid velocity at the point. - * \param[in] val_gridvel - value array of the grid velocities. - */ - void SetGridVel_n(su2double *val_gridvel); - - /*! - * \brief Set the values of the grid velocity to the ones in GridVel_n. - */ - void SetGridVel_n1(void); - - /*! - * \brief Set the value of the grid velocity at the point. - * \param[in] val_gridvel - value array of the grid velocities. - */ - void SetGridVel_n1(su2double *val_gridvel); - - /*! + /*! * \brief Set the gradient of the grid velocity. * \param[in] val_var - Index of the variable. * \param[in] val_dim - Index of the dimension. diff --git a/Common/include/dual_grid_structure.inl b/Common/include/dual_grid_structure.inl index 877058e7552b..48c6103a6155 100755 --- a/Common/include/dual_grid_structure.inl +++ b/Common/include/dual_grid_structure.inl @@ -118,12 +118,6 @@ inline su2double *CPoint::GetCoord_Sum(void) { return Coord_Sum; } inline su2double *CPoint::GetGridVel(void) { return GridVel; } -inline su2double *CPoint::GetGridVel_n(void) { return GridVel_n; } - -inline su2double *CPoint::GetGridVel_n1(void) { return GridVel_n1; } - -inline su2double *CPoint::GetGridVel_Old(void) { return GridVel_Old; } - inline su2double **CPoint::GetGridVel_Grad(void) { return GridVel_Grad; } inline void CPoint::SetCoord_Old(su2double *val_coord_old) { @@ -182,36 +176,6 @@ inline void CPoint::SetGridVel(su2double *val_gridvel) { GridVel[iDim] = val_gridvel[iDim]; } -inline void CPoint::SetGridVel_Old(su2double *val_gridvel) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_Old[iDim] = val_gridvel[iDim]; -} - -inline void CPoint::SetGridVel_Old(void) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_Old[iDim] = GridVel[iDim]; -} - -inline void CPoint::SetGridVel_n(void) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_n[iDim] = GridVel[iDim]; -} - -inline void CPoint::SetGridVel_n(su2double *val_gridvel) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_n[iDim] = val_gridvel[iDim]; -} - -inline void CPoint::SetGridVel_n1(void) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_n1[iDim] = GridVel_n[iDim]; -} - -inline void CPoint::SetGridVel_n1(su2double *val_gridvel) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - GridVel_n1[iDim] = val_gridvel[iDim]; -} - inline void CPoint::SetVolume_n (void) { Volume[1] = Volume[0]; } inline void CPoint::SetVolume_nM1 (void) { Volume[2] = Volume[1]; } @@ -220,11 +184,6 @@ inline su2double CPoint::GetVolume_n (void) { return Volume[1]; } inline su2double CPoint::GetVolume_nM1 (void) { return Volume[2]; } -inline void CPoint::SetCoord_Old (void) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Coord_Old[iDim] = Coord[iDim]; -} - inline void CPoint::SetCoord_n (void) { for (unsigned short iDim = 0; iDim < nDim; iDim++) Coord_n[iDim] = Coord[iDim]; diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 7f0a76b81169..e9de6dcc23ca 100755 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -4041,9 +4041,9 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Restart_Flow = false; - //if (GetGrid_Movement()) { - // SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION); - //} + if (GetGrid_Movement()) { + SU2_MPI::Error("Dynamic mesh movement currently not supported for the discrete adjoint solver.", CURRENT_FUNCTION); + } if (Unst_AdjointIter- long(nExtIter) < 0){ SU2_MPI::Error(string("Invalid iteration number requested for unsteady adjoint.\n" ) + @@ -4261,12 +4261,6 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } } - /*--- Grid motion is not yet supported with the incompressible solver. ---*/ - - /*if ((Kind_Regime == INCOMPRESSIBLE) && (GetGrid_Movement())) { - SU2_MPI::Error("Support for grid movement not yet implemented for incompressible flows.", CURRENT_FUNCTION); - }*/ - /*--- Assert that there are two markers being analyzed if the pressure drop objective function is selected. ---*/ diff --git a/Common/src/dual_grid_structure.cpp b/Common/src/dual_grid_structure.cpp index 907d31a28c4f..c550c58bca82 100755 --- a/Common/src/dual_grid_structure.cpp +++ b/Common/src/dual_grid_structure.cpp @@ -55,8 +55,7 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; - GridVel_n = NULL; GridVel_n1 = NULL; + GridVel = NULL; GridVel_Grad = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -108,9 +107,6 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * if ( config->GetGrid_Movement() ) { GridVel = new su2double[nDim]; - GridVel_Old = new su2double[nDim]; - GridVel_n = new su2double[nDim]; - GridVel_n1 = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim++) GridVel[iDim] = 0.0; @@ -130,7 +126,6 @@ CPoint::CPoint(unsigned short val_nDim, unsigned long val_globalindex, CConfig * Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; - Coord_Old = new su2double[nDim]; } } @@ -154,8 +149,7 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; - GridVel_n = NULL; GridVel_n1 = NULL; + GridVel = NULL; GridVel_Grad = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ @@ -208,10 +202,6 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g /*--- Storage of grid velocities for dynamic meshes ---*/ if ( config->GetGrid_Movement() ) { GridVel = new su2double[nDim]; - GridVel_Old = new su2double[nDim]; - GridVel_n = new su2double[nDim]; - GridVel_n1 = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) GridVel[iDim] = 0.0; @@ -229,7 +219,6 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, unsigned long val_g Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; - Coord_Old = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim ++) { Coord_p1[iDim] = Coord[iDim]; Coord_n[iDim] = Coord[iDim]; @@ -258,9 +247,7 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Volume = NULL; Vertex = NULL; Coord = NULL; Coord_Old = NULL; Coord_Sum = NULL; Coord_n = NULL; Coord_n1 = NULL; Coord_p1 = NULL; - GridVel = NULL; GridVel_Grad = NULL; GridVel_Old = NULL; - GridVel_n = NULL; GridVel_n1 = NULL; - + GridVel = NULL; GridVel_Grad = NULL; /*--- Volume (0 -> Vol_nP1, 1-> Vol_n, 2 -> Vol_nM1 ) and coordinates of the control volume ---*/ if ( config->GetUnsteady_Simulation() == NO ) { @@ -314,10 +301,6 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord if (config->GetGrid_Movement()) { GridVel = new su2double[nDim]; - GridVel_Old = new su2double[nDim]; - GridVel_n = new su2double[nDim]; - GridVel_n1 = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim ++) GridVel[iDim] = 0.0; @@ -335,7 +318,6 @@ CPoint::CPoint(su2double val_coord_0, su2double val_coord_1, su2double val_coord Coord_p1 = new su2double[nDim]; Coord_n = new su2double[nDim]; Coord_n1 = new su2double[nDim]; - Coord_Old = new su2double[nDim]; for (iDim = 0; iDim < nDim; iDim ++) { Coord_p1[iDim] = Coord[iDim]; Coord_n[iDim] = Coord[iDim]; @@ -363,9 +345,6 @@ CPoint::~CPoint() { if (Coord_n1 != NULL) delete[] Coord_n1; if (Coord_p1 != NULL) delete[] Coord_p1; if (GridVel != NULL) delete[] GridVel; - if (GridVel_Old != NULL) delete[] GridVel_Old; - if (GridVel_n != NULL) delete[] GridVel_n; - if (GridVel_n1 != NULL) delete[] GridVel_n1; if (GridVel_Grad != NULL) { for (unsigned short iDim = 0; iDim < nDim; iDim++) delete [] GridVel_Grad[iDim]; diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp index bb3105757918..60b251f9393c 100755 --- a/Common/src/geometry_structure.cpp +++ b/Common/src/geometry_structure.cpp @@ -16059,7 +16059,7 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config) { if (compressible) { skipVar += skipMult*(nDim+2); } if (sst && !frozen_visc) { skipVar += skipMult*2;} if (sa && !frozen_visc) { skipVar += skipMult*1;} - //CVC: Debug: if (grid_movement) { skipVar += nDim;} + if (grid_movement) { skipVar += nDim;} } else if (Kind_Solver == DISC_ADJ_HEAT) { skipVar += 1; diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp index a0d95f6b85c3..e3e807630597 100755 --- a/Common/src/grid_movement_structure.cpp +++ b/Common/src/grid_movement_structure.cpp @@ -1904,7 +1904,7 @@ void CVolumetricMovement::Rigid_Rotation(CGeometry *geometry, CConfig *config, su2double dtheta, dphi, dpsi, cosTheta, sinTheta; su2double cosPhi, sinPhi, cosPsi, sinPsi; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); + bool adjoint = config->GetContinuous_Adjoint(); /*--- Problem dimension and physical time step ---*/ @@ -2067,7 +2067,7 @@ void CVolumetricMovement::Rigid_Pitching(CGeometry *geometry, CConfig *config, u unsigned short nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); + bool adjoint = config->GetContinuous_Adjoint(); /*--- Retrieve values from the config file ---*/ @@ -2214,7 +2214,7 @@ void CVolumetricMovement::Rigid_Plunging(CGeometry *geometry, CConfig *config, u unsigned short iDim, nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); + bool adjoint = config->GetContinuous_Adjoint(); /*--- Retrieve values from the config file ---*/ @@ -2346,7 +2346,7 @@ void CVolumetricMovement::Rigid_Translation(CGeometry *geometry, CConfig *config unsigned short iDim, nDim = geometry->GetnDim(); unsigned long iPoint; bool harmonic_balance = (config->GetUnsteady_Simulation() == HARMONIC_BALANCE); - bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); // This should be consistent over all Rigid_* routines + bool adjoint = config->GetContinuous_Adjoint(); /*--- Retrieve values from the config file ---*/ @@ -6441,7 +6441,7 @@ void CSurfaceMovement::SetExternal_Deformation(CGeometry *geometry, CConfig *con string DV_Filename, UnstExt, text_line; ifstream surface_positions; bool unsteady = config->GetUnsteady_Simulation(); - bool adjoint = (config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint()); + bool adjoint = config->GetContinuous_Adjoint(); /*--- Load stuff from config ---*/ diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 9c2d138e7d60..bdae734f1f79 100755 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2127,7 +2127,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, unsigned short iMesh; int Direct_Iter; bool heat = config[val_iZone]->GetWeakly_Coupled_Heat(); - bool grid_movement_bool = config[val_iZone]->GetGrid_Movement(); /*--- Read the target pressure for inverse design. ---------------------------------------------*/ if (config[val_iZone]->GetInvDesign_Cp() == YES) @@ -2154,6 +2153,7 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if (dual_time_2nd) { /*--- Load solution at timestep n-2 ---*/ + LoadUnsteady_Solution(geometry, solver,config, val_iZone, val_iInst, Direct_Iter-2); /*--- Push solution back to correct array ---*/ @@ -2162,13 +2162,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(); solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n1(); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n1(); - - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n1(); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(); solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n1(); @@ -2183,6 +2176,7 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if (dual_time) { /*--- Load solution at timestep n-1 ---*/ + LoadUnsteady_Solution(geometry, solver,config, val_iZone, val_iInst, Direct_Iter-1); /*--- Push solution back to correct array ---*/ @@ -2190,10 +2184,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(); } @@ -2213,12 +2203,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, if ((ExtIter > 0) && dual_time){ - /*--- - Here the primal solutions (only working variables) are loaded and put in the correct order - into containers. For ALE the mesh coordinates and the Grid-Velocities have to be put into the - correct containers as well, i.e. follow the same logic for the solution. - ---*/ - /*--- Load solution timestep n-1 | n-2 for DualTimestepping 1st | 2nd order ---*/ if (dual_time_1st){ LoadUnsteady_Solution(geometry, solver,config, val_iInst, val_iZone, Direct_Iter - 1); @@ -2232,10 +2216,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_OldSolution(); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_Old(); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_Old(); - } if (turbulent){ solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_OldSolution(); } @@ -2250,10 +2230,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->SetSolution(solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_time_n()); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_n()); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_n()); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->SetSolution(solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_time_n()); } @@ -2267,10 +2243,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_Old()); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_Old()); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_Old()); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_Old()); } @@ -2285,10 +2257,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n(solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_time_n1()); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_n1()); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_n1()); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n(solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_time_n1()); } @@ -2301,10 +2269,6 @@ void CDiscAdjFluidIteration::Preprocess(COutput *output, for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { for(iPoint=0; iPointGetnPoint();iPoint++) { solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->Set_Solution_time_n1(solver[val_iZone][val_iInst][iMesh][FLOW_SOL]->node[iPoint]->GetSolution_Old()); - if (grid_movement_bool) { - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetCoord_n1(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetCoord_Old()); - geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->SetGridVel_n1(geometry[val_iZone][val_iInst][iMesh]->node[iPoint]->GetGridVel_Old()); - } if (turbulent) { solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->Set_Solution_time_n1(solver[val_iZone][val_iInst][iMesh][TURB_SOL]->node[iPoint]->GetSolution_Old()); } diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index 7deb754347e3..0bf1abc05f1b 100644 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -13031,7 +13031,7 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- Update the old geometry (coordinates n and n-1) in dual time-stepping strategy ---*/ - if (dual_time && grid_movement && (config->GetKind_GridMovement() != RIGID_MOTION)) + if (dual_time && grid_movement) Restart_OldGeometry(geometry[MESH_0], config); delete [] Coord; From 02592d4d5707f14cb38ae4185948638b9e78cd5b Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 22 Aug 2019 09:57:11 +0200 Subject: [PATCH 19/31] Reverted all excecutable file permissions for cpp/hpp/inl files. --- Common/include/dual_grid_structure.hpp | 0 Common/include/dual_grid_structure.inl | 0 Common/src/config_structure.cpp | 0 Common/src/dual_grid_structure.cpp | 0 Common/src/geometry_structure.cpp | 0 Common/src/grid_movement_structure.cpp | 0 SU2_CFD/include/solver_structure.hpp | 0 SU2_CFD/src/drivers/CDriver.cpp | 0 SU2_CFD/src/iteration_structure.cpp | 0 SU2_CFD/src/numerics_direct_mean_inc.cpp | 0 SU2_CFD/src/output_su2.cpp | 0 SU2_CFD/src/solver_structure.cpp | 0 12 files changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 Common/include/dual_grid_structure.hpp mode change 100755 => 100644 Common/include/dual_grid_structure.inl mode change 100755 => 100644 Common/src/config_structure.cpp mode change 100755 => 100644 Common/src/dual_grid_structure.cpp mode change 100755 => 100644 Common/src/geometry_structure.cpp mode change 100755 => 100644 Common/src/grid_movement_structure.cpp mode change 100755 => 100644 SU2_CFD/include/solver_structure.hpp mode change 100755 => 100644 SU2_CFD/src/drivers/CDriver.cpp mode change 100755 => 100644 SU2_CFD/src/iteration_structure.cpp mode change 100755 => 100644 SU2_CFD/src/numerics_direct_mean_inc.cpp mode change 100755 => 100644 SU2_CFD/src/output_su2.cpp mode change 100755 => 100644 SU2_CFD/src/solver_structure.cpp diff --git a/Common/include/dual_grid_structure.hpp b/Common/include/dual_grid_structure.hpp old mode 100755 new mode 100644 diff --git a/Common/include/dual_grid_structure.inl b/Common/include/dual_grid_structure.inl old mode 100755 new mode 100644 diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp old mode 100755 new mode 100644 diff --git a/Common/src/dual_grid_structure.cpp b/Common/src/dual_grid_structure.cpp old mode 100755 new mode 100644 diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp old mode 100755 new mode 100644 diff --git a/Common/src/grid_movement_structure.cpp b/Common/src/grid_movement_structure.cpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/include/solver_structure.hpp b/SU2_CFD/include/solver_structure.hpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/src/output_su2.cpp b/SU2_CFD/src/output_su2.cpp old mode 100755 new mode 100644 diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp old mode 100755 new mode 100644 From e5d456e5f9abf357eda031425cfb972dd822005f Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 22 Aug 2019 23:50:57 +0200 Subject: [PATCH 20/31] Little cleaning of CSourceIncRotatingFrame_Flow. --- SU2_CFD/include/numerics_structure.hpp | 5 ++ SU2_CFD/src/numerics_direct_mean_inc.cpp | 103 +++++++++++------------ SU2_CFD/src/solver_direct_mean_inc.cpp | 31 ++++--- 3 files changed, 68 insertions(+), 71 deletions(-) diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index 6b62fb0157ac..815d44f919d8 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5249,6 +5249,11 @@ class CSourceIncBodyForce : public CNumerics { * \ingroup SourceDiscr */ class CSourceIncRotatingFrame_Flow : public CNumerics { + +private: + su2double Omega[3]; + bool implicit; + public: /*! diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 9bd365a82747..8c4f312f618a 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -981,78 +981,71 @@ 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) { - - Gamma = config->GetGamma(); - Gamma_Minus_One = Gamma - 1.0; - -} - -CSourceIncRotatingFrame_Flow::~CSourceIncRotatingFrame_Flow(void) { } - -void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { +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); - unsigned short iDim, iVar, jVar; - su2double Omega[3] = {0,0,0}, Momentum[3] = {0,0,0}, Velocity_i[3] = {0,0,0}; + Gamma = config->GetGamma(); + Gamma_Minus_One = Gamma - 1.0; - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); + /*--- Retrieve the angular velocity vector from config. ---*/ + for (unsigned short iDim = 0; iDim < 3; iDim++) + Omega[iDim] = config->GetRotation_Rate(iDim)/config->GetOmega_Ref(); - /*--- Retrieve the angular velocity vector from config. ---*/ +} - for (iDim = 0; iDim < 3; iDim++){ - Omega[iDim] = config->GetRotation_Rate(iDim)/config->GetOmega_Ref(); - } +CSourceIncRotatingFrame_Flow::~CSourceIncRotatingFrame_Flow(void) { } - /*--- Primitive variables at point i and j ---*/ +void CSourceIncRotatingFrame_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, CConfig *config) { - DensityInc_i = V_i[nDim+2]; + unsigned short iDim, iVar, jVar; + su2double Momentum[3] = {0,0,0}, + Velocity_i[3] = {0,0,0}; - for (iDim = 0; iDim < nDim; iDim++) { - Velocity_i[iDim] = V_i[iDim+1]; - } + /*--- Primitive variables plus momentum at the node (point i) ---*/ - /*--- Get the momentum vector at the current node. ---*/ + DensityInc_i = V_i[nDim+2]; - for (iDim = 0; iDim < nDim; iDim++) { - Momentum[iDim] = DensityInc_i*Velocity_i[iDim]; + 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 as ( Omega X Rho-U ) ---*/ + /*--- 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; + 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; + 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; - } - } + /*--- 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) { diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 07d94540b941..800bdefd222a 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -2091,39 +2091,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 primitive variables ---*/ - + numerics->SetPrimitive(node[iPoint]->GetPrimitive(), NULL); /*--- Set incompressible density ---*/ - numerics->SetDensity(node[iPoint]->GetDensity(), - node[iPoint]->GetDensity()); - + 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 ---*/ From 2cabc267bf259fd3d7741631d143339f6ae4b894 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Fri, 23 Aug 2019 00:13:38 +0200 Subject: [PATCH 21/31] Remove trailing whitespaces from CSourceIncRotatingFrame_Flow class declaration. --- SU2_CFD/include/numerics_structure.hpp | 58 +++++++++++++------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp index 815d44f919d8..3709b4f03f71 100644 --- a/SU2_CFD/include/numerics_structure.hpp +++ b/SU2_CFD/include/numerics_structure.hpp @@ -5244,38 +5244,38 @@ class CSourceIncBodyForce : public CNumerics { }; /*! - * \class CSourceIncRotatingFrame_Flow - * \brief Class for a rotating frame source term. - * \ingroup SourceDiscr - */ -class CSourceIncRotatingFrame_Flow : public CNumerics { + * \class CSourceIncRotatingFrame_Flow + * \brief Class for a rotating frame source term. + * \ingroup SourceDiscr + */ +class CSourceIncRotatingFrame_Flow : public CNumerics { private: - su2double Omega[3]; - bool implicit; + 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); -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); + /*! + * \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); }; From 4129bacd144ac6769e4b7d8da5f4c48349b75a00 Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 11:49:50 +0100 Subject: [PATCH 22/31] Added new incompressible ale pitching airfoil case --- .../config_incomp_turb_sa.cfg | 126 ++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg 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..7a77312b3f7d --- /dev/null +++ b/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg @@ -0,0 +1,126 @@ +% ------------------------- PHYSICAL PROBLEM ----------------------------------% +% +PHYSICAL_PROBLEM= NAVIER_STOKES +REGIME_TYPE= INCOMPRESSIBLE +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 +UNST_INT_ITER= 201 + +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 +EXT_ITER= 10000 + +% ------------------------ 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 + From 6a9c87d08269cfee0b15b82aaa6984b916c5603f Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 14:10:25 +0100 Subject: [PATCH 23/31] Updated grid velocity correction comments --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 8c4f312f618a..7cf127f5eef8 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -230,7 +230,7 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } } - /*--- Jacobian contributions due to grid motion ---*/ + /*--- Corrections due to grid motion ---*/ if (grid_movement) { /*--- Recompute conservative variables ---*/ @@ -244,9 +244,12 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_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++){ @@ -385,7 +388,7 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } - /*--- Jacobian contributions due to grid motion ---*/ + /*--- Corrections due to grid motion ---*/ if (grid_movement) { /*--- Recompute conservative variables ---*/ @@ -399,9 +402,12 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ 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++){ @@ -594,7 +600,7 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } - /*--- Jacobian contributions due to grid motion ---*/ + /*--- Corrections due to grid motion ---*/ if (grid_movement) { /*--- Recompute conservative variables ---*/ @@ -607,9 +613,12 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ 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++){ From 3781d7cb72a54bb824f77fa81d927bedad6a16ee Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 15:40:06 +0100 Subject: [PATCH 24/31] Call SetPreconditioner for Jacobian matrix for static and dynamic problems --- SU2_CFD/src/solver_direct_mean_inc.cpp | 133 +++---------------------- 1 file changed, 14 insertions(+), 119 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 800bdefd222a..14a566aa3f85 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -5775,12 +5775,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Local variables ---*/ - unsigned short iVar, jVar, iMarker, iDim, jDim; + unsigned short iVar, jVar, iMarker, iDim; unsigned long iPoint, jPoint, iEdge, iVertex; su2double Density, Cp; - su2double BetaInc2, dRhodT, Temperature, oneOverCp; - su2double Velocity[3] = {0.0,0.0,0.0}; su2double *V_time_nM1, *V_time_n, *V_time_nP1; su2double U_time_nM1[5], U_time_n[5], U_time_nP1[5]; su2double Volume_nM1, Volume_nP1, TimeStep; @@ -5822,27 +5820,11 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = node[iPoint]->GetSolution_time_n(); V_time_nP1 = node[iPoint]->GetSolution(); - /*--- Access the primitive variables at this node. ---*/ + /*--- Access the density and Cp at this node (constant for now). ---*/ Density = node[iPoint]->GetDensity(); - BetaInc2 = node[iPoint]->GetBetaInc2(); Cp = node[iPoint]->GetSpecificHeatCp(); - oneOverCp = 1.0/Cp; - 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; - } - + /*--- Compute the conservative variable vector for all time levels. ---*/ U_time_nM1[0] = Density; @@ -5885,61 +5867,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) @@ -6001,26 +5935,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = node[iPoint]->GetSolution_time_n(); - /*--- Access the primitive variables at this node. ---*/ + /*--- Access the density and Cp at this node (constant for now). ---*/ Density = node[iPoint]->GetDensity(); - BetaInc2 = node[iPoint]->GetBetaInc2(); Cp = node[iPoint]->GetSpecificHeatCp(); - oneOverCp = 1.0/Cp; - 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; - } /*--- Compute the conservative variable vector for all time levels. ---*/ @@ -6174,36 +6092,13 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver if (!energy) Residual[nDim+1] = 0.0; LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - - /*--- 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; + SetPreconditioner(config, iPoint); + for (iVar = 0; iVar < nVar; iVar++) { + for (jVar = 0; jVar < nVar; jVar++) { + Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar]; } - Jacobian_i[nDim+1][jDim+1] = 0.0; } - 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) From e3951df8a2d050ff5b4a7ce8c02ef5fe9a3216e2 Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 17:21:09 +0100 Subject: [PATCH 25/31] Updated to new dynamic grid boolean and added preaccumulation --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 37 +++++++++++++++++------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index de347fe3489b..6aad6c23816c 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -91,7 +91,7 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+9); AD::SetPreaccIn(V_j, nDim+9); AD::SetPreaccIn(Normal, nDim); - if (grid_movement) { + if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); } @@ -129,7 +129,7 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J /*--- Projected velocity adjustment due to mesh motion ---*/ - if (grid_movement) { + if (dynamic_grid) { ProjGridVel = 0.0; for (iDim = 0; iDim < nDim; iDim++) { ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; @@ -232,7 +232,7 @@ void CUpwFDSInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_J } /*--- Corrections due to grid motion ---*/ - if (grid_movement) { + if (dynamic_grid) { /*--- Recompute conservative variables ---*/ @@ -326,9 +326,15 @@ CCentJSTInc_Flow::~CCentJSTInc_Flow(void) { void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_Jacobian_i, su2double **val_Jacobian_j, CConfig *config) { - //Preaccumulation? 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); + } /*--- Primitive variables at point i and j ---*/ @@ -391,7 +397,7 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } /*--- Corrections due to grid motion ---*/ - if (grid_movement) { + if (dynamic_grid) { /*--- Recompute conservative variables ---*/ @@ -438,7 +444,7 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ /*--- Projected velocity adjustment due to mesh motion ---*/ - if (grid_movement) { + if (dynamic_grid) { ProjGridVel = 0.0; for (iDim = 0; iDim < nDim; iDim++) { ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; @@ -492,7 +498,9 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - + + AD::SetPreaccOut(val_residual, nVar); + AD::EndPreacc(); } CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { @@ -541,6 +549,13 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ 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; + 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); + } + /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -604,7 +619,7 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } /*--- Corrections due to grid motion ---*/ - if (grid_movement) { + if (dynamic_grid) { /*--- Recompute conservative variables ---*/ @@ -651,7 +666,7 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ /*--- Projected velocity adjustment due to mesh motion ---*/ - if (grid_movement) { + if (dynamic_grid) { ProjGridVel = 0.0; for (iDim = 0; iDim < nDim; iDim++) { ProjGridVel += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; @@ -699,7 +714,9 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - + + AD::SetPreaccOut(val_residual, nVar); + AD::EndPreacc(); } CAvgGradInc_Flow::CAvgGradInc_Flow(unsigned short val_nDim, From 4051b40e0d266ede0b941383a1ba4a3aabc02b47 Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 17:21:58 +0100 Subject: [PATCH 26/31] Updated SOLVER definition in test case --- .../pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 index 7a77312b3f7d..e050a61645e2 100644 --- a/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg +++ b/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg @@ -1,7 +1,6 @@ % ------------------------- PHYSICAL PROBLEM ----------------------------------% % -PHYSICAL_PROBLEM= NAVIER_STOKES -REGIME_TYPE= INCOMPRESSIBLE +SOLVER= INC_NAVIER_STOKES KIND_TURB_MODEL= SA MATH_PROBLEM= DIRECT RESTART_SOL= NO From 16f6ee9f28483058eeaa49ae299b60e110fdd2e2 Mon Sep 17 00:00:00 2001 From: cvencro Date: Fri, 30 Aug 2019 23:28:08 +0100 Subject: [PATCH 27/31] revert preaccumulation to pass JST regression tests --- SU2_CFD/src/numerics_direct_mean_inc.cpp | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/SU2_CFD/src/numerics_direct_mean_inc.cpp b/SU2_CFD/src/numerics_direct_mean_inc.cpp index 6aad6c23816c..95df77910bcc 100644 --- a/SU2_CFD/src/numerics_direct_mean_inc.cpp +++ b/SU2_CFD/src/numerics_direct_mean_inc.cpp @@ -329,13 +329,6 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ 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); - } - /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -498,9 +491,6 @@ void CCentJSTInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - - AD::SetPreaccOut(val_residual, nVar); - AD::EndPreacc(); } CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_nVar, CConfig *config) : CNumerics(val_nDim, val_nVar, config) { @@ -549,13 +539,6 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ 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; - 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); - } - /*--- Primitive variables at point i and j ---*/ Pressure_i = V_i[0]; Pressure_j = V_j[0]; @@ -714,9 +697,6 @@ void CCentLaxInc_Flow::ComputeResidual(su2double *val_residual, su2double **val_ } } } - - AD::SetPreaccOut(val_residual, nVar); - AD::EndPreacc(); } CAvgGradInc_Flow::CAvgGradInc_Flow(unsigned short val_nDim, From 92b58d73d5e8ce9fdd11fe682e5e915ba0764277 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Tue, 3 Sep 2019 09:49:36 +0200 Subject: [PATCH 28/31] PR767: Added serial and parallel regression test. --- .travis.yml | 2 +- TestCases/parallel_regression.py | 12 ++++++++++++ TestCases/serial_regression.py | 12 ++++++++++++ .../config_incomp_turb_sa.cfg | 15 +++++++++------ 4 files changed, 34 insertions(+), 7 deletions(-) diff --git a/.travis.yml b/.travis.yml index 07b1bda0b8d3..c72874d4d439 100644 --- a/.travis.yml +++ b/.travis.yml @@ -72,7 +72,7 @@ install: before_script: # Get the test cases - - git clone -b develop https://github.com/su2code/TestCases.git ./TestData + - git clone -b feature_incompressible_ale https://github.com/su2code/TestCases.git ./TestData - cp -R ./TestData/* ./TestCases/ # Get the tutorial cases 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 index e050a61645e2..7905138ebaef 100644 --- a/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg +++ b/TestCases/unsteady/pitching_naca0015_rans_inc/config_incomp_turb_sa.cfg @@ -10,12 +10,16 @@ RESTART_SOL= NO 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 - -SINGLEZONE_DRIVER= YES -TIME_DOMAIN= YES -TIME_ITER= 151 -ITER= 201 +% +% New driver +%SINGLEZONE_DRIVER= YES +%TIME_DOMAIN= YES +%TIME_ITER= 151 +%ITER= 201 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % GRID_MOVEMENT= RIGID_MOTION @@ -52,7 +56,6 @@ MARKER_MONITORING= ( airfoil ) NUM_METHOD_GRAD= GREEN_GAUSS CFL_NUMBER= 10.0 CFL_ADAPT= NO -EXT_ITER= 10000 % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % From fa587e2936ebc207472db34af4adfb8c81dc398a Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Tue, 3 Sep 2019 13:01:13 +0200 Subject: [PATCH 29/31] Revert .travis.yml to develop. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index c72874d4d439..07b1bda0b8d3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -72,7 +72,7 @@ install: before_script: # Get the test cases - - git clone -b feature_incompressible_ale https://github.com/su2code/TestCases.git ./TestData + - git clone -b develop https://github.com/su2code/TestCases.git ./TestData - cp -R ./TestData/* ./TestCases/ # Get the tutorial cases From a91c843d777417136d10166d075e1bf073acef18 Mon Sep 17 00:00:00 2001 From: cvencro Date: Tue, 3 Sep 2019 17:12:34 +0100 Subject: [PATCH 30/31] Error catch for rotating frame in incompressible solver --- Common/src/config_structure.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp index 7f58410c0feb..154b872abe20 100644 --- a/Common/src/config_structure.cpp +++ b/Common/src/config_structure.cpp @@ -4297,6 +4297,12 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } } + /*--- Rotating frame is not yet supported with the incompressible solver. ---*/ + + 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 pressure drop objective function is selected. ---*/ From a7d10a1c3b89e3bddd0523c3122d089c0600d78f Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Wed, 4 Sep 2019 13:54:56 +0200 Subject: [PATCH 31/31] Remove Unused Variable which throws compiler warning. --- SU2_CFD/src/solver_direct_mean_inc.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 513c0f7ccacf..d45d0aaafdce 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -5784,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 ---*/