Skip to content

Commit 0e7844e

Browse files
WallyMaierWallyMaier
authored
Merge pull request #1627 from su2code/feature_NEMO_newton
Newton-Raphson for NEMO temperature computations
2 parents 7a57336 + 3588426 commit 0e7844e

8 files changed

Lines changed: 73 additions & 35 deletions

File tree

‎SU2_CFD/include/fluid/CMutationTCLib.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -125,7 +125,7 @@ class CMutationTCLib : public CNEMOGas {
125125
/*!
126126
* \brief Compute translational and vibrational temperatures vector.
127127
*/
128-
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel) final;
128+
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel, su2double Tve_old) final;
129129

130130
/*!
131131
* \brief Get species molar mass.

‎SU2_CFD/include/fluid/CNEMOGas.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -177,7 +177,7 @@ class CNEMOGas : public CFluidModel {
177177
/*!
178178
* \brief Compute translational and vibrational temperatures vector.
179179
*/
180-
virtual vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel) = 0;
180+
virtual vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel, su2double Tve_old) = 0;
181181

182182
/*!
183183
* \brief Compute speed of sound.

‎SU2_CFD/include/fluid/CSU2TCLib.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -186,7 +186,7 @@ class CSU2TCLib : public CNEMOGas {
186186
/*!
187187
* \brief Compute translational and vibrational temperatures vector.
188188
*/
189-
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel) final;
189+
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel, su2double Tve_old) final;
190190

191191
private:
192192

‎SU2_CFD/src/fluid/CMutationTCLib.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -185,7 +185,7 @@ vector<su2double>& CMutationTCLib::GetThermalConductivities(){
185185
return ThermalConductivities;
186186
}
187187

188-
vector<su2double>& CMutationTCLib::ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel){
188+
vector<su2double>& CMutationTCLib::ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel, su2double Tve_old){
189189

190190
rhos = val_rhos;
191191

‎SU2_CFD/src/fluid/CSU2TCLib.cpp‎

Lines changed: 60 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1986,9 +1986,8 @@ void CSU2TCLib::ThermalConductivitiesGY(){
19861986

19871987
}
19881988

1989-
vector<su2double>& CSU2TCLib::ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel){
1989+
vector<su2double>& CSU2TCLib::ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel, su2double Tve_old) {
19901990

1991-
vector<su2double> val_eves;
19921991
rhos = val_rhos;
19931992

19941993
/*----------Translational temperature----------*/
@@ -2004,39 +2003,77 @@ vector<su2double>& CSU2TCLib::ComputeTemperatures(vector<su2double>& val_rhos, s
20042003
T = (rhoE - rhoEve - rhoE_f + rhoE_ref - rhoEvel) / rhoCvtr;
20052004

20062005
/*--- Set temperature clipping values ---*/
2007-
su2double Tmin = 50.0; su2double Tmax = 8E4;
2008-
su2double Tve_o = 50.0; su2double Tve2 = 8E4;
2006+
const su2double Tmin = 50.0; const su2double Tmax = 8E4;
2007+
const su2double Tvemin = 50.0; const su2double Tvemax = 8E4;
2008+
su2double Tve_o = 50.0; su2double Tve2 = 8E4;
20092009

20102010
/* Determine if the temperature lies within the acceptable range */
2011-
if (T < Tmin) T = Tmin;
2012-
else if (T > Tmax) T = Tmax;
2011+
if (Tve_old < 1) Tve_old = T; //For first fluid iteration
2012+
if (T < Tmin) T = Tmin; else if (T > Tmax) T = Tmax;
2013+
if (Tve_old<Tvemin) Tve_old = Tvemin; else if (Tve_old>Tvemax) Tve_old = Tvemax;
20132014

20142015
/*--- Set vibrational temperature algorithm parameters ---*/
2015-
su2double Btol = 1.0E-6; // Tolerance for the Bisection method
2016-
unsigned short maxBIter = 50; // Maximum Bisection method iterations
2016+
const su2double NRtol = 1.0E-6; // Tolerance for the Newton-Raphson method
2017+
const su2double Btol = 1.0E-6; // Tolerance for the Bisection method
2018+
const unsigned short maxBIter = 50; // Maximum Bisection method iterations
2019+
const unsigned short maxNIter = 50; // Maximum Newton-Raphson iterations
2020+
const su2double scale = 0.9; // Scaling factor for Newton-Raphson step
20172021

2022+
/*--- Execute a Newton-Raphson root-finding method for Tve ---*/
20182023
//Initialize solution
2019-
Tve = T;
2024+
Tve = Tve_old;
20202025

2021-
// Execute the root-finding method
20222026
bool Bconvg = false;
2023-
su2double rhoEve_t = 0.0;
2024-
2025-
for (unsigned short iIter = 0; iIter < maxBIter; iIter++) {
2026-
Tve = (Tve_o+Tve2)/2.0;
2027-
val_eves = ComputeSpeciesEve(Tve);
2028-
rhoEve_t = 0.0;
2029-
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
2030-
if (fabs(rhoEve_t - rhoEve) < Btol) {
2031-
Bconvg = true;
2032-
break;
2027+
bool NRconvg = false;
2028+
su2double rhoEve_t = 0.0, rhoCvve = 0.0;
2029+
2030+
/*--- Newton-Raphson Method --*/
2031+
for (unsigned short iIter = 0; iIter < maxNIter; iIter++) {
2032+
rhoEve_t = rhoCvve = 0.0;
2033+
const auto& val_eves = ComputeSpeciesEve(Tve);
2034+
const auto& val_cvves = ComputeSpeciesCvVibEle(Tve);
2035+
2036+
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++){
2037+
rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
2038+
rhoCvve += rhos[iSpecies] * val_cvves[iSpecies];
2039+
}
2040+
2041+
/*--- Find the roots ---*/
2042+
su2double f = rhoEve - rhoEve_t;
2043+
su2double df = -rhoCvve;
2044+
Tve2 = Tve - (f/df)*scale;
2045+
2046+
/*--- Check for convergence ---*/
2047+
if ((fabs(Tve2-Tve) < NRtol) && (Tve > Tvemin) && (Tve < Tvemax)) {
2048+
NRconvg = true;
2049+
Tve = Tve2;
2050+
break;
20332051
} else {
2034-
if (rhoEve_t > rhoEve) Tve2 = Tve;
2035-
else Tve_o = Tve;
2052+
Tve = Tve2;
20362053
}
20372054
}
2055+
2056+
// If the Newton-Raphson method has converged, assign the value of Tve.
2057+
// Otherwise, execute a bisection root-finding method
2058+
Tve_o = Tvemin; Tve2 = Tvemax;
2059+
if (!NRconvg) {
2060+
for (unsigned short iIter = 0; iIter < maxBIter; iIter++) {
2061+
Tve = (Tve_o+Tve2)/2.0;
2062+
const auto& val_eves = ComputeSpeciesEve(Tve);
2063+
rhoEve_t = 0.0;
2064+
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
2065+
if (fabs(rhoEve_t - rhoEve) < Btol) {
2066+
Bconvg = true;
2067+
break;
2068+
} else {
2069+
if (rhoEve_t > rhoEve) Tve2 = Tve;
2070+
else Tve_o = Tve;
2071+
}
2072+
}
2073+
}
2074+
20382075
// If absolutely no convergence, then assign to the TR temperature
2039-
if (!Bconvg) {
2076+
if (!NRconvg && !Bconvg ) {
20402077
Tve = T;
20412078
cout <<"Warning: temperatures did not converge, error= "<< fabs(rhoEve_t-rhoEve)<<endl;
20422079
}

‎SU2_CFD/src/variables/CNEMOEulerVariable.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -219,7 +219,8 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V,
219219
}
220220

221221
/*--- Assign temperatures ---*/
222-
const auto& T = fluidmodel->ComputeTemperatures(rhos, rhoE, rhoEve, 0.5*rho*sqvel);
222+
const su2double Tve_old = V[TVE_INDEX];
223+
const auto& T = fluidmodel->ComputeTemperatures(rhos, rhoE, rhoEve, 0.5*rho*sqvel, Tve_old);
223224

224225
/*--- Temperatures ---*/
225226
V[T_INDEX] = T[0];

‎TestCases/parallel_regression.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ def main():
5959
ionized.cfg_dir = "nonequilibrium/thermalbath/finitechemistry"
6060
ionized.cfg_file = "weakly_ionized.cfg"
6161
ionized.test_iter = 10
62-
ionized.test_vals = [-29.322805, -10.246260, -11.382786, -16.183264, -17.165896, -13.928855, -24.658131, -32.000000, -4.541637, 0.000000, 0.000000]
62+
ionized.test_vals = [-29.806157, -11.130797, -11.337264, -17.235059, -17.578729, -15.190274, -25.013626, -32.000000, -5.174887, 0.000000, 0.000000]
6363
ionized.su2_exec = "mpirun -n 2 SU2_CFD"
6464
ionized.timeout = 1600
6565
ionized.new_output = True
@@ -71,7 +71,7 @@ def main():
7171
thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen"
7272
thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg"
7373
thermalbath_frozen.test_iter = 10
74-
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.039251, -12.171781, -32.000000, 10.013545]
74+
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -11.962477, -11.962477, -32.000000, 10.013545]
7575
thermalbath_frozen.su2_exec = "mpirun -n 2 SU2_CFD"
7676
thermalbath_frozen.timeout = 1600
7777
thermalbath_frozen.new_output = True
@@ -83,7 +83,7 @@ def main():
8383
invwedge.cfg_dir = "nonequilibrium/invwedge"
8484
invwedge.cfg_file = "invwedge.cfg"
8585
invwedge.test_iter = 10
86-
invwedge.test_vals = [-1.042843, -1.567606, -18.300689, -18.628064, -18.574092, 2.275192, 1.879772, 5.319420, 0.873699]
86+
invwedge.test_vals = [-1.042842, -1.567605, -18.300680, -18.628055, -18.574084, 2.275192, 1.879772, 5.319421, 0.873699]
8787
invwedge.su2_exec = "mpirun -n 2 SU2_CFD"
8888
invwedge.timeout = 1600
8989
invwedge.new_output = True
@@ -95,7 +95,7 @@ def main():
9595
visc_cone.cfg_dir = "nonequilibrium/axi_visccone"
9696
visc_cone.cfg_file = "axi_visccone.cfg"
9797
visc_cone.test_iter = 10
98-
visc_cone.test_vals = [-5.222001, -5.746254, -20.569422, -20.633784, -20.547640, -1.928394, -2.246909, 1.255970, -3.208248]
98+
visc_cone.test_vals = [-5.222278, -5.746529, -20.569425, -20.633787, -20.547644, -1.928717, -2.247306, 1.255759, -3.208374]
9999
visc_cone.su2_exec = "mpirun -n 2 SU2_CFD"
100100
visc_cone.timeout = 1600
101101
visc_cone.new_output = True

‎TestCases/serial_regression.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ def main():
5757
thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen"
5858
thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg"
5959
thermalbath_frozen.test_iter = 10
60-
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.018022, -12.217291, -32.000000, 10.013545]
60+
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.018022, -11.978730, -32.000000, 10.013545]
6161
thermalbath_frozen.su2_exec = "SU2_CFD"
6262
thermalbath_frozen.timeout = 1600
6363
thermalbath_frozen.new_output = True
@@ -69,7 +69,7 @@ def main():
6969
invwedge.cfg_dir = "nonequilibrium/invwedge"
7070
invwedge.cfg_file = "invwedge.cfg"
7171
invwedge.test_iter = 10
72-
invwedge.test_vals = [-1.046323, -1.571086, -18.300667, -18.628064, -18.574092, 2.271777, 1.875687, 5.315769, 0.870008]
72+
invwedge.test_vals = [-1.046323, -1.571086, -18.300667, -18.628064, -18.574092, 2.271778, 1.875687, 5.315769, 0.870008]
7373
invwedge.su2_exec = "SU2_CFD"
7474
invwedge.timeout = 1600
7575
invwedge.new_output = True
@@ -81,7 +81,7 @@ def main():
8181
visc_cone.cfg_dir = "nonequilibrium/axi_visccone"
8282
visc_cone.cfg_file = "axi_visccone.cfg"
8383
visc_cone.test_iter = 10
84-
visc_cone.test_vals = [-5.215288, -5.739428, -20.545050, -20.618702, -20.502532, -1.917680, -2.239596, 1.262771, -3.205521]
84+
visc_cone.test_vals = [-5.215236, -5.739371, -20.545048, -20.618700, -20.502532, -1.917617, -2.239664, 1.262783, -3.205463]
8585
visc_cone.su2_exec = "SU2_CFD"
8686
visc_cone.timeout = 1600
8787
visc_cone.new_output = True

0 commit comments

Comments
 (0)