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.

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

Conversation

WallyMaier
Copy link
Contributor

Proposed Changes

This replaces the bisection method in NEMO to compute the Tve values. The NR solver should be faster and offer some speed up, I believe the ComputeTemperatures function is the slowest in NEMO.

So I have some cases running to highlight the speedup (hopefully not a slowdown). From a quick glance, it seems the iteration time is very sensitive to the scaling factor.

Lastly, this will probably destroy the regression tests for NEMO, since our temperatures will be slightly different now.

Related Work

Resolve any issues (bug fix or feature request), note any related PRs, or mention interactions with the work of others, if any.

PR Checklist

  • I am submitting my contribution to the develop branch.
  • My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags, or simply --warnlevel=2 when using meson).
  • My contribution is commented and consistent with SU2 style.
  • I have added a test case that demonstrates my contribution, if necessary.
  • I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.

su2double Btol = 1.0E-6; // Tolerance for the Bisection method
unsigned short maxBIter = 50; // Maximum Bisection method iterations
unsigned short maxNIter = 50; // Maximum Newton-Raphson iterations
su2double scale = 0.5; // Scaling factor for Newton-Raphson step
Copy link
Contributor Author

Choose a reason for hiding this comment

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

This scaling factor of 0.5 would result in a slow down (~25%). But 0.9 seems to be a 20% speed increase.....

SU2_CFD/src/fluid/CSU2TCLib.cpp Outdated Show resolved Hide resolved
SU2_CFD/src/fluid/CSU2TCLib.cpp Outdated Show resolved Hide resolved
SU2_CFD/src/fluid/CSU2TCLib.cpp Outdated Show resolved Hide resolved
SU2_CFD/src/fluid/CSU2TCLib.cpp Outdated Show resolved Hide resolved

// 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

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

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.

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).

/*--- 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.

// 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.


/* 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.

@jtneedels jtneedels self-requested a review May 27, 2022 04:14
@WallyMaier WallyMaier merged commit 0e7844e into develop May 27, 2022
@WallyMaier WallyMaier deleted the feature_NEMO_newton branch May 27, 2022 04:40
@pr-triage pr-triage bot added the PR: merged label May 27, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants