Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Newton-Raphson for NEMO temperature computations #1627

Merged
merged 14 commits into from
May 27, 2022
Merged

Newton-Raphson for NEMO temperature computations #1627

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
3115641
add in newton-raphson root finder
WallyMaier May 7, 2022
8d7762f
debugging issues
WallyMaier May 7, 2022
9cce43a
missing braces
WallyMaier May 7, 2022
274c34e
needs testing, but runs
WallyMaier May 8, 2022
8368b55
Update SU2_CFD/src/fluid/CSU2TCLib.cpp
WallyMaier May 8, 2022
faab298
constifying variables
WallyMaier May 8, 2022
6537a36
Initialize newton method with Tve
WallyMaier May 9, 2022
837bf4f
update regression values
WallyMaier May 9, 2022
6ba0643
fix serial regression format issue
WallyMaier May 9, 2022
c2e6445
Initialization fix
WallyMaier May 16, 2022
9e7dc31
Merge branch 'develop' into feature_NEMO_newton
WallyMaier May 16, 2022
9c0b4fe
Merge branch 'develop' into feature_NEMO_newton
WallyMaier May 19, 2022
ec73ef8
Merge branch 'develop' into feature_NEMO_newton
WallyMaier May 24, 2022
3588426
update regression values
WallyMaier May 24, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion SU2_CFD/include/fluid/CMutationTCLib.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ class CMutationTCLib : public CNEMOGas {
/*!
* \brief Compute translational and vibrational temperatures vector.
*/
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel) final;
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoE, su2double rhoEve, su2double rhoEvel, su2double Tve_old) final;

/*!
* \brief Get species molar mass.
Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/include/fluid/CNEMOGas.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,7 @@ class CNEMOGas : public CFluidModel {
/*!
* \brief Compute translational and vibrational temperatures vector.
*/
virtual vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel) = 0;
virtual vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel, su2double Tve_old) = 0;

/*!
* \brief Compute speed of sound.
Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/include/fluid/CSU2TCLib.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -186,7 +186,7 @@ class CSU2TCLib : public CNEMOGas {
/*!
* \brief Compute translational and vibrational temperatures vector.
*/
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel) final;
vector<su2double>& ComputeTemperatures(vector<su2double>& val_rhos, su2double rhoEmix, su2double rhoEve, su2double rhoEvel, su2double Tve_old) final;

private:

Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/src/fluid/CMutationTCLib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,7 @@ vector<su2double>& CMutationTCLib::GetThermalConductivities(){
return ThermalConductivities;
}

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

rhos = val_rhos;

Expand Down
83 changes: 60 additions & 23 deletions SU2_CFD/src/fluid/CSU2TCLib.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1986,9 +1986,8 @@ void CSU2TCLib::ThermalConductivitiesGY(){

}

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

vector<su2double> val_eves;
rhos = val_rhos;

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

/*--- Set temperature clipping values ---*/
su2double Tmin = 50.0; su2double Tmax = 8E4;
su2double Tve_o = 50.0; su2double Tve2 = 8E4;
const su2double Tmin = 50.0; const su2double Tmax = 8E4;
const su2double Tvemin = 50.0; const su2double Tvemax = 8E4;
su2double Tve_o = 50.0; su2double Tve2 = 8E4;

/* Determine if the temperature lies within the acceptable range */
if (T < Tmin) T = Tmin;
else if (T > Tmax) T = Tmax;
if (Tve_old < 1) Tve_old = T; //For first fluid iteration
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

During the first fluid iteration, the Tve_old is zero. This causes the newton solver to fail.

This line fixes that issue, and seemed to be least intrusive fix, though open to other ideas.

if (T < Tmin) T = Tmin; else if (T > Tmax) T = Tmax;
if (Tve_old<Tvemin) Tve_old = Tvemin; else if (Tve_old>Tvemax) Tve_old = Tvemax;

