Skip to content

Commit

Permalink
Merge pull request #950 from su2code/multiphysics_discadj_tweaks
Browse files Browse the repository at this point in the history
Small discrete adjoint tweaks and other fixes
  • Loading branch information
talbring authored Apr 29, 2020
2 parents 3380878 + f9efa1c commit 53fb61f
Show file tree
Hide file tree
Showing 24 changed files with 124 additions and 93 deletions.
7 changes: 3 additions & 4 deletions Common/include/CConfig.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -570,7 +570,7 @@ class CConfig {
unsigned short Linear_Solver_ILU_n; /*!< \brief ILU fill=in level. */
su2double SemiSpan; /*!< \brief Wing Semi span. */
su2double Roe_Kappa; /*!< \brief Relaxation of the Roe scheme. */
su2double Relaxation_Factor_AdjFlow; /*!< \brief Relaxation coefficient of the linear solver adjoint mean flow. */
su2double Relaxation_Factor_Adjoint; /*!< \brief Relaxation coefficient for variable updates of adjoint solvers. */
su2double Relaxation_Factor_CHT; /*!< \brief Relaxation coefficient for the update of conjugate heat variables. */
su2double AdjTurb_Linear_Error; /*!< \brief Min error of the turbulent adjoint linear solver for the implicit formulation. */
su2double EntropyFix_Coeff; /*!< \brief Entropy fix coefficient. */
Expand Down Expand Up @@ -3998,10 +3998,9 @@ class CConfig {
su2double GetLinear_Solver_Smoother_Relaxation(void) const { return Linear_Solver_Smoother_Relaxation; }

/*!
* \brief Get the relaxation coefficient of the linear solver for the implicit formulation.
* \return relaxation coefficient of the linear solver for the implicit formulation.
* \brief Get the relaxation factor for solution updates of adjoint solvers.
*/
su2double GetRelaxation_Factor_AdjFlow(void) const { return Relaxation_Factor_AdjFlow; }
su2double GetRelaxation_Factor_Adjoint(void) const { return Relaxation_Factor_Adjoint; }

/*!
* \brief Get the relaxation coefficient of the CHT coupling.
Expand Down
21 changes: 5 additions & 16 deletions Common/src/CConfig.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1670,8 +1670,8 @@ void CConfig::SetConfig_Options() {
addDoubleOption("LINEAR_SOLVER_SMOOTHER_RELAXATION", Linear_Solver_Smoother_Relaxation, 1.0);
/* DESCRIPTION: Custom number of threads used for additive domain decomposition for ILU and LU_SGS (0 is "auto"). */
addUnsignedLongOption("LINEAR_SOLVER_PREC_THREADS", Linear_Solver_Prec_Threads, 0);
/* DESCRIPTION: Relaxation of the flow equations solver for the implicit formulation */
addDoubleOption("RELAXATION_FACTOR_ADJFLOW", Relaxation_Factor_AdjFlow, 1.0);
/* DESCRIPTION: Relaxation factor for updates of adjoint variables. */
addDoubleOption("RELAXATION_FACTOR_ADJOINT", Relaxation_Factor_Adjoint, 1.0);
/* DESCRIPTION: Relaxation of the CHT coupling */
addDoubleOption("RELAXATION_FACTOR_CHT", Relaxation_Factor_CHT, 1.0);
/* DESCRIPTION: Roe coefficient */
Expand Down Expand Up @@ -2832,20 +2832,9 @@ void CConfig::SetConfig_Parsing(char case_filename[MAX_STRING_SIZE]) {
newString.append(": invalid option name");
newString.append(". Check current SU2 options in config_template.cfg.");
newString.append("\n");
if (!option_name.compare("EXT_ITER")) newString.append("Option EXT_ITER is deprecated as of v7.0. Please use TIME_ITER, OUTER_ITER or ITER \n"
"to specify the number of time iterations, outer multizone iterations or iterations, respectively.");
if (!option_name.compare("UNST_TIMESTEP")) newString.append("UNST_TIMESTEP is now TIME_STEP.\n");
if (!option_name.compare("UNST_TIME")) newString.append("UNST_TIME is now MAX_TIME.\n");
if (!option_name.compare("UNST_INT_ITER")) newString.append("UNST_INT_ITER is now INNER_ITER.\n");
if (!option_name.compare("RESIDUAL_MINVAL")) newString.append("RESIDUAL_MINVAL is now CONV_RESIDUAL_MINVAL.\n");
if (!option_name.compare("STARTCONV_ITER")) newString.append("STARTCONV_ITER is now CONV_STARTITER.\n");
if (!option_name.compare("CAUCHY_ELEMS")) newString.append("CAUCHY_ELEMS is now CONV_CAUCHY_ELEMS.\n");
if (!option_name.compare("CAUCHY_EPS")) newString.append("CAUCHY_EPS is now CONV_CAUCHY_EPS.\n");
if (!option_name.compare("OUTPUT_FORMAT")) newString.append("OUTPUT_FORMAT is now TABULAR_FORMAT.\n");
if (!option_name.compare("PHYSICAL_PROBLEM")) newString.append("PHYSICAL_PROBLEM is now SOLVER.\n");
if (!option_name.compare("REGIME_TYPE")) newString.append("REGIME_TYPE has been removed.\n "
"If you want use the incompressible solver, \n"
"use INC_EULER, INC_NAVIER_STOKES or INC_RANS as value of the SOLVER option.");
if (!option_name.compare("RELAXATION_FACTOR_ADJFLOW"))
newString.append("Option RELAXATION_FACTOR_ADJFLOW is now RELAXATION_FACTOR_ADJOINT, "
"and it also applies to discrete adjoint problems\n.");
errorString.append(newString);
err_count++;
continue;
Expand Down
6 changes: 6 additions & 0 deletions SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,12 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
*/
void Add_External_To_Solution(unsigned short iZone);

/*!
* \brief Puts Solution into SolutionOld.
* \param[in] iZone - Zone index.
*/
void Set_SolutionOld_To_Solution(unsigned short iZone);

/*!
* \brief Extract contribution of iZone to jZone with BGS relaxation.
* \param[in] iZone - Source zone (the one that was initialized).
Expand Down
4 changes: 2 additions & 2 deletions SU2_CFD/include/solvers/CFEASolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@ class CFEASolver : public CSolver {
su2double Conv_Check[3]; /*!< \brief Current values for convergence check: UTOL, RTOL, ETOL. */
su2double FSI_Conv[2]; /*!< \brief Values to check the convergence of the FSI problem. */

unsigned long idxIncrement; /*!< \brief Index of the current load increment */
su2double loadIncrement; /*!< \brief Coefficient that determines the amount of load which is applied. */
unsigned long idxIncrement = 0; /*!< \brief Index of the current load increment */
su2double loadIncrement = 1.0; /*!< \brief Coefficient that determines the amount of load which is applied. */

su2double WAitken_Dyn; /*!< \brief Aitken's dynamic coefficient. */
su2double WAitken_Dyn_tn1; /*!< \brief Aitken's dynamic coefficient in the previous iteration. */
Expand Down
19 changes: 17 additions & 2 deletions SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,7 @@ void CDiscAdjMultizoneDriver::Run() {
Add_External_To_Solution(iZone);
}
else {
/*--- If we restarted, Solution already has all contribution,
/*--- If we restarted, Solution already has all contributions,
* we run only one inner iter to compute the cross terms. ---*/
eval_transfer = true;
}
Expand All @@ -274,7 +274,12 @@ void CDiscAdjMultizoneDriver::Run() {
solver_container, numerics_container, config_container,
surface_movement, grid_movement, FFDBox, iZone, INST_0);

/*--- Print out the convergence data to screen and history file ---*/
/*--- This is done explicitly here for multizone cases, only in inner iterations and not when
* extracting cross terms so that the adjoint residuals in each zone still make sense. ---*/

Set_SolutionOld_To_Solution(iZone);

/*--- Print out the convergence data to screen and history file. ---*/

bool converged = iteration_container[iZone][INST_0]->Monitor(output_container[iZone], integration_container,
geometry_container, solver_container, numerics_container,
Expand All @@ -293,6 +298,7 @@ void CDiscAdjMultizoneDriver::Run() {
/*--- Extracting adjoints for solvers in jZone w.r.t. to the output of all solvers in iZone,
* that is, for the cases iZone != jZone we are evaluating cross derivatives between zones. ---*/

config_container[jZone]->SetInnerIter(0);
iteration_container[jZone][INST_0]->Iterate(output_container[jZone], integration_container, geometry_container,
solver_container, numerics_container, config_container,
surface_movement, grid_movement, FFDBox, jZone, INST_0);
Expand Down Expand Up @@ -898,6 +904,15 @@ void CDiscAdjMultizoneDriver::Add_External_To_Solution(unsigned short iZone) {
}
}

void CDiscAdjMultizoneDriver::Set_SolutionOld_To_Solution(unsigned short iZone) {

for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {
auto solver = solver_container[iZone][INST_0][MESH_0][iSol];
if (solver && solver->GetAdjoint())
solver->GetNodes()->Set_OldSolution();
}
}

void CDiscAdjMultizoneDriver::Update_Cross_Term(unsigned short iZone, unsigned short jZone) {

for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) {
Expand Down
7 changes: 4 additions & 3 deletions SU2_CFD/src/iteration_structure.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -473,11 +473,12 @@ void CFluidIteration::Iterate(COutput *output,
RUNTIME_RADIATION_SYS, val_iZone, val_iInst);
}

/*--- Adapt the CFL number using an exponential progression
with under-relaxation approach. ---*/
/*--- Adapt the CFL number using an exponential progression with under-relaxation approach. ---*/

if (config[val_iZone]->GetCFL_Adapt() == YES) {
solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone]);
SU2_OMP_PARALLEL
solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->AdaptCFLNumber(geometry[val_iZone][val_iInst],
solver[val_iZone][val_iInst], config[val_iZone]);
}

/*--- Call Dynamic mesh update if AEROELASTIC motion was specified ---*/
Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/src/solvers/CAdjEulerSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2301,7 +2301,7 @@ void CAdjEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **sol

for (iPoint = 0; iPoint < nPointDomain; iPoint++)
for (iVar = 0; iVar < nVar; iVar++) {
nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_AdjFlow()*LinSysSol[iPoint*nVar+iVar]);
nodes->AddSolution(iPoint,iVar, config->GetRelaxation_Factor_Adjoint()*LinSysSol[iPoint*nVar+iVar]);
}

/*--- MPI solution ---*/
Expand Down
2 changes: 0 additions & 2 deletions SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -667,8 +667,6 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co
}
}

if(multizone) nodes->Set_OldSolution();

SetResidual_RMS(geometry, config);

}
Expand Down
55 changes: 25 additions & 30 deletions SU2_CFD/src/solvers/CDiscAdjSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -476,26 +476,26 @@ void CDiscAdjSolver::SetAdj_ObjFunc(CGeometry *geometry, CConfig *config) {

void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config){

bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
bool multizone = config->GetMultizone_Problem();
const bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND;
const bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed;
const bool multizone = config->GetMultizone_Problem();

unsigned short iVar;
unsigned long iPoint;
su2double residual;
const su2double relax = (config->GetInnerIter()==0)? 1.0 : config->GetRelaxation_Factor_Adjoint();

/*--- Set Residuals to zero ---*/

for (iVar = 0; iVar < nVar; iVar++) {
SetRes_RMS(iVar,0.0);
SetRes_Max(iVar,0.0,0);
for (auto iVar = 0u; iVar < nVar; iVar++) {
SetRes_RMS(iVar,0.0);
SetRes_Max(iVar,0.0,0);
}

/*--- Set the old solution ---*/
/*--- Set the old solution and compute residuals. ---*/

if(!multizone) nodes->Set_OldSolution();

for (iPoint = 0; iPoint < nPoint; iPoint++) {
for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {

const su2double isdomain = (iPoint < nPointDomain)? 1.0 : 0.0;

/*--- Extract the adjoint solution ---*/

Expand All @@ -506,13 +506,22 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);
}

/*--- Store the adjoint solution ---*/
/*--- Relax and store the adjoint solution, compute the residuals. ---*/

nodes->SetSolution(iPoint,Solution);
for (auto iVar = 0u; iVar < nVar; iVar++) {
su2double residual = relax*(Solution[iVar]-nodes->GetSolution_Old(iPoint,iVar));
nodes->AddSolution(iPoint, iVar, residual);

residual *= isdomain;
AddRes_RMS(iVar,pow(residual,2));
AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
}
}

SetResidual_RMS(geometry, config);

if (time_n_needed) {
for (iPoint = 0; iPoint < nPoint; iPoint++) {
for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {

/*--- Extract the adjoint solution at time n ---*/

Expand All @@ -523,8 +532,9 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
nodes->Set_Solution_time_n(iPoint,Solution);
}
}

if (time_n1_needed) {
for (iPoint = 0; iPoint < nPoint; iPoint++) {
for (auto iPoint = 0u; iPoint < nPoint; iPoint++) {

/*--- Extract the adjoint solution at time n-1 ---*/

Expand All @@ -536,21 +546,6 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
}
}

/*--- Set the residuals ---*/

for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
for (iVar = 0; iVar < nVar; iVar++) {
residual = nodes->GetSolution(iPoint,iVar) - nodes->GetSolution_Old(iPoint,iVar);

AddRes_RMS(iVar,residual*residual);
AddRes_Max(iVar,fabs(residual),geometry->node[iPoint]->GetGlobalIndex(),geometry->node[iPoint]->GetCoord());
}
}

if(multizone) nodes->Set_OldSolution();

SetResidual_RMS(geometry, config);

}

void CDiscAdjSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) {
Expand Down
4 changes: 0 additions & 4 deletions SU2_CFD/src/solvers/CFEASolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,6 @@ CFEASolver::CFEASolver(bool mesh_deform_mode) : CSolver(mesh_deform_mode) {
Total_CFEA = 0.0;
WAitken_Dyn = 0.0;
WAitken_Dyn_tn1 = 0.0;
idxIncrement = 0;
loadIncrement = 1.0;

element_container = new CElement** [MAX_TERMS]();
for (unsigned short iTerm = 0; iTerm < MAX_TERMS; iTerm++)
Expand Down Expand Up @@ -124,8 +122,6 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() {
Total_CFEA = 0.0;
WAitken_Dyn = 0.0;
WAitken_Dyn_tn1 = 0.0;
idxIncrement = 0;
loadIncrement = 0.0;

SetFSI_ConvValue(0,0.0);
SetFSI_ConvValue(1,0.0);
Expand Down
60 changes: 46 additions & 14 deletions SU2_CFD/src/solvers/CSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2320,8 +2320,6 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
CSolver ***solver_container,
CConfig *config) {

/// TODO: Add OpenMP stuff to this method.

/* Adapt the CFL number on all multigrid levels using an
exponential progression with under-relaxation approach. */

Expand All @@ -2330,6 +2328,7 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
const su2double CFLFactorIncrease = config->GetCFL_AdaptParam(1);
const su2double CFLMin = config->GetCFL_AdaptParam(2);
const su2double CFLMax = config->GetCFL_AdaptParam(3);
const bool fullComms = (config->GetComm_Level() == COMM_FULL);

for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {

Expand Down Expand Up @@ -2365,6 +2364,9 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
/* Check that we are meeting our nonlinear residual reduction target
over time so that we do not get stuck in limit cycles. */

SU2_OMP_MASTER
{ /* Only the master thread updates the shared variables. */

Old_Func = New_Func;
unsigned short Res_Count = 100;
if (NonLinRes_Series.size() == 0) NonLinRes_Series.resize(Res_Count,0.0);
Expand Down Expand Up @@ -2408,17 +2410,30 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
Reset the array so that we delay the next decrease for some iterations. */

if (fabs(NonLinRes_Value) < 0.1*New_Func) {
reduceCFL = true;
NonLinRes_Counter = 0;
for (unsigned short iCounter = 0; iCounter < Res_Count; iCounter++)
NonLinRes_Series[iCounter] = New_Func;
}

} /* End SU2_OMP_MASTER, now all threads update the CFL number. */
SU2_OMP_BARRIER

if (fabs(NonLinRes_Value) < 0.1*New_Func) {
reduceCFL = true;
}

/* Loop over all points on this grid and apply CFL adaption. */

su2double myCFLMin = 1e30;
su2double myCFLMax = 0.0;
su2double myCFLSum = 0.0;
su2double myCFLMin = 1e30, myCFLMax = 0.0, myCFLSum = 0.0;

SU2_OMP_MASTER
if ((iMesh == MESH_0) && fullComms) {
Min_CFL_Local = 1e30;
Max_CFL_Local = 0.0;
Avg_CFL_Local = 0.0;
}

SU2_OMP_FOR_STAT(roundUpDiv(geometry[iMesh]->GetnPointDomain(),omp_get_max_threads()))
for (unsigned long iPoint = 0; iPoint < geometry[iMesh]->GetnPointDomain(); iPoint++) {

/* Get the current local flow CFL number at this point. */
Expand Down Expand Up @@ -2476,20 +2491,37 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
solverTurb->GetNodes()->SetLocalCFL(iPoint, CFL);
}

/* Store min and max CFL for reporting on fine grid. */
/* Store min and max CFL for reporting on the fine grid. */

myCFLMin = min(CFL,myCFLMin);
myCFLMax = max(CFL,myCFLMax);
myCFLSum += CFL;
if ((iMesh == MESH_0) && fullComms) {
myCFLMin = min(CFL,myCFLMin);
myCFLMax = max(CFL,myCFLMax);
myCFLSum += CFL;
}

}

/* Reduce the min/max/avg local CFL numbers. */

SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
if ((iMesh == MESH_0) && fullComms) {
SU2_OMP_CRITICAL
{ /* OpenMP reduction. */
Min_CFL_Local = min(Min_CFL_Local,myCFLMin);
Max_CFL_Local = max(Max_CFL_Local,myCFLMax);
Avg_CFL_Local += myCFLSum;
}
SU2_OMP_BARRIER

SU2_OMP_MASTER
{ /* MPI reduction. */
myCFLMin = Min_CFL_Local; myCFLMax = Max_CFL_Local; myCFLSum = Avg_CFL_Local;
SU2_MPI::Allreduce(&myCFLMin, &Min_CFL_Local, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
SU2_MPI::Allreduce(&myCFLMax, &Max_CFL_Local, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
SU2_MPI::Allreduce(&myCFLSum, &Avg_CFL_Local, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
Avg_CFL_Local /= su2double(geometry[iMesh]->GetGlobal_nPointDomain());
}
SU2_OMP_BARRIER
}

}

Expand Down
2 changes: 1 addition & 1 deletion TestCases/cont_adj_euler/wedge/inv_wedge_ROE.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 )
TIME_DISCRE_ADJFLOW= EULER_IMPLICIT
%
% Relaxation coefficient
RELAXATION_FACTOR_ADJFLOW= 1.0
RELAXATION_FACTOR_ADJOINT= 1.0
%
% Reduction factor of the CFL coefficient in the adjoint problem
CFL_REDUCTION_ADJFLOW= 0.8
Expand Down
Loading

0 comments on commit 53fb61f

Please sign in to comment.