diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp index e8a747f09ffb..701c517f8b9f 100644 --- a/SU2_CFD/src/solver_direct_mean.cpp +++ b/SU2_CFD/src/solver_direct_mean.cpp @@ -3293,7 +3293,7 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container bool disc_adjoint = config->GetDiscrete_Adjoint(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); - bool limiter = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); bool center_jst = center && (config->GetKind_Centered_Flow() == JST); bool engine = ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)); @@ -3682,8 +3682,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain unsigned long ExtIter = config->GetExtIter(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); - bool disc_adjoint = config->GetDiscrete_Adjoint(); - bool limiter = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool grid_movement = config->GetGrid_Movement(); bool roe_turkel = (config->GetKind_Upwind_Flow() == TURKEL); bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS ); @@ -5819,9 +5818,20 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { unsigned long iEdge, iPoint, jPoint; unsigned short iVar, iDim; su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Primitive, *Primitive_i, *Primitive_j, *LocalMinPrimitive, *LocalMaxPrimitive, - *GlobalMinPrimitive, *GlobalMaxPrimitive, + *Primitive, *Primitive_i, *Primitive_j, + *LocalMinPrimitive = NULL, *LocalMaxPrimitive = NULL, + *GlobalMinPrimitive = NULL, *GlobalMaxPrimitive = NULL, dave, LimK, eps2, eps1, dm, dp, du, y, limiter; + +#ifdef CODI_REVERSE_TYPE + bool TapeActive = false; + + if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { + /*--- If limiters are frozen do not record the computation ---*/ + TapeActive = AD::globalTape.isActive(); + AD::StopRecording(); + } +#endif dave = config->GetRefElemLength(); LimK = config->GetVenkat_LimiterCoeff(); @@ -5965,33 +5975,34 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { if ((config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN) || (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - /*--- Allocate memory for the max and min primitive value --*/ - - LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; - LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; - - /*--- Compute the max value and min value of the solution ---*/ - - Primitive = node[0]->GetPrimitive(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = Primitive[iVar]; - LocalMaxPrimitive[iVar] = Primitive[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + + /*--- Allocate memory for the max and min primitive value --*/ - /*--- Get the primitive variables ---*/ + LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; + LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; - Primitive = node[iPoint]->GetPrimitive(); + /*--- Compute the max value and min value of the solution ---*/ + Primitive = node[0]->GetPrimitive(); for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); - LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); + LocalMinPrimitive[iVar] = Primitive[iVar]; + LocalMaxPrimitive[iVar] = Primitive[iVar]; + } + + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + + /*--- Get the primitive variables ---*/ + + Primitive = node[iPoint]->GetPrimitive(); + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); + LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); + } + } - } - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { #ifdef HAVE_MPI SU2_MPI::Allreduce(LocalMinPrimitive, GlobalMinPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); SU2_MPI::Allreduce(LocalMaxPrimitive, GlobalMaxPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); @@ -6011,16 +6022,22 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { Gradient_j = node[jPoint]->GetGradient_Primitive(); Coord_i = geometry->node[iPoint]->GetCoord(); Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nPrimVarGrad, nDim); - AD::SetPreaccIn(Gradient_j, nPrimVarGrad, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - AD::SetPreaccIn(eps2); for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + + AD::StartPreacc(); + AD::SetPreaccIn(Gradient_i[iVar], nDim); + AD::SetPreaccIn(Gradient_j[iVar], nDim); + AD::SetPreaccIn(Coord_i, nDim); + AD::SetPreaccIn(Coord_j, nDim); + AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + AD::SetPreaccIn(GlobalMaxPrimitive[iVar]); + AD::SetPreaccIn(GlobalMinPrimitive[iVar]); eps1 = LimK * (GlobalMaxPrimitive[iVar] - GlobalMinPrimitive[iVar]); eps2 = eps1*eps1; } @@ -6029,11 +6046,6 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { eps2 = eps1*eps1*eps1; } - AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); - /*--- Calculate the interface left gradient, delta- (dm) ---*/ dm = 0.0; @@ -6068,14 +6080,14 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { AD::SetPreaccOut(node[jPoint]->GetLimiter_Primitive()[iVar]); } + AD::EndPreacc(); } - - AD::EndPreacc(); - } - delete [] LocalMinPrimitive; delete [] GlobalMinPrimitive; - delete [] LocalMaxPrimitive; delete [] GlobalMaxPrimitive; + if (LocalMinPrimitive != NULL) delete [] LocalMinPrimitive; + if (LocalMaxPrimitive != NULL) delete [] LocalMaxPrimitive; + if (GlobalMinPrimitive != NULL) delete [] GlobalMinPrimitive; + if (GlobalMaxPrimitive != NULL) delete [] GlobalMaxPrimitive; } @@ -6091,6 +6103,9 @@ void CEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) { InitiateComms(geometry, config, PRIMITIVE_LIMITER); CompleteComms(geometry, config, PRIMITIVE_LIMITER); +#ifdef CODI_REVERSE_TYPE + if (TapeActive) AD::StartRecording(); +#endif } void CEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { @@ -15374,8 +15389,8 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); bool center_jst = center && config->GetKind_Centered_Flow() == JST; - bool limiter_flow = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); - bool limiter_turb = ((config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter_flow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); + bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool limiter_adjflow = (cont_adjoint && (config->GetKind_SlopeLimit_AdjFlow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter())); bool fixed_cl = config->GetFixed_CL_Mode(); bool engine = ((config->GetnMarker_EngineInflow() != 0) || (config->GetnMarker_EngineExhaust() != 0)); diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp index 3e64d30b3cd2..9494a3d699ec 100644 --- a/SU2_CFD/src/solver_direct_mean_inc.cpp +++ b/SU2_CFD/src/solver_direct_mean_inc.cpp @@ -1528,7 +1528,7 @@ void CIncEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contai bool disc_adjoint = config->GetDiscrete_Adjoint(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); - bool limiter = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool center = ((config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED)); bool center_jst = center && (config->GetKind_Centered_Flow() == JST); bool fixed_cl = config->GetFixed_CL_Mode(); @@ -1902,8 +1902,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont unsigned long ExtIter = config->GetExtIter(); bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); - bool disc_adjoint = config->GetDiscrete_Adjoint(); - bool limiter = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool grid_movement = config->GetGrid_Movement(); bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; @@ -3778,10 +3777,21 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) unsigned long iEdge, iPoint, jPoint; unsigned short iVar, iDim; su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Primitive, *Primitive_i, *Primitive_j, *LocalMinPrimitive, *LocalMaxPrimitive, - *GlobalMinPrimitive, *GlobalMaxPrimitive, + *Primitive, *Primitive_i, *Primitive_j, + *LocalMinPrimitive = NULL, *LocalMaxPrimitive = NULL, + *GlobalMinPrimitive = NULL, *GlobalMaxPrimitive = NULL, dave, LimK, eps2, eps1, dm, dp, du, y, limiter; +#ifdef CODI_REVERSE_TYPE + bool TapeActive = false; + + if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { + /*--- If limiters are frozen do not record the computation ---*/ + TapeActive = AD::globalTape.isActive(); + AD::StopRecording(); + } +#endif + dave = config->GetRefElemLength(); LimK = config->GetVenkat_LimiterCoeff(); @@ -3925,33 +3935,34 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) if ((config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN) || (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - /*--- Allocate memory for the max and min primitive value --*/ - - LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; - LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; - - /*--- Compute the max value and min value of the solution ---*/ - - Primitive = node[0]->GetPrimitive(); - for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = Primitive[iVar]; - LocalMaxPrimitive[iVar] = Primitive[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { - /*--- Get the primitive variables ---*/ + /*--- Allocate memory for the max and min primitive value --*/ - Primitive = node[iPoint]->GetPrimitive(); - + LocalMinPrimitive = new su2double [nPrimVarGrad]; GlobalMinPrimitive = new su2double [nPrimVarGrad]; + LocalMaxPrimitive = new su2double [nPrimVarGrad]; GlobalMaxPrimitive = new su2double [nPrimVarGrad]; + + /*--- Compute the max value and min value of the solution ---*/ + + Primitive = node[0]->GetPrimitive(); for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); - LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); + LocalMinPrimitive[iVar] = Primitive[iVar]; + LocalMaxPrimitive[iVar] = Primitive[iVar]; } - } + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + + /*--- Get the primitive variables ---*/ + + Primitive = node[iPoint]->GetPrimitive(); + + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + LocalMinPrimitive[iVar] = min (LocalMinPrimitive[iVar], Primitive[iVar]); + LocalMaxPrimitive[iVar] = max (LocalMaxPrimitive[iVar], Primitive[iVar]); + } + + } - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { #ifdef HAVE_MPI SU2_MPI::Allreduce(LocalMinPrimitive, GlobalMinPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); SU2_MPI::Allreduce(LocalMaxPrimitive, GlobalMaxPrimitive, nPrimVarGrad, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); @@ -3971,15 +3982,22 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) Gradient_j = node[jPoint]->GetGradient_Primitive(); Coord_i = geometry->node[iPoint]->GetCoord(); Coord_j = geometry->node[jPoint]->GetCoord(); - - AD::StartPreacc(); - AD::SetPreaccIn(Gradient_i, nPrimVarGrad, nDim); - AD::SetPreaccIn(Gradient_j, nPrimVarGrad, nDim); - AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); - + for (iVar = 0; iVar < nPrimVarGrad; iVar++) { + + AD::StartPreacc(); + AD::SetPreaccIn(Gradient_i[iVar], nDim); + AD::SetPreaccIn(Gradient_j[iVar], nDim); + AD::SetPreaccIn(Coord_i, nDim); + AD::SetPreaccIn(Coord_j, nDim); + AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + AD::SetPreaccIn(GlobalMaxPrimitive[iVar]); + AD::SetPreaccIn(GlobalMinPrimitive[iVar]); eps1 = LimK * (GlobalMaxPrimitive[iVar] - GlobalMinPrimitive[iVar]); eps2 = eps1*eps1; } @@ -3988,11 +4006,6 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) eps2 = eps1*eps1*eps1; } - AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); - /*--- Calculate the interface left gradient, delta- (dm) ---*/ dm = 0.0; @@ -4027,14 +4040,14 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) AD::SetPreaccOut(node[jPoint]->GetLimiter_Primitive()[iVar]); } + AD::EndPreacc(); } - - AD::EndPreacc(); - } - delete [] LocalMinPrimitive; delete [] GlobalMinPrimitive; - delete [] LocalMaxPrimitive; delete [] GlobalMaxPrimitive; + if (LocalMinPrimitive != NULL) delete [] LocalMinPrimitive; + if (LocalMaxPrimitive != NULL) delete [] LocalMaxPrimitive; + if (GlobalMinPrimitive != NULL) delete [] GlobalMinPrimitive; + if (GlobalMaxPrimitive != NULL) delete [] GlobalMaxPrimitive; } @@ -4049,7 +4062,10 @@ void CIncEulerSolver::SetPrimitive_Limiter(CGeometry *geometry, CConfig *config) InitiateComms(geometry, config, PRIMITIVE_LIMITER); CompleteComms(geometry, config, PRIMITIVE_LIMITER); - + +#ifdef CODI_REVERSE_TYPE + if (TapeActive) AD::StartRecording(); +#endif } void CIncEulerSolver::SetFarfield_AoA(CGeometry *geometry, CSolver **solver_container, @@ -7437,8 +7453,8 @@ void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); bool center = ((config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED)); bool center_jst = center && config->GetKind_Centered_Flow() == JST; - bool limiter_flow = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); - bool limiter_turb = ((config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter_flow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); + bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()); bool limiter_adjflow = (cont_adjoint && (config->GetKind_SlopeLimit_AdjFlow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter())); bool fixed_cl = config->GetFixed_CL_Mode(); bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; diff --git a/SU2_CFD/src/solver_direct_turbulent.cpp b/SU2_CFD/src/solver_direct_turbulent.cpp index 6915b54de9f3..68891f3937ce 100644 --- a/SU2_CFD/src/solver_direct_turbulent.cpp +++ b/SU2_CFD/src/solver_direct_turbulent.cpp @@ -1148,10 +1148,7 @@ CTurbSASolver::~CTurbSASolver(void) { void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { unsigned long iPoint; - unsigned long ExtIter = config->GetExtIter(); - bool disc_adjoint = config->GetDiscrete_Adjoint(); - bool limiter_flow = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); - bool limiter_turb = ((config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (config->GetExtIter() <= config->GetLimiterIter()); unsigned short kind_hybridRANSLES = config->GetKind_HybridRANSLES(); su2double** PrimGrad_Flow = NULL; su2double* Vorticity = NULL; @@ -1176,8 +1173,6 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe if (limiter_turb) SetSolution_Limiter(geometry, config); - if (limiter_flow) solver_container[FLOW_SOL]->SetPrimitive_Limiter(geometry, config); - if (kind_hybridRANSLES != NO_HYBRIDRANSLES){ /*--- Set the vortex tilting coefficient at every node if required ---*/ @@ -3490,10 +3485,7 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain unsigned long iPoint; - unsigned long ExtIter = config->GetExtIter(); - bool disc_adjoint = config->GetDiscrete_Adjoint(); - bool limiter_flow = ((config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); - bool limiter_turb = ((config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (ExtIter <= config->GetLimiterIter()) && !(disc_adjoint && config->GetFrozen_Limiter_Disc())); + bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (config->GetExtIter() <= config->GetLimiterIter()); for (iPoint = 0; iPoint < nPoint; iPoint ++) { @@ -3514,8 +3506,6 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) SetSolution_Gradient_LS(geometry, config); if (limiter_turb) SetSolution_Limiter(geometry, config); - - if (limiter_flow) solver_container[FLOW_SOL]->SetPrimitive_Limiter(geometry, config); } diff --git a/SU2_CFD/src/solver_structure.cpp b/SU2_CFD/src/solver_structure.cpp index f16d6c772bd7..2a17d1524e4b 100644 --- a/SU2_CFD/src/solver_structure.cpp +++ b/SU2_CFD/src/solver_structure.cpp @@ -3140,10 +3140,21 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { unsigned long iEdge, iPoint, jPoint; unsigned short iVar, iDim; su2double **Gradient_i, **Gradient_j, *Coord_i, *Coord_j, - *Solution, *Solution_i, *Solution_j, *LocalMinSolution, *LocalMaxSolution, - *GlobalMinSolution, *GlobalMaxSolution, + *Solution, *Solution_i, *Solution_j, + *LocalMinSolution = NULL, *LocalMaxSolution = NULL, + *GlobalMinSolution = NULL, *GlobalMaxSolution = NULL, dave, LimK, eps1, eps2, dm, dp, du, ds, y, limiter, SharpEdge_Distance; +#ifdef CODI_REVERSE_TYPE + bool TapeActive = false; + + if (config->GetDiscrete_Adjoint() && config->GetFrozen_Limiter_Disc()) { + /*--- If limiters are frozen do not record the computation ---*/ + TapeActive = AD::globalTape.isActive(); + AD::StopRecording(); + } +#endif + dave = config->GetRefElemLength(); LimK = config->GetVenkat_LimiterCoeff(); @@ -3286,33 +3297,34 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { if ((config->GetKind_SlopeLimit() == VENKATAKRISHNAN) || (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG)) { - /*--- Allocate memory for the max and min solution value --*/ - - LocalMinSolution = new su2double [nVar]; GlobalMinSolution = new su2double [nVar]; - LocalMaxSolution = new su2double [nVar]; GlobalMaxSolution = new su2double [nVar]; - - /*--- Compute the max value and min value of the solution ---*/ - - Solution = node[iPoint]->GetSolution(); - for (iVar = 0; iVar < nVar; iVar++) { - LocalMinSolution[iVar] = Solution[iVar]; - LocalMaxSolution[iVar] = Solution[iVar]; - } - - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + + /*--- Allocate memory for the max and min solution value --*/ - /*--- Get the solution variables ---*/ + LocalMinSolution = new su2double [nVar]; GlobalMinSolution = new su2double [nVar]; + LocalMaxSolution = new su2double [nVar]; GlobalMaxSolution = new su2double [nVar]; - Solution = node[iPoint]->GetSolution(); + /*--- Compute the max value and min value of the solution ---*/ + Solution = node[iPoint]->GetSolution(); for (iVar = 0; iVar < nVar; iVar++) { - LocalMinSolution[iVar] = min (LocalMinSolution[iVar], Solution[iVar]); - LocalMaxSolution[iVar] = max (LocalMaxSolution[iVar], Solution[iVar]); + LocalMinSolution[iVar] = Solution[iVar]; + LocalMaxSolution[iVar] = Solution[iVar]; + } + + for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + + /*--- Get the solution variables ---*/ + + Solution = node[iPoint]->GetSolution(); + + for (iVar = 0; iVar < nVar; iVar++) { + LocalMinSolution[iVar] = min (LocalMinSolution[iVar], Solution[iVar]); + LocalMaxSolution[iVar] = max (LocalMaxSolution[iVar], Solution[iVar]); + } + } - } - - if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { #ifdef HAVE_MPI SU2_MPI::Allreduce(LocalMinSolution, GlobalMinSolution, nVar, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); SU2_MPI::Allreduce(LocalMaxSolution, GlobalMaxSolution, nVar, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); @@ -3339,8 +3351,20 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { AD::SetPreaccIn(Coord_i, nDim); AD::SetPreaccIn(Coord_j, nDim); for (iVar = 0; iVar < nVar; iVar++) { + + AD::StartPreacc(); + AD::SetPreaccIn(Gradient_i[iVar], nDim); + AD::SetPreaccIn(Gradient_j[iVar], nDim); + AD::SetPreaccIn(Coord_i, nDim); + AD::SetPreaccIn(Coord_j, nDim); + AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); + AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); if (config->GetKind_SlopeLimit_Flow() == VENKATAKRISHNAN_WANG) { + AD::SetPreaccIn(GlobalMaxSolution[iVar]); + AD::SetPreaccIn(GlobalMinSolution[iVar]); eps1 = LimK * (GlobalMaxSolution[iVar] - GlobalMinSolution[iVar]); eps2 = eps1*eps1; } @@ -3348,11 +3372,6 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { eps1 = LimK*dave; eps2 = eps1*eps1*eps1; } - - AD::SetPreaccIn(node[iPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[iPoint]->GetSolution_Min(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Max(iVar)); - AD::SetPreaccIn(node[jPoint]->GetSolution_Min(iVar)); /*--- Calculate the interface left gradient, delta- (dm) ---*/ @@ -3387,14 +3406,15 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { node[jPoint]->SetLimiter(iVar, limiter); AD::SetPreaccOut(node[jPoint]->GetLimiter()[iVar]); } + + AD::EndPreacc(); } - - AD::EndPreacc(); - } - delete [] LocalMinSolution; delete [] GlobalMinSolution; - delete [] LocalMaxSolution; delete [] GlobalMaxSolution; + if (LocalMinSolution != NULL) delete [] LocalMinSolution; + if (LocalMaxSolution != NULL) delete [] LocalMaxSolution; + if (GlobalMinSolution != NULL) delete [] GlobalMinSolution; + if (GlobalMaxSolution != NULL) delete [] GlobalMaxSolution; } @@ -3554,6 +3574,9 @@ void CSolver::SetSolution_Limiter(CGeometry *geometry, CConfig *config) { InitiateComms(geometry, config, SOLUTION_LIMITER); CompleteComms(geometry, config, SOLUTION_LIMITER); +#ifdef CODI_REVERSE_TYPE + if (TapeActive) AD::StartRecording(); +#endif } void CSolver::Gauss_Elimination(su2double** A, su2double* rhs, unsigned short nVar) {