Skip to content

Refactor code in fusion module to avoid machine-precision issues #3229

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

Merged
merged 3 commits into from
Jul 13, 2022
Merged

Refactor code in fusion module to avoid machine-precision issues #3229

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d91eda9
Refactor code in fusion module to avoid machine-precision issues
RemiLehe Jul 12, 2022
e0eaee4
Update formula
RemiLehe Jul 12, 2022
19411dd
Merge branch 'development' into nuclear_avoid_negative
RemiLehe Jul 12, 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
4 changes: 2 additions & 2 deletions Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -338,8 +338,8 @@ def E_com_to_p_sq_com(m1, m2, E):
## E is the total (kinetic+mass) energy of a two particle (with mass m1 and m2) system in
## its center of mass frame, in J.
## Returns the square norm of the momentum of each particle in that frame.
return E**2/(4.*scc.c**2) - (m1**2 + m2**2)*scc.c**2/2. + \
scc.c**6/(4.*E**2)*((m1**2 - m2**2)**2)
E_ratio = E/((m1+m2)*scc.c**2)
return m1*m2*scc.c**2 * (E_ratio**2 - 1) + (m1-m2)**2*scc.c**2/4 * (E_ratio - 1./E_ratio)**2

def compute_relative_v_com(E):
## E is the kinetic energy of proton+boron in the center of mass frame, in keV
Expand Down
17 changes: 9 additions & 8 deletions Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,6 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
const amrex::ParticleReal w_max = amrex::max(w1, w2);

constexpr auto one_pr = amrex::ParticleReal(1.);
constexpr auto inv_two_pr = amrex::ParticleReal(1./2.);
constexpr auto inv_four_pr = amrex::ParticleReal(1./4.);
constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
constexpr amrex::ParticleReal inv_csq = one_pr / ( c_sq );
Expand Down Expand Up @@ -107,7 +106,8 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq;

// Kinetic energy in the center of mass frame
const amrex::ParticleReal E_kin_star = std::sqrt(E_star_sq) - (m1 + m2)*c_sq;
const amrex::ParticleReal E_star = std::sqrt(E_star_sq);
const amrex::ParticleReal E_kin_star = E_star - (m1 + m2)*c_sq;

// Compute fusion cross section as a function of kinetic energy in the center of mass frame
auto fusion_cross_section = amrex::ParticleReal(0.);
Expand All @@ -118,14 +118,15 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part

// Square of the norm of the momentum of one of the particles in the center of mass frame
// Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
auto constexpr pow3 = [](double const x) { return x*x*x; };
const amrex::ParticleReal p_star_sq =
E_star_sq*inv_four_pr*inv_csq - (m1_sq + m2_sq)*c_sq*inv_two_pr +
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);
Copy link
Contributor

@NeilZaim NeilZaim Jul 12, 2022

Choose a reason for hiding this comment

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

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.

Copy link
Member Author

@RemiLehe RemiLehe Jul 12, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

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?

Copy link
Contributor

@NeilZaim NeilZaim Jul 13, 2022

Choose a reason for hiding this comment

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

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.

const amrex::ParticleReal p_star_sq = m1*m2*c_sq * ( pow2(E_ratio) - one_pr )
+ pow2(m1 - m2)*c_sq*inv_four_pr * pow2( E_ratio - 1._prt/E_ratio );

// Lorentz factors in the center of mass frame
const amrex::ParticleReal g1_star = std::sqrt(1._prt + p_star_sq / (m1_sq*c_sq));
const amrex::ParticleReal g2_star = std::sqrt(1._prt + p_star_sq / (m2_sq*c_sq));
const amrex::ParticleReal g1_star = std::sqrt(one_pr + p_star_sq / (m1_sq*c_sq));
const amrex::ParticleReal g2_star = std::sqrt(one_pr + p_star_sq / (m2_sq*c_sq));

// relative velocity in the center of mass frame
const amrex::ParticleReal v_rel = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) +
Expand Down
9 changes: 5 additions & 4 deletions Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionUtil.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,10 +108,11 @@ namespace {

// Square of the norm of the momentum of the products in the center of mass frame
// Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
auto constexpr pow3 = [](double const x) { return x*x*x; };
const amrex::ParticleReal p_star_f_sq =
E_star_f_sq*0.25_prt*inv_csq - (m1_out*m1_out + m2_out*m2_out)*c_sq*0.5_prt +
pow3(c_sq)*0.25_prt * pow2(m1_out*m1_out - m2_out*m2_out) / E_star_f_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 = std::sqrt(E_star_f_sq)/((m1_out + m2_out)*c_sq);
const amrex::ParticleReal p_star_f_sq = m1_out*m2_out*c_sq * ( pow2(E_ratio) - 1._prt )
+ pow2(m1_out - m2_out)*c_sq*0.25_prt * pow2( E_ratio - 1._prt/E_ratio );

// Compute momentum of first product in the center of mass frame, assuming isotropic
// distribution
Expand Down