/*--- Set vibrational temperature algorithm parameters ---*/
su2double Btol = 1.0E-6; // Tolerance for the Bisection method
unsigned short maxBIter = 50; // Maximum Bisection method iterations
const su2double NRtol = 1.0E-6; // Tolerance for the Newton-Raphson method
const su2double Btol = 1.0E-6; // Tolerance for the Bisection method
const unsigned short maxBIter = 50; // Maximum Bisection method iterations
const unsigned short maxNIter = 50; // Maximum Newton-Raphson iterations
const su2double scale = 0.9; // Scaling factor for Newton-Raphson step
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have a sense of the stability of the method wrt to learning rate? 0.9 ~ 1, seems like you could get away without this parameter. Are there cases for which the method fails to converge for large learning rates, but converges for small ones? Would it make sense to have this be an adaptive parameter? (i.e. scale = 1, if NR does not converge, scale = 0.2, break?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From my limited test cases, I haven't seen any issue with stability. I think I can test once I find a case where we seen temperature non-convergence.

I think having the parameter adaptive could be interesting, though may not be necessary if 0.9 (or 1 is sufficient).


/*--- Execute a Newton-Raphson root-finding method for Tve ---*/
//Initialize solution
Tve = T;
Tve = Tve_old;

// Execute the root-finding method
bool Bconvg = false;
su2double rhoEve_t = 0.0;

for (unsigned short iIter = 0; iIter < maxBIter; iIter++) {
Tve = (Tve_o+Tve2)/2.0;
val_eves = ComputeSpeciesEve(Tve);
rhoEve_t = 0.0;
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
if (fabs(rhoEve_t - rhoEve) < Btol) {
Bconvg = true;
break;
bool NRconvg = false;
su2double rhoEve_t = 0.0, rhoCvve = 0.0;

/*--- Newton-Raphson Method --*/
for (unsigned short iIter = 0; iIter < maxNIter; iIter++) {
rhoEve_t = rhoCvve = 0.0;
const auto& val_eves = ComputeSpeciesEve(Tve);
const auto& val_cvves = ComputeSpeciesCvVibEle(Tve);

for (iSpecies = 0; iSpecies < nSpecies; iSpecies++){
rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
rhoCvve += rhos[iSpecies] * val_cvves[iSpecies];
}

/*--- Find the roots ---*/
su2double f = rhoEve - rhoEve_t;
su2double df = -rhoCvve;
Tve2 = Tve - (f/df)*scale;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to the question on learning rate above, do we know any properties about the minimization problem we are solving here that help us guarantee convergence, or maybe can help us demonstrate robustness?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it a backtracking method like e.g. the Armijo rule used for gradient descent/line search what you have in mind? Such a thing could at least prevent you from divergence.
Maybe adding a solver class for Newton-like methods to somewhere in Common/ could make sense at some point :)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree we could have a common folder of these methods, though Im not sure where else in the code they are used.


/*--- Check for convergence ---*/
if ((fabs(Tve2-Tve) < NRtol) && (Tve > Tvemin) && (Tve < Tvemax)) {
NRconvg = true;
Tve = Tve2;
break;
} else {
if (rhoEve_t > rhoEve) Tve2 = Tve;
else Tve_o = Tve;
Tve = Tve2;
}
}

// If the Newton-Raphson method has converged, assign the value of Tve.
// Otherwise, execute a bisection root-finding method
Tve_o = Tvemin; Tve2 = Tvemax;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is having the bisection as a backup still necessary? Is this how it was implemented in feature_TNE2? Have any of the test cases you tried required the use of the bisection to solve for value? Are there specific circumstances where bisection will work where NR will fail? The only examples I can think of are bad initialization or if there is some weirdness/non-smoothness in the function.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was how the original TNE2 was implemented, and generally doesn't seem to get activated except the first iteration.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So at the first iteration, NR fails and it switches to bisection, then it switches back to NR? That sounds like a problematic initialization.

Copy link
Contributor Author

@WallyMaier WallyMaier May 12, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, the issue is the first iteration where the primitive vector isn't set ie:
Tve_old = V[TVE_INDEX] in the solver is actually just doing Tve_old = 0.

So the two option in my head would be set the Primitive temperature early in the process, or have a boolean check. The boolean check works fine.

Copy link
Contributor Author

@WallyMaier WallyMaier May 16, 2022 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added the boolean above and there hasn't been a case where the bisection was needed.

if (!NRconvg) {
for (unsigned short iIter = 0; iIter < maxBIter; iIter++) {
Tve = (Tve_o+Tve2)/2.0;
const auto& val_eves = ComputeSpeciesEve(Tve);
rhoEve_t = 0.0;
for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) rhoEve_t += rhos[iSpecies] * val_eves[iSpecies];
if (fabs(rhoEve_t - rhoEve) < Btol) {
Bconvg = true;
break;
} else {
if (rhoEve_t > rhoEve) Tve2 = Tve;
else Tve_o = Tve;
}
}
}

// If absolutely no convergence, then assign to the TR temperature
if (!Bconvg) {
if (!NRconvg && !Bconvg ) {
Tve = T;
cout <<"Warning: temperatures did not converge, error= "<< fabs(rhoEve_t-rhoEve)<<endl;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you tried this on any cases where we previously didn't get convergence? Does it help with that, or does it just speed up the solver?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have a case with the Fire2 waiting to be run, I will report back on that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new temperature routine avoid the "did not converge" errors, but did not improve the solution quality.

}
Expand Down
3 changes: 2 additions & 1 deletion SU2_CFD/src/variables/CNEMOEulerVariable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,8 @@ bool CNEMOEulerVariable::Cons2PrimVar(su2double *U, su2double *V,
}

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

/*--- Temperatures ---*/
V[T_INDEX] = T[0];
Expand Down
8 changes: 4 additions & 4 deletions TestCases/parallel_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ def main():
ionized.cfg_dir = "nonequilibrium/thermalbath/finitechemistry"
ionized.cfg_file = "weakly_ionized.cfg"
ionized.test_iter = 10
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]
ionized.test_vals = [-29.806157, -11.130797, -11.337264, -17.235059, -17.578729, -15.190274, -25.070169, -32.000000, -5.174887, 0.000000, 0.000000]
ionized.su2_exec = "mpirun -n 2 SU2_CFD"
ionized.timeout = 1600
ionized.new_output = True
Expand All @@ -71,7 +71,7 @@ def main():
thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen"
thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg"
thermalbath_frozen.test_iter = 10
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.039251, -12.171781, -32.000000, 10.013545]
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -11.962477, -11.962477, -32.000000, 10.013545]
thermalbath_frozen.su2_exec = "mpirun -n 2 SU2_CFD"
thermalbath_frozen.timeout = 1600
thermalbath_frozen.new_output = True
Expand All @@ -83,7 +83,7 @@ def main():
invwedge.cfg_dir = "nonequilibrium/invwedge"
invwedge.cfg_file = "invwedge.cfg"
invwedge.test_iter = 10
invwedge.test_vals = [-1.042843, -1.567606, -18.300689, -18.628064, -18.574092, 2.275192, 1.879772, 5.319420, 0.873699]
invwedge.test_vals = [-1.042842, -1.567605, -18.299898, -18.627275, -18.573299, 2.275192, 1.879772, 5.319421, 0.873699]
invwedge.su2_exec = "mpirun -n 2 SU2_CFD"
invwedge.timeout = 1600
invwedge.new_output = True
Expand All @@ -95,7 +95,7 @@ def main():
visc_cone.cfg_dir = "nonequilibrium/axi_visccone"
visc_cone.cfg_file = "axi_visccone.cfg"
visc_cone.test_iter = 10
visc_cone.test_vals = [-5.222001, -5.746254, -20.569422, -20.633784, -20.547640, -1.928394, -2.246909, 1.255970, -3.208248]
visc_cone.test_vals = [-5.222267, -5.746522, -20.569425, -20.633787, -20.547644, -1.928714, -2.247316, 1.255761, -3.208360]
visc_cone.su2_exec = "mpirun -n 2 SU2_CFD"
visc_cone.timeout = 1600
visc_cone.new_output = True
Expand Down
6 changes: 3 additions & 3 deletions TestCases/serial_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ def main():
thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen"
thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg"
thermalbath_frozen.test_iter = 10
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.018022, -12.217291, -32.000000, 10.013545]
thermalbath_frozen.test_vals = [-32.000000, -32.000000, -12.018022, -11.978730, -32.000000, 10.013545]
thermalbath_frozen.su2_exec = "SU2_CFD"
thermalbath_frozen.timeout = 1600
thermalbath_frozen.new_output = True
Expand All @@ -69,7 +69,7 @@ def main():
invwedge.cfg_dir = "nonequilibrium/invwedge"
invwedge.cfg_file = "invwedge.cfg"
invwedge.test_iter = 10
invwedge.test_vals = [-1.046323, -1.571086, -18.300667, -18.628064, -18.574092, 2.271777, 1.875687, 5.315769, 0.870008]
invwedge.test_vals = [-1.046323, -1.571086, -18.299885, -18.627285, -18.573308, 2.271778, 1.875687, 5.315769, 0.870008]
invwedge.su2_exec = "SU2_CFD"
invwedge.timeout = 1600
invwedge.new_output = True
Expand All @@ -81,7 +81,7 @@ def main():
visc_cone.cfg_dir = "nonequilibrium/axi_visccone"
visc_cone.cfg_file = "axi_visccone.cfg"
visc_cone.test_iter = 10
visc_cone.test_vals = [-5.215288, -5.739428, -20.545050, -20.618702, -20.502532, -1.917680, -2.239596, 1.262771, -3.205521]
visc_cone.test_vals = [-5.215236, -5.739371, -20.545048, -20.618700, -20.502532, -1.917617, -2.239664, 1.262783, -3.205463]
visc_cone.su2_exec = "SU2_CFD"
visc_cone.timeout = 1600
visc_cone.new_output = True
Expand Down