Proposal: Add BaikalETK to the Toolkit

[Zach Etienne]

Update (May 2020): The original description below is a bit outdated, as BaikalETK was restructured (for the better!) a bit after ETK telecon discussions. For one, BaikalETK is now two thorns: Baikal and BaikalVacuum, where the former includes stress-energy source terms and the latter does not. Also the former contains 2nd and 4th order finite-difference derivative C code kernels and the latter includes 4th, 6th, and 8th order kernels. These may be specified at runtime, so the user need not rerun NRPy+ to generate the most commonly used kernels. Still for maximum flexibility one would need to run NRPy+ to generate kernels that are not already provided. For more details please read the documentation linked to below, and do check out the contents of the doc/ directories of each thorn.

BaikalETK is an open-source NRPy+-based BSSN thorn, which has been demonstrated to yield excellent agreement with ML_BSSN when evolving BBHs and BNSs, and in certain cases (e.g., ADM mass volume integral over entire grid vs time in a pre-merged BNS simulation) better results than ML_BSSN.

BaikalETK is BSD 2-clause open-source licensed and available in the NRPy+ github page (https://github.com/zachetienne/nrpytutorial). One can browse its source code & documentation with nbviewer: https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb

BaikalETK has the following neat features:

1. The core BSSN kernel is generated by NRPy+/SymPy, which is entirely free and open source, just like BaikalETK. According to benchmarks by Erik, BaikalETK runs more than 10% faster than the Kranc/Mathematica -generated ML_BSSN currently part of the Toolkit. Work is ongoing to improve BaikalETK's performance even more.
2. It is documented in a Jupyter notebook, which itself is linked to other Jupyter notebooks that document the BSSN formalism & gauge choices.
3. It can be generated by either running the aforementioned Jupyter notebook or just importing a Python module within NRPy+. The Jupyter notebook contains a self-validation test that will fail if the Python module is updated without a corresponding update to the Jupyter notebook.
4. Its scheduling and general structure are closely based on the very nice Lean code.
5. It supports arbitrary finite-differencing order and a variety of gauge conditions, which can be extended further. It implements Brown’s covariant formulation of the BSSN equations.
6. It supports as an option nonzero BSSN source terms, with appropriate linkage to TmunuBase.

Possible complications:

1. Unlike ML_BSSN, which (due to Mathematica licensing constraints and consistency with other monolithic thorns) needed to contain all gauge choices, finite difference orders, etc., BaikalETK will only generate a BSSN thorn with a particular gauge and finite-differencing order chosen at the time of thorn generation. That is to say, choices about finite difference order, etc must be made prior to compile-time, which is not conventional at the moment.

2. BaikalETK is not a complete ML_BSSN replacement yet as:

   1. BaikalETK does not yet support symmetries, but this should be rather straightforward to fix in an automated way, as NRPy+ requires information about gridfunctions and gridfunction symmetries.
   2. BaikalETK does not support Llama grids, and I’m not sure the path forward in this sense

‌

[Ian Hinder]

Hi Zach,

I was interested to read about BaikalETK. It's great to see development in this area, and while of course I have quite an attachment to Kranc (and less so to ML_BSSN), it's good to see a modern automatically generated BSSN solver that doesn't rely on proprietary software!

BaikalETK is an open-source NRPy+-based BSSN thorn, which has been demonstrated to yield excellent agreement with ML_BSSN when evolving BBHs and BNSs, and in certain cases (e.g., ADM mass volume integral over entire grid vs time in a pre-merged BNS simulation) better results than ML_BSSN. BaikalETK is BSD 2-clause open-source licensed and available in the NRPy+ github page (https://github.com/zachetienne/nrpytutorial). One can browse its source code & documentation with nbviewer: https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb BaikalETK has the following neat features: • The core BSSN kernel is generated by NRPy+/SymPy, which is entirely free and open source, just like BaikalETK. According to benchmarks by Erik, BaikalETK runs more than 10% faster than the Kranc/Mathematica -generated ML_BSSN currently part of the Toolkit. Work is ongoing to improve BaikalETK's performance even more.

Do you know whether this is because the same equations are implemented more efficiently, or because the formulation is inherently more efficient? I'd be interested to hear more about Erik's benchmark results.

Of course, in production BBH simulations, you usually see a depressingly small percentage of the runtime in the BSSN solver, but that's a topic for a different discussion!

Possible complications: • Unlike ML_BSSN, which (due to Mathematica licensing constraints and consistency with other monolithic thorns) needed to contain all gauge choices, finite difference orders, etc., BaikalETK will only generate a BSSN thorn with a particular gauge and finite-differencing order chosen at the time of thorn generation. That is to say, choices about finite difference order, etc must be made prior to compile-time, which is not conventional at the moment.

I didn't understand this. What sort of Mathematica licensing constraints would affect the requirement for ML_BSSN to contain all gauge choices? I thought that this was a performance-based decision by Erik in the early days. The idea was that you want to be able to use quantities computed already in the other RHSs also in the gauge condition, so the gauge condition needs to be in the same loop. Since you don't want to have multiple main evolution RHSs, one for each gauge condition, this means that all gauge conditions need to be in the main evolution loop. The Penn State BSSN code that I wrote, which was the evolution of the code I wrote for my PhD, had the gauge conditions as separate calculations, and could easily have split them out into other thorns if necessary. For finite difference orders, the original ML_BSSN used to have different thorns for the different orders. It was a relatively "recent" innovation to be able to generate multiple finite difference orders within the same thorn in Kranc, which is obviously a lot cleaner.

• BaikalETK is not a complete ML_BSSN replacement yet as: • BaikalETK does not yet support symmetries, but this should be rather straightforward to fix in an automated way, as NRPy+ requires information about gridfunctions and gridfunction symmetries. • BaikalETK does not support Llama grids, and I’m not sure the path forward in this sense

Essentially, for Llama, you need to support symmetries first, as Llama interpatch boundaries are implemented as symmetries. Then, you need to replace all derivatives with derivatives multiplied by the Jacobian matrix provided by the Coordinates thorn. You can have a look at the generated ML_BSSN code to see how this can be done.

-- Ian Hinder Research Software Engineer University of Manchester, UK

[Erik Schnetter]

On Thu, Nov 7, 2019 at 8:16 PM Ian Hinder [e-mail] wrote:

I compared three codes: (1) ML_BSSN as in the ET, (2) BaikalETK, and (3) a newer version of ML_BSSN that uses explicit loop tiling, inspired by our work on DG methods. I found that (2) was faster than (1) by about 10%-20%, and that (3) was faster than (2) by about 10%-20% in the one benchmark I ran (on Intel Skylake). I assume (but I did not check) that BaikalETK splits the loops differently than ML_BSSN, which is then faster either due to instruction or data cache effects.

I have shared my findings with others (but apparently not on the ET mailing list – sorry!), and I now think that this loop tiling is the way to go. It is also significantly faster for hydro, and improves performance on Intel's KNL. I've sketched the implementation to Zach, and I hope that we can have an even faster BaikalETK in the future. I'm also looking forward to advertising some of the optimizations we built into Kranc to Zach – he has a similar optimization stage already, so we'll have to see whether important ones are missing.

[helvi witek]

Adding my original report to Zach and our respective replies:

Report on Baikal and BaikalVacuum (13 May 2020)

• Two thorns "Baikal" (Einstein's equations coupled to Tmn) and "BaikalVacuum" (Einstein's equations in vacuum)
• very clean and readable implementation, both the the jupyter notebook and the auto-generated code.
• exceptionally well documented
• benchmarked against McLachlan (and Lean test is underway): reproduces standard qc0 run (BBH with 1 orbit before merger)

Baikal (13 May 2020):

• README:

  • "...in parfile_examples/" should be "... in par/" to be consistent with source code structure

• interface.ccl

  • definition of BSSN grid functions: you might consider using tags for the tensortype and tensorweight.
    (Q: how does it impact the evolution or is it only information?).
    Doing so would require to split the evolution variables into groups of same type and weight.

• param.ccl:

  • [CHECK:] no floor_value for conformal factor or lapse -- how is the zero at the puncture treated?

• schedule.ccl:

  • Why do you store the ADMBase gridfunctions (again) in Baikal? This is done in ADMBase/schedule.ccl, so why is it necessary here?
  • [CHECK:] BSSN-to-ADM conversion before calculating Ricci and rhs. Why?

• src/Symmetry_registration_oldCartGrid3D.c:

  • [CORRECTION:] default symmetry setting is hardwired that of a scalar in every tensor component.
    This should be corrected, otherwise using symmetries to reduce run costs might yield incorrect results.
    For example, gammax has symmetry (-1,1,1).

• src/BSSN_to_ADM.c

  • [CHECK:] trafo is $\gamma_{xx} = \frac{1}{cf^2} (\tilde{\gamma}_{xx}+1 )$ (why "+1"?) and same for yy and zz components

• src/BoundaryCondition_NewRad.c

  • [CORRECTION:] notation for NewRad is
    NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower)
    where var0 is the background value of the function. For $\alpha$, and diagonal components of $\tilde_{\gamma}_{ij}$ this is 1.
    In "src/BoundaryCondition_NewRad.c" they are registered as $0.0$ and should be corrected.

• [COSMETICS:] src/ADM_to_BSSN.c and other files: typo in comment, l.58 of source code: "Step 2 of 1"

BaikalVacuum (15 May 2020):

• same comments as above, except for schedule.ccl (BSSN->ADM conversion seems correct)
• you might consider adding a documenting notebook that you used to generate the vacuum code (similar to the "doc" folder in Baikal)
• [COSMETICS] param.ccl: order finite difference order from lowest to highest
• [OK] in "par": what is difference between bench-qc0-...par and qc0-...par?
• [OK] example parameter files do not use symmetry? Why? not implemented or not used?

NRPy+ Jupyter notebook (15 May 2020)

• in Step 2, before input [4]: link "BSSN.BSS_RHSs", "BSSN.BSSN_gauge_RHSs" and "finite_difference" lead to "404: not found" error.
  Please either fix it or revert to normal text. The tutorials load fine.
• Eq. (1) in BSSN_RHS tutorial (=11a in 1712.07658), i.e., evolution equation for $\bar{\gamma}{ij}$:
  you include a term $\bar{A}^{k}{k}$, i.e., the trace of the tracefree part of the extrinsic curvature.
  Should this not be zero by construction, or is it a way to enforce the trA = 0 constraint?
• What is the default choice for the reference metric, i.e., in the pre-generated Baikal/BaikalVacuum thorns?
  Do you state it in the jupyter notebook/source code and could point me to the information? If not, it may be useful to add.
  Btw, do you allow for different options at the level of code generation or should the user add it themselves?
• Step 3.c (schedule.ccl): see above
• Step 4.a, input [18]: check for symmetries (see above).

‌

[helvi witek]

[Zach’s response]

Hi Helvi and Roland.

Here is my response on the Baikal ETK thorns "Baikal" and "BaikalVacuum".

Thanks again Helvi for this amazing review. The Baikal ETK thorns are much better as a result! I believe I have completed all requested changes and made all appropriate fixes. Please see comments below.

Report on Baikal ETK (13 May 2020)
--------------------------------
- Two thorns "Baikal" (Einstein's equations coupled to Tmn) and "BaikalVacuum" (Einstein's equations in vacuum)
- very clean and readable implementation, both the the jupyter notebook and the auto-generated code.
- exceptionally well documented
- benchmarked against McLachlan (and Lean test is underway): reproduces standard qc0 run (BBH with 1 orbit before merger)

BaikalVacuum (15 May 2020):
---------------
- same comments as above, except for schedule.ccl (BSSN->ADM conversion seems correct)

Yep, see responses to Baikal.

- you might consider adding a documenting notebook that you used to generate the vacuum code (similar to the "doc" folder in Baikal)

Both thorns are generated using the same notebook, simultaneously. Regardless this is a good point (BaikalVacuum currently has no doc directory), so I added documentation.tex files along the lines Roland suggested (include a short description and links to the full documentation).

- [OK] in "par": what is difference between bench-qc0-...par and qc0-...par?

bench-* is for benchmarking purposes only; it only outputs to a bare minimum of files to avoid file I/O from causing noisy benchmarks. qc0-*.par is for science, outputting all the usual diagnostics files.

- [OK] example parameter files do not use symmetry? Why? not implemented or not used?

Symmetries aren't yet fully tested so I cannot vouch that they are 100% functional. I have added documentation to this effect in the thorns' doc/documentation.tex files.

NRPy+ Jupyter notebook (15 May 2020)
-------------------------------------
- in Step 2, before input [4]:  link "BSSN.BSS_RHSs", "BSSN.BSSN_gauge_RHSs" and "finite_difference" lead to "404: not found" error.
  Please either fix it or revert to normal text. The tutorials load fine.

These links work fine in the live Jupyter session, and point to a syntax-highlighed (editable) version of the Python file. Indeed in the read-only nbviewer (https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb),,) the editor is disabled and the user will find a 404 error instead. I updated the links to just open the files as read-only text files (without syntax highlighting, oh well).

- Eq. (1) in BSSN_RHS tutorial (=11a in 1712.07658), i.e., evolution equation for $\bar{\gamma}_{ij}$:
  you include a term $\bar{A}^{k}_{k}$, i.e., the trace of the tracefree part of the extrinsic curvature.
  Should this not be zero by construction, or is it a way to enforce the trA = 0 constraint?

This is how we enforce trA=0. I have added comments to this effect to the BSSN_RHSs tutorial, with an equation reference to Brown 2009:

First, note that the term containing $\bar{A}{k}^{k}$ is there to drive $\bar{A}{k}^{k}$ to zero. It was implemented in this form in T. Baumgarte's BSSN-in-spherical-coordinates code, against which NRPy+ and SENR were originally validated. You will find this term in Brown's covariant formulation of the BSSN time-evolution equations (see third term in Eq 21a).

- What is the default choice for the reference metric, i.e., in the pre-generated Baikal/BaikalVacuum thorns?
  Do you state it in the jupyter notebook/source code and could point me to the information? If not, it may be useful to add.

Sure, we use the Cartesian reference metric, since we're solving BSSN in Cartesian coordinates. This info is in the subtitle of the notebook "This notebook generates Baikal and BaikalVacuum, Einstein Toolkit thorns for solving Einstein's equations in the BSSN formalism, in Cartesian coordinates." (added emphasis). Support for other reference metrics would require appropriate boundary condition support e.g., provided by CarpetSlab. Note that V. Mewes' SphericalNR code does just this (with NRPy), but it's closed source at the moment.

  Btw, do you allow for different options at the level of code generation or should the user add it themselves?

Yes, NRPy+ supports a few other gauge conditions and arbitrary finite-differencing order (so long as the order is even). It's up to the user to extend as needed. :)

- Step 3.c (schedule.ccl): see comments on Baikal

See responses below.

- Step 4.a, input [18]: check for symmetries (see comments on Baikal).

Yep, see responses below.

Baikal (13 May 2020):
--------
- README:
  - "...in parfile_examples/" should be "... in par/" to be consistent with source code structure

Fixed and pushed (both thorns).

- interface.ccl
  - definition of BSSN grid functions: you might consider using tags for the tensortype and tensorweight.
    (Q: how does it impact the evolution or is it only information?).
    Doing so would require to split the evolution variables into groups of same type and weight.

My understanding is that tensortype & tensorweight parameters are only used by Llama for multiblock simulations. As stated in Baikal*/doc/documentation.tex, we don't yet support Llama/multiblock.

- param.ccl:
  - [CHECK:] no floor_value for conformal factor or lapse -- how is the zero at the puncture treated?

Fixed and pushed (both thorns). See floor_the_lapse.c; the function floor_the_lapse() does the dirtywork.

- schedule.ccl:
  - Why do you store the ADMBase gridfunctions (again) in Baikal? This is done in ADMBase/schedule.ccl, so why is it necessary here?

Fixed according to Roland's suggestion (both thorns). As we write to ADM gridfunctions, it is necessary to at least store 1 timelevel to ensure we're not writing to unallocated memory. Allocating 1 timelevel here ensures that this will not happen, while at the same time it allows for other thorns (e.g., ADMBase) to request allocation of more timelevels of storage. As Roland told me, Cactus always allocates the largest requested number of timelevels to a given gridfunction group if different thorns request different number of timelevels.

  - [CHECK:] BSSN-to-ADM conversion before calculating Ricci and rhs. Why?

Fixed! Turns out this was due to a bug in IllinoisGRMHD. I fixed IllinoisGRMHD and was able to remove this conversion.

- src/Symmetry_registration_oldCartGrid3D.c:
  - [CORRECTION:] default symmetry setting is hardwired that of a scalar in every tensor component.
                  This should be corrected, otherwise using symmetries to reduce run costs might yield incorrect results.
                  For example, gammax has symmetry (-1,1,1).

Fixed! (both thorns)

- src/BSSN_to_ADM.c
  - [CHECK:] trafo is $\gamma_{xx} = \frac{1}{cf^2} (\tilde{\gamma}_{xx}+1 )$ (why "+1"?) and same for yy and zz components
- src/BoundaryCondition_NewRad.c
  - [CORRECTION:] notation for NewRad is
                NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower)
        where var0 is the background value of the function. For $\alpha$, and diagonal components of $\tilde_{\gamma}_{ij}$ this is 1.
        In "src/BoundaryCondition_NewRad.c" they are registered as $0.0$ and should be corrected.

The two comments are actually related; we evolve h_{ij} = gamma_{ij} - delta_{ij}, where delta_{ij} is the Kronecker delta. Therefore the value of h_{ij} should indeed be zero at infinity.

Regarding the lapse having the wrong value at infinity -- nice catch! I have fixed and pushed that (both thorns).

- [COSMETICS:] src/ADM_to_BSSN.c and other files: typo in comment, l.58 of source code: "Step 2 of 1"

This has been fixed and pushed (both thorns)!

P.S.,

Per Roland's suggestion, I also added the NRPy+ git hash used to generate the thorns in each thorn's README file.

Thanks again everyone!

[helvi witek]

Updated report -- response to Zach's response (18 May 2020)

SUMMARY: All my comments have been addressed and corrected. So, I'm giving my OK to include the code.

**** BaikalVacuum
** Yep, see responses to Baikal.

HW: OK

‌

** Both thorns are generated using the same notebook, simultaneously. Regardless this is a good point
** (BaikalVacuum currently has no doc directory), so I added documentation.tex files along the lines
** Roland suggested (include a short description and links to the full documentation).

HW: OK

** bench-* is for benchmarking purposes only; it only outputs to a bare minimum of files to avoid
** file I/O from causing noisy benchmarks. qc0-*.par is for science, outputting all the usual diagnostics files.

HW: that's what I guessed. It may be worthwhile noting this in the script or README.

** Symmetries aren't yet fully tested so I cannot vouch that they are 100% functional.
** I have added documentation to this effect in the thorns' doc/documentation.tex files.

HW: OK

‌

**** NRPy+ Jupyter notebook (15 May 2020)

** These links work fine in the live Jupyter session, and point to a syntax-highlighed (editable) version of the Python file.
** Indeed in the read-only nbviewer (https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb),,)
** the editor is disabled and the user will find a 404 error instead. I updated the links to just open the files as read-only text files
** (without syntax highlighting, oh well).

HW: OK

‌

** This is how we enforce trA=0. I have added comments to this effect to the BSSN_RHSs tutorial, with an equation reference to Brown 2009:
** First, note that the term containing $\bar{A}{k}^{k}$ is there to drive $\bar{A}{k}^{k}$ to zero. It was implemented in this form in
** T. Baumgarte's BSSN-in-spherical-coordinates code, against which NRPy+ and SENR were originally validated. You will find this
** term in Brown's covariant formulation of the BSSN time-evolution equations (see third term in Eq 21a).

HW: perfect. Thanks! OK

‌

** Sure, we use the Cartesian reference metric, since we're solving BSSN in Cartesian coordinates.
** This info is in the subtitle of the notebook "This notebook generates Baikal and BaikalVacuum, Einstein Toolkit thorns for solving
** Einstein's equations in the BSSN formalism, in Cartesian coordinates." (added emphasis). Support for other reference metrics would
** require appropriate boundary condition support e.g., provided by CarpetSlab. Note that V. Mewes' SphericalNR code does just this (with NRPy),
** but it's closed source at the moment.

HW: Thanks for the highlight. OK.

** Yes, NRPy+ supports a few other gauge conditions and arbitrary finite-differencing order (so long as the order is even).
** It's up to the user to extend as needed. :)

HW: OK :)

**** - Step 3.c (schedule.ccl): see comments on Baikal
** See responses below.

HW: OK

‌

**** - Step 4.a, input [18]: check for symmetries (see comments on Baikal).
** Yep, see responses below.

HW: OK

‌

**** - README:
**** - "...in parfile_examples/" should be "... in par/" to be consistent with source code structure
** Fixed and pushed (both thorns).

HW: OK

** My understanding is that tensortype & tensorweight parameters are only used by Llama for multiblock simulations.
** As stated in Baikal*/doc/documentation.tex, we don't yet support Llama/multiblock.

HW: OK

** Fixed and pushed (both thorns). See floor_the_lapse.c; the function floor_the_lapse() does the dirtywork.

HW: OK. Thanks!

** Fixed according to Roland's suggestion (both thorns). As we write to ADM gridfunctions, it is necessary to at least store 1 timelevel
** to ensure we're not writing to unallocated memory. Allocating 1 timelevel here ensures that this will not happen, while at the same
** time it allows for other thorns (e.g., ADMBase) to request allocation of more timelevels of storage. As Roland told me, Cactus always
** allocates the largest requested number of timelevels to a given gridfunction group if different thorns request different number of timelevels.

HW: OK. Thanks for the clarification.

‌

**** - [CHECK:] BSSN-to-ADM conversion before calculating Ricci and rhs. Why?
** Fixed! Turns out this was due to a bug in IllinoisGRMHD. I fixed IllinoisGRMHD and was able to remove this conversion.

HW: OK

‌

**** - src/Symmetry_registration_oldCartGrid3D.c:
** Fixed! (both thorns)

HW: OK. Thanks!

‌

**** - src/BSSN_to_ADM.c
**** - [CHECK:] trafo is $\gamma_{xx} = \frac{1}{cf^2} (\tilde{\gamma}{xx}+1 )$ (why "+1"?) and same for yy and zz components
**** - src/BoundaryCondition_NewRad.c
**** - [CORRECTION:] notation for NewRad is
**** NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower)
**** where var0 is the background value of the function. For $\alpha$, and diagonal components of $\tilde{\gamma}{ij}$ this is 1.
**** In "src/BoundaryCondition_NewRad.c" they are registered as $0.0$ and should be corrected.
** The two comments are actually related; we evolve h{ij} = gamma_{ij} - delta_{ij}, where delta_{ij} is the Kronecker delta.
** Therefore the value of h_{ij} should indeed be zero at infinity.

HW: OK. Thanks for the clarification.

‌

** Regarding the lapse having the wrong value at infinity -- nice catch! I have fixed and pushed that (both thorns).

HW: OK. Thanks!

‌

**** - [COSMETICS:] src/ADM_to_BSSN.c and other files: typo in comment, l.58 of source code: "Step 2 of 1"
** This has been fixed and pushed (both thorns)!

HW: OK

‌

** Per Roland's suggestion, I also added the NRPy+ git hash used to generate the thorns in each thorn's README file.

HW: OK. Great.

[Roland Haas]

@helvi witek @Zach Etienne is there anything left to do here? BaikalETK is part of the ET as of ET_2020_05. If this is all done, if reviewer and reviewee are in agreement, please close this ticket as resolved.

[Zach Etienne]

The BaikalETK thorns Baikal & BaikalVacuum were reviewed, revised, approved, and merged into the ET_2020_05 of the Einstein Toolkit.

[Roland Haas]

For the record, there is currently no officially requested or suggested citation for Baikal (on http://einsteintoolkit.org/citation.html). @Zach Etienne suggested using the original NRPy+ paper https://arxiv.org/abs/1712.07658 until one is added.