Add 3D tests for proton boron fusion

[NeilZaim]

This PR is based on #2330, which is based on #2315, so it should probably be reviewed once these PRs have been merged.

This adds a 3D automated test for proton boron fusion (2D incoming in a follow-up PR). I tried to do all the checks I could think of, but that ended up being quite long, I'm sorry about that (at least, I'm now quite confident that the module works as expected 😄). If that is too hard to review, I can make the effort to split it into somewhat smaller PRs.

Since this PR is very long, this is the order that I'd recommend reading it:

1. The long comment at the start of the analysis script that explains the 5 different tests that are done.
2. The input file.
3. The analysis script, starting by the main function and reading the python functions as they are used.

For now I'm putting the [WIP] tag because I haven't tested this on GPU yet.
Edit: I've tested on GPU and it looks like the tests are working, so I'm removing the [WIP]

[lgtm-com]

This pull request introduces 1 alert when merging f1027f9 into fcf8872 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 1 alert when merging dd64205 into 02096d0 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 1 alert when merging 1f00d4c into 02096d0 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[lgtm-com]

This pull request introduces 1 alert when merging edf52f4 into 02096d0 - view on LGTM.com

new alerts:

• 1 for Except block handles 'BaseException'

[NeilZaim]

I see that the lgtm-bot raised an alert because I've added a try-except block. Is this not a good idea in an automated test? If so, I'll have to find a way to check if a species still exists in yt before reading the data.

[dpgrote]

The alert may be because it is a bare except:. It would be better to specify the error you are checking against. Perhaps except KeyError:?

[NeilZaim]

I think you are right, I've added the type of the exception and there was no alert raised.
Thanks for the hint!

[ax3l]

@NeilZaim ready for rebase :)

[ax3l]

This looks pretty solid to me and is well documented. Thank you! ✨

[ax3l]

the line continuation \ is not needed if the RHS expression is currently in an open bracket (, as in the cases here

[ax3l]

(optional change, just a note)

[ax3l]

This might need to be generalized, please see #2385 / #2683

[ax3l]

actually still works, as it only replaces the last 4 digits here to get the first step

[ax3l]

Current analysis python output in CI:

./analysis_proton_boron_fusion.py  Proton_Boron_Fusion_3D_plt00001
./analysis_proton_boron_fusion.py:303: RuntimeWarning: divide by zero encountered in true_divide
  return astrophysical_factor(E)/E_MeV * np.exp(-np.sqrt(E_Gamow_MeV / E_MeV))
./analysis_proton_boron_fusion.py:303: RuntimeWarning: invalid value encountered in multiply
  return astrophysical_factor(E)/E_MeV * np.exp(-np.sqrt(E_Gamow_MeV / E_MeV))
Traceback (most recent call last):
  File "./analysis_proton_boron_fusion.py", line 700, in <module>
    main()
  File "./analysis_proton_boron_fusion.py", line 690, in main
    eval("specific_check"+str(i)+"(data)")
  File "<string>", line 1, in <module>
  File "./analysis_proton_boron_fusion.py", line 629, in specific_check1
    check_alpha_yield(data, expected_fusion_number, E_com, proton_density = 1.,
  File "./analysis_proton_boron_fusion.py", line 400, in check_alpha_yield
    assert(np.all(is_close(alpha_weight_theory, alpha_weight_simulation,
AssertionError

[RemiLehe]

This looks good.
@NeilZaim Could you address the remaining inline comments?
Also, is it possible that some of these tests may fail sporadically because they rely on a random number generator? Or is the seed fixed here, so that the tests are reproducible?

[ax3l]

I see this in CI now:

./analysis_proton_boron_fusion.py:306: RuntimeWarning: divide by zero encountered in true_divide
  return astrophysical_factor(E)/E_MeV * np.exp(-np.sqrt(E_Gamow_MeV / E_MeV))
./analysis_proton_boron_fusion.py:306: RuntimeWarning: invalid value encountered in multiply
  return astrophysical_factor(E)/E_MeV * np.exp(-np.sqrt(E_Gamow_MeV / E_MeV))

[ax3l]

This one currently fails:

Traceback (most recent call last):
  File "./analysis_proton_boron_fusion.py", line 703, in <module>
    main()
  File "./analysis_proton_boron_fusion.py", line 693, in main
    eval("specific_check"+str(i)+"(data)")
  File "<string>", line 1, in <module>
  File "./analysis_proton_boron_fusion.py", line 632, in specific_check1
    check_alpha_yield(data, expected_fusion_number, E_com, proton_density = 1.,
  File "./analysis_proton_boron_fusion.py", line 403, in check_alpha_yield
    assert(np.all(is_close(alpha_weight_theory, alpha_weight_simulation,

[ax3l]

ping @NeilZaim as bespoken: did you fix a couple of div by zeros here already or can you push such a change please? :)

[NeilZaim]

Hi @ax3l and @RemiLehe. Thanks for the reviews and sorry for the long delay in responding.

Point by point response:

• The division by zero warning was expected, it's because when I wrote

    conditions = [E_com == 0, E_com > 0]
    choices = [0., pb_cross_section(E_com)*compute_relative_v_com(E_com)]
    sigma_times_vrel = np.select(conditions, choices)

The function pb_cross_section(0) was evaluated and then discarded. I've written a small fix to remove the warning (somewhat dirty so if there's a cleaner solution I'm interested)

• The failing test was because the analysis script did not yet take into account the high energy modification of the cross section done in Add high energy asymptotic fit for Proton-Boron total cross section #2408.

• Also, is it possible that some of these tests may fail sporadically because they rely on a random number generator? Or is the seed fixed here, so that the tests are reproducible?

Yes that's a very good point, a lot of the tolerances have been set empirically so we may definitely expect sporadic failures (even if we fix the seed a change in another part of the code that uses random number generation may still break the test). I've run a series of tests (~25) locally with various seeds and updated the tolerances accordingly. Hopefully that should make this CI test more robust in the future.