num test failures beyond bit-for-bit differences with intrinsic math

[warrickball]

If you try to build MESA with USE_CRMATH = NO in utils/makefile_header, some of the differences in num's test results are surprisingly large. You'll need to skip auto_diff's tests with touch auto_diff/skip_test.

e.g. for the Cash–Karp method, the expected output is

                                             calculated           1    1.7633919298826672D+00
                                             calculated           2   -8.3553976467637270D-01

whereas we get

                                             calculated           1    1.7634625993645798D+00
                                             calculated           2   -8.3545009921766566D-01

which differs already at the ~0.001 level. There are also various similarly large differences when calling the implicit solvers with numerical Jacobians (though differences seem to remain small with analytic Jacobians).

I'm personally not worried by differences within a few orders of magnitude of machine precision (say, <1e-10) but I don't understand why these are so different. There are very few calls to functions that are replaced by CRMATH and they seem to be confined to adjusting stepsizes anyway.

[rjfarmer]

I guess the first question is usually do we even use cash-karp? There's a lot in num and mtx we don't use anymore (or maybe ever)

[warrickball]

Cash–Karp isn't the only one but FWIW

$ grep cash_karp star/p*/*
star/private/hydro_rotation.f90:         call cash_karp_work_sizes(1,liwork,lwork)
star/private/hydro_rotation.f90:            call cash_karp( &

That call is in eval_fp_ft.

[adamjermyn]

Numerical jacobians will be computed with finite differences, which loses ~most precision right off the bat. Small differences then get magnified quickly...

…

-Adam

On Tue, Oct 05, 2021 at 9:54 AM, Warrick Ball ***@***.***> wrote: Cash–Karp isn't the only one but FWIW $ grep cash_karp star/p*/* star/private/hydro_rotation.f90: call cash_karp_work_sizes(1,liwork,lwork) star/private/hydro_rotation.f90: call cash_karp( & That call is in eval_fp_ft. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#331 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABPR5H7VA6PHUQPS3ESHKN3UFL7RJANCNFSM5FLV3ZYA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[fxt44]

the jacobian "only" guides the solution (e.g., how big a step can be taken for some accuracy criteria) and should not change the solution obtained. yes?

[warrickball]

Numerical jacobians will be computed with finite differences, which loses ~most precision right off the bat. Small differences then get magnified quickly...

On one hand, I totally agree. That's is also where I expect the problem is. On the other hand, even the finite difference Jacobians appear to only use things like +, -, *, / etc, rather than more troublesome functions like log, exp, pow, sin, cos, etc. I don't think CRMATH replaces those basic operations, so I expect the Jacobians to only differ to within a few orders of machine precision. Will that really be amplified?

the jacobian "only" guides the solution (e.g., how big a step can be taken for some accuracy criteria) and should not change the solution obtained. yes?

The tests still meet their tolerance requirements, though they're usually quite loose (e.g. 1e-4). My question is whether that should be the deciding limit when we switch from CRMATH to intrinsic operations.

[warrickball]

e.g., I presume (wrongly?) that CRMATH and intrinsic math both give the same imprecise difference of two similar floats, so even if that's what leads to an imprecise Jacobian, I still expect the same imprecise Jacobian.

[adamjermyn]

Right, but imagine if those two large numbers are different. Then the subtraction is doing the same thing, but it's amplifying a small difference still.

…

-Adam

On Tue, Oct 05, 2021 at 12:06 PM, Warrick Ball ***@***.***> wrote: e.g., I presume (wrongly?) that CRMATH and intrinsic math both give the same imprecise difference of two large numbers, so even if that's what leads to an imprecise Jacobian, I still expect the same imprecise Jacobian. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#331 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABPR5H3CPMQ6C3PUYRGQYZTUFMO6VANCNFSM5FLV3ZYA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[rjfarmer]

There are several calls to pow in mod_cash_karp.f

mesa/num/private/mod_cash_karp.f

Line 200 in 0434563

delta = SAFETY*pow(maxerr,-0.2d0)

Not sure how significant they are to the final result but still.

[orlox]

The calls to cash_karp in hydro_rotation come from the old way to compute the rotation corrections that is no longer the default. Kept those as an alternative when we introduced the new method which uses fits for the fp and ft corrections. By this point I say we can just strip that old functionality, the new stuff is working much better and no one has complained about it.

[warrickball]

I think I failed to provide enough context for what I'm getting at in this issue. I don't think it needs to be a priority and I certainly don't think we need to start changing things outside of num. Cash–Karp isn't the only problem but I'm not even sure any of the problematic cases are actually being used. We just ruled out Cash–Karp anyway. Diffusion uses implicit solvers but as far as I can see only with analytic Jacobians, which aren't a problem. It's probably been a while since we last ran the test suite with intrinsic math but I don't remember there being any catastrophes.

Regardless of whether any of these routines are being used, I think the first thing to establish is whether or not this level of failure in the num module is acceptable, which it might be. I in principle agree with @adamjermyn: anything involving finite differences of similar-sized numbers is likely to lose precision and amplify even small differences such that the limiting factor is the loose internal tolerances on these tests.

The reason I'm not immediately convinced by this is that these routines use very few of the functions that are (as far as I know) replaced by CRMATH. Even if x ≈ y, I expect to get the same (imprecise) x-y whether I use CRMATH or intrinsic. There are a few calls to pow sprinkled around:

$ grep -i pow\( num/p*/*
num/private/mod_cash_karp.f:            delta = SAFETY*pow(maxerr,-0.2d0)
num/private/mod_cash_karp.f:            delta = SAFETY*pow(maxerr,-0.25d0)
num/private/mod_dop853.f:      fac11=pow(err,expo1)
num/private/mod_dop853.f:      fac=fac11/pow(facold,beta)
num/private/mod_dop853.f:         h1=pow(0.01d0/der12,1.d0/iord) 
num/private/mod_dopri5.f:      fac11=pow(err,expo1)
num/private/mod_dopri5.f:      fac=fac11/pow(facold,beta)
num/private/mod_dopri5.f:         h1=pow(0.01d0/der12,1.d0/iord) 
num/private/mod_rosenbrock.f:      fac=max(fac2,min(fac1,pow(err,eloi)/safe))
num/private/mod_rosenbrock.f:               facgus=(hacc/h)*pow(err**2/erracc,eloi)/safe
num/private/mod_simplex.f90:                     weight = pow(weight,centroid_weight_power)

so maybe it's just those few calls (I think mostly involving the step adjustments) getting amplified by the Jacobians.

Happy to put this aside for a "rainy day".

[adamjermyn]

Note that those pow calls aren't to integer powers. Under the hood CRMATH implement those as calls to log, a multiplication, and a call to exp (more or less). So they're not any more innocuous than ordinary log calls!

[adamjermyn]

Just FYI I've removed cash_karp now that we no longer use it anywhere.

[adamjermyn]

Added labels 'rainy day' and 'good first issue' (subject to some experience with numerical analysis).