Refactor code in fusion module to avoid machine-precision issues #3229
Conversation
There was a problem hiding this comment.
Thanks for this PR. The formula looks good to me (note that there is a square missing at the end of the last parenthesis in the image).
I'm not yet 100% sure that the number can never be negative (see comment below), what do you think?
Also, be aware that the old formula is also currently used here: https://github.com/ECP-WarpX/WarpX/blob/32a5989d1765853a9bdac38d93c1c863790c65cc/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H#L109-L114
| pow3(c_sq) * inv_four_pr * pow2(m1_sq - m2_sq) / E_star_sq; | ||
| // The expression below is specifically written in a form that avoids returning | ||
| // small negative numbers due to machine precision errors, for low-energy particles | ||
| const amrex::ParticleReal E_ratio = E_star/((m1 + m2)*c_sq); |
There was a problem hiding this comment.
Do we have a complete guarantee that this number will never be strictly smaller than 1? E.g. if all the momenta are 0, this corresponds to the ratio between sqrt(E_rest**2) and E_rest so it's not obvious which of the two numbers will be smaller within machine precision.
There was a problem hiding this comment.
Agreed. This could in principle be negative.
However, I introduced this formula also in this PR #3153, where it did suppress an earlier failure due to invalid values.
(See this failure: https://dev.azure.com/ECP-WarpX/WarpX/_build/results?buildId=9168&view=logs&jobId=fc138a38-7e4f-5d28-e3f1-923e4d419807&j=fc138a38-7e4f-5d28-e3f1-923e4d419807&t=7e91da1c-aa6d-5e37-9c36-add2c6c921e8)
So my (very incomplete) understanding at this point is that the effect that you describe is maybe possible, but is maybe not as common as the previous failure.
If we encounter issues in the future with this new formula, we could certainly reconsider its expression. As an alternative, we could also intentionally enforce that E_star is larger than (m1 + m2)*c_sq, in the code, by e.g. setting E_star = max( E_star, (m1 + m2)*c_sq ). What do you think?
There was a problem hiding this comment.
That definitely makes sense to me.
I think I'm fine with both leaving it as it is (and reconsider if we encounter issues in the future) and with anticipating possible future issues by adding a line like E_star = max( E_star, (m1 + m2)*c_sq ).
Feel free to do what you think is best, I'll approve the PR in the meantime.
|
Thanks for checking the formula. I just updated it in the PR description. |
|
Also: thanks for pointing out the other part of the code where this is used. I updated the formula accordingly. |
…ST-WarpX#3229) * Refactor code in fusion module to avoid machine-precision issues * Update formula
4 participants
The formula
p_star_sq(which should always be a positive number) could occasionally return small negative numbers, due to machine-precision errors. (in the case where energy in the center-of-mass frame is very close to the mass energy of the two particles, i.e. when the particles have negligible kinetic energy in the center-of-mass frame). This is a problem because we later take the square root ofp_star_sqin the code.Therefore, this PR refactors this formula, so as to write it as the sum of 2 intrinsically positive numbers: