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 root find for Helmholtz EOS #304

Merged
merged 10 commits into from
Oct 3, 2023
Merged

Conversation

AlexHls
Copy link
Collaborator

@AlexHls AlexHls commented Sep 12, 2023

PR Summary

Adds a Newton-Raphson root find, specifically for use with the Helmholtz EOS. Although the current regula falsi produces on average more precise results, it is significantly more expensive. Moreover it is the method used in the originally published paper by Timmes.

Plot comparing both root finds

Overview

This benchmark was generated evaluating the EOS at random input values. It is probably not super scientific/ accurate, but should suffice the advantage of the Newton-Raphson root find.

A few points may (or may not) need to be addressed still:

  • The current implementation expects a different function to be passed compared to the other root finds, i.e. it expects a tuple containing the derivative. This works straightforwards with the Helmholtz EOS, but might not work in general. However since this is only an internal function call this should be fine.
  • I choose to make the root find method a runtime option, but I'm not sure if this is the way to go. So far I think it is better to have this as a runtime option so one can easily switch between the two methods, but I might be mistaken here.
  • The default is now the Newton-Raphson root find, purely based on the better performance.

PR Checklist

  • Adds a test for any bugs fixed. Adds tests for new features.
  • Format your changes by using the make format command after configuring with cmake.
  • Document any new features, update documentation for changes made.
  • Make sure the copyright notice on any files you modified is up to date.
  • After creating a pull request, note it in the CHANGELOG.md file
  • If preparing for a new release, update the version in cmake.

Sorry, something went wrong.

@AlexHls AlexHls added the enhancement New feature or request label Sep 12, 2023
@jhp-lanl
Copy link
Collaborator

Is there a possibility that the Newton method won't converge? If so, it might make sense to structure the root find option as use RF is NR fails or if the user specifically requests not to use NR.

@AlexHls
Copy link
Collaborator Author

AlexHls commented Sep 13, 2023

Is there a possibility that the Newton method won't converge? If so, it might make sense to structure the root find option as use RF is NR fails or if the user specifically requests not to use NR.

Thanks for the suggestion! I've added RF as fallback to the NR in cases the NR root find fails.

As far as convergence is concerned, this a bit more difficult to answer. In my naive test program, RF was actually more likely to not converge, but I know that in astrophysical simulations the NR also has trouble converging in certain parameter spaces. The question is if the RF would actually converge better in those cases (e.g. there is nothing that can be done once the root find runs into the edges of the underlying table) and, more importantly, if the non-convergence actually matters in these cases (vs the increased computational cost). To my knowledge, most astrophysical simulations have 'ignored' this and just bound the result of the root find back to the tabulated values. I might be mistaken though, but I guess this is something that one would first need to investigate extensively.

For now I have implemented RF as a fallback if NR does not converge, in my opinion this is a sensible default.

@jhp-lanl
Copy link
Collaborator

As far as convergence is concerned, this a bit more difficult to answer. In my naive test program, RF was actually more likely to not converge, but I know that in astrophysical simulations the NR also has trouble converging in certain parameter spaces. The question is if the RF would actually converge better in those cases (e.g. there is nothing that can be done once the root find runs into the edges of the underlying table) and, more importantly, if the non-convergence actually matters in these cases (vs the increased computational cost). To my knowledge, most astrophysical simulations have 'ignored' this and just bound the result of the root find back to the tabulated values. I might be mistaken though, but I guess this is something that one would first need to investigate extensively.

Interesting! I suppose I'm used to using the RF method on problems that have an easy way to bracket the solution a priori via some known bounds. In those cases, the RF method is guaranteed to converge albeit slowly sometimes. I assume that when you say that RF didn't converge, it means that you were unable to bracket the solution? Or was it simply converging so slowly that it ran out of iterations?

Another thing we could consider in the future is to improve the RF method to use either the Pegasus or Illinois algorithms (http://paulklein.se/newsite/teaching/rootfinding.pdf) to get superlinear convergence.

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
@AlexHls
Copy link
Collaborator Author

AlexHls commented Sep 14, 2023

As far as convergence is concerned, this a bit more difficult to answer. In my naive test program, RF was actually more likely to not converge, but I know that in astrophysical simulations the NR also has trouble converging in certain parameter spaces. The question is if the RF would actually converge better in those cases (e.g. there is nothing that can be done once the root find runs into the edges of the underlying table) and, more importantly, if the non-convergence actually matters in these cases (vs the increased computational cost). To my knowledge, most astrophysical simulations have 'ignored' this and just bound the result of the root find back to the tabulated values. I might be mistaken though, but I guess this is something that one would first need to investigate extensively.

Interesting! I suppose I'm used to using the RF method on problems that have an easy way to bracket the solution a priori via some known bounds. In those cases, the RF method is guaranteed to converge albeit slowly sometimes. I assume that when you say that RF didn't converge, it means that you were unable to bracket the solution? Or was it simply converging so slowly that it ran out of iterations?

Another thing we could consider in the future is to improve the RF method to use either the Pegasus or Illinois algorithms (http://paulklein.se/newsite/teaching/rootfinding.pdf) to get superlinear convergence.

Thanks for the interesting reference! I definitely want to take a deeper dive into this at some point.

In case of the RF failures, as far as I can tell the root find reaches the maximum number of iterations. (Note MAX_ITER_RF = 1000 vs MAX_ITER_NR = 100.

Ultimately I think that this needs a thorough investigation to make an informed decision / analysis on the advantages of different root finding algorithms. The point of this MR is mainly to make the NR root find available for the Helmholtz EOS, as it is the tried and testend method for this EOS (not to speak of the significant speedup, which is important for the applications I'm working with) that so far has not been questioned. In the end it is probably an extremely efficient algorithm for this use case since F. Timmes fine tuned the underlying table to work with NR (i.e. taking care that the derivatives are accurate).

@AlexHls AlexHls marked this pull request as ready for review September 14, 2023 09:24
@jhp-lanl
Copy link
Collaborator

The point of this MR is mainly to make the NR root find available for the Helmholtz EOS, as it is the tried and testend method for this EOS (not to speak of the significant speedup, which is important for the applications I'm working with) that so far has not been questioned. In the end it is probably an extremely efficient algorithm for this use case since F. Timmes fine tuned the underlying table to work with NR (i.e. taking care that the derivatives are accurate).

Seems fair! I'll take a closer look when I get a chance but @Yurlungur should also review

Copy link
Collaborator

@jhp-lanl jhp-lanl left a comment

Choose a reason for hiding this comment

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

Minor changes/questions

Comment on lines 240 to 246
// solves for f(x,params) - ytarget = 0
template <typename T>
PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget,
const Real guess, const Real a,
const Real b, const Real ytol, Real &xroot,
const RootCounts *counts = nullptr,
const bool &verbose = false) {
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think it's okay to use the same function signature as the regula falsi method, but I think you should at least add a comment about the behavior of the callable f being different (i.e. needing to return both function and derivative). Anybody using it now will quickly run into a compile time error if they try to just return the function, but why not save them the work of compiling to figure that out?

In the future I could see possibly wanting to develop an overload for this signature that uses the same behavior for f as the regula falsi and then wraps that callable in a function to perform a finite difference of some sort when an analytic derivative isn't available. Regardless of whether that happens though, comments will always help in understanding the subtle differences in template behavior between the root finding methods.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I agree with @jhp-lanl here

Comment on lines 268 to 282
// check if we are out of bounds
if (_x <= a && _xold <= a) {
if (verbose) {
printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess,
a, b);
}
break;
}
if (_x >= b && _xold >= b) {
if (verbose) {
printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess,
a, b);
}
break;
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

Maybe I'm missing the difference, but it looks like these conditions could be combined to make this code more DRY

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Thanks for the catch! This is a bit of a tricky one. In the original reference implementation, in each condition the root value was set to the appropriate boundary value, hence the two conditions. This has been moved to the Helmholtz call itself (since it is also necessary for the RF roots). I'm not sure if this is the smart way to do, i.e. if this bounding should be handled by the root finding methods or by the function calling the root find. For now I have decided that this should be handled by the function calling the root find (since one might not always want to bound the value to the boundary values) and added a comment in the appropriate location.

This is also the reason that this out-of-bounds case does not return a FAILURE status. One could argue that this is not really a successful root, but at least in the context of the Helmholtz EOS a root bound to the edges of the table is (at least in the reference implementation) taken as a success. Most likely this is a peculiarity for the Helmholtz EOS due to the fine tuning and won't translate to other application. I am however hesitant to implement a FAILURE return here since this would trigger an additional root find using RF. The bound result should be good enough, at least good enough as to not trigger a new attempt.

Copy link
Collaborator

Choose a reason for hiding this comment

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

You could add a fail_on_bound_root flag and default it to true and then just set it to false for the Helmholtz EOS. I'm a bit hesitant to have a root finder in the code whose default behavior could silently (if verbose isn't on) fail to find a root.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Fair point. Thanks for the suggestion! I've added the relevant option to the root find.

singularity-eos/base/root-finding-1d/root_finding.hpp Outdated Show resolved Hide resolved
Copy link
Collaborator

@Yurlungur Yurlungur left a comment

Choose a reason for hiding this comment

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

I approve. A few minor nitpicks/comments below that match @jhp-lanl 's thoughts.

Comment on lines 240 to 246
// solves for f(x,params) - ytarget = 0
template <typename T>
PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget,
const Real guess, const Real a,
const Real b, const Real ytol, Real &xroot,
const RootCounts *counts = nullptr,
const bool &verbose = false) {
Copy link
Collaborator

Choose a reason for hiding this comment

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

I agree with @jhp-lanl here

Comment on lines 269 to 273
if (_x <= a && _xold <= a) {
if (verbose) {
printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess,
a, b);
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

This is fine for now, but I strongly suspect this can be made much more robust with, e.g., a line search, or even a simple bounds enforcement. NR is prone to overshooting, and simply forcing it back in the domain is often enough to get it to converge.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I agree. There is most likely several ways to improve the convergence properties. For now the point is more to implement the NR as given by the reference implementation. I'm planning to conduct some more extensive benchmarking in the future to establish the performance vs. convergence in this context (since I suspect that this implementation is tuned for performance with just "good enough" convergence properties), but there is no timeline yet.

@AlexHls AlexHls requested a review from jhp-lanl October 2, 2023 19:59
Copy link
Collaborator

@jhp-lanl jhp-lanl left a comment

Choose a reason for hiding this comment

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

I approve, but @Yurlungur should take a final look before this gets merged.

@Yurlungur Yurlungur merged commit c186451 into main Oct 3, 2023
@Yurlungur Yurlungur deleted the hls/newton-raphson branch October 3, 2023 18:18
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants