Add NRPyEllipticET to the Einstein Toolkit

[Former user]

NRPyElliptic is a new, hyperbolic relaxation elliptic solver (Assumpcao et al, https://doi.org/10.1103/PhysRevD.105.104037). NRPyEllipticET is the associated Einstein Toolkit thorn, which we propose for inclusion.

In its first application, NRPyElliptic sets up conformally flat, binary black hole (BBH) puncture initial data (ID) on a single numerical domain, similar to TwoPunctures. While it is about 6x slower than TwoPunctures, NRPyElliptic requires only about 0.3% of the runtime for a full BBH simulation in the Einstein Toolkit.

As consumers of numerical relativity ID generally already possess expertise in solving hyperbolic PDEs, they will generally find NRPyElliptic easier to tweak and extend than other elliptic solvers. We are actively working to extend NRPyElliptic to solve other types of ID, and these new ID types will be added to the ET thorn.

[Zach Etienne]

Here is the link to the GitHub repository containing NRPyElliptic & NRPyEllipticET:

https://github.com/assumpcaothiago/NRPyElliptic/

[Roland Haas]

Cheng-Hsin Cheng and Giuseppe Ficarra agreed to review this.

[Roland Haas]

Draft review guidelines are at https://docs.einsteintoolkit.org/et-docs/How_to_Review_a_new_component

[Roland Haas]

@Zach Etienne has there been progress on this?

[Zach Etienne]

@Leonardo Rosa Werneck has an update.

[Leonardo Werneck]

Hi @Roland Haas ! Yes, I’ve already completed most of the major changes we were planning for the thorn. I’m working with Thiago Assumpção (lead developer of NRPyElliptic) to add an extra functionality to the thorn where the damping parameter of the hyperbolic system is computed automatically, yielding optimized relaxations for all grid resolutions (unless the user chooses to override it). This should be complete soon, at which point we will be in touch with the reviewers.

[Roland Haas]

“Soon” is a vague statement. Note that the initial review of harmlessness must finish before 2022-09-01 for NrPyElliptic to be included. The reviewers will also require some time to review which the proposers should provide them with.

Also note that major changes to the code may pose an issue wrt to stating that the code that is being reviewed is actually the one used in a prior publication (ie code quality may actually go down due to newly introduced bugs and loss of correctness compared to the code used in publication).

Usually the ET does not include “bleeding edge” components but components that have proven to themselves already in real world applications (ie used in a scientirfic publication) and it may be desirable to postpone inclusion if that is still to be seen (without passing judgement on the component of course).

[Zach Etienne]

I don’t think this “new” functionality qualifies as “bleeding edge”, as it was hashed out in the publication. It’s basically an empirically derived means to set a more optimal value for the damping parameter.

Also I believe “soon” means “within a few days”.

[Roland Haas]

I would be quite worried about having any sort extra functionality (in particular automatisms) that has not yet been used in publication in the new contribution. There is the worry that one may end up with a new contribution to the ET that claims to have been used in publication X only for a ET users trying to reproduce X and finding out that this is not possible due to the change. Obviously this would be quite embarrassing.

While I trust the NRPyElliptic group to carefully run tests, I wonder why take the risk? Submit the version that was actually used for publication, then add the new functionality in the next release cycle.

Also note that potentially giving the reviewers only a single week for the “does no harm” review is really pushing it. They may be out on a conference, vacation or unavailable for other reasons.

Note that we had called for new contributions in June already (https://docs.einsteintoolkit.org/et-docs/Meeting_agenda#2022-06-02) so starting to adapt the contribution now is somewhat late in my opinion. Before the June call should have been the time for adaptation.

[Roland Haas]

Please be reminded of the deadline to provide the “does no harm” initial review verdict before the deadline of Sept 1st 2022.

[Cheng-Hsin Cheng]

Hi all, I finished checking NRPyEllipticET and I have some questions and comments about the code. To keep it brief I am leaving out the “conformally_flat_BBH” when referring to the file names in the src folder.

---

In NRPyEllipticET folder

• Please add a README in the Cactus format with the author, maintainer, license, and brief description.

• interface.ccl

  • uuGF: Please add a description for it. Although it might be familiar for someone who has used TwoPunctures, a new user might not know what is the purpose of it.

• param.ccl:

  • info_output_freq: I think a better description would be e.g. "Print progress of relaxation time evolution every info_output_freq iterations"
  • N2: allows a "forbidden value" of -1 to make sure it is set explicitly in the parfile, but the default is set to 16 rather than -1. Was the value of -1 to be removed (since N0 and N1 don't have a forbidden value like N2)?
  • initialize_offdiagonal_metric_to_zero: This looks to be unused in the code. Should it be removed from the param.ccl?

---

C source files

• Hyberbolic_Relaxation.c

  • Using CCTK_ERROR (or CCTK_INFO) to error out might be preferable to using printf and calling exit(1).
  • Instead of hard-coding the value for integration_radius, could you add it as a thorn parameter?

• conformally_flat_BBH_driver_bcstruct.c

  • Line 23: It took me a while to understand what the comment meant. I think it might be clearer to rewrite as "in the case that a boundary point is..."
  • Lines 59-61 are commented out. I don't know if set_bcstruct handles grids without the outer boundary, but if it doesn't, should this be left uncommented?

• set_bcstruct.c

  • I really like how at the top there is a brief outline of the routine's (in addition to the comments within the body). It could be useful to add them for other longer routines, e.g. Hyperbolic_Relaxation() or InitializeADMBase()

• find_timestep.c

  • The argument "CFL_FACTOR" isn't used in the function, moreover the variable is already defined in param.ccl, so is this required or can it be removed?

• wavespeed_gf_all_points.c

  • There is duplicate code from find_timestep.c to compute ds_dirn0, ds_dir1, ds_dir2. Is it possible to implement this only once, e.g. write a separate function to compute min(ds_dirn0, ds_dir1, ds_dir2) for a single grid point, and then call it from both find_timestep.c and wavespeed_gf_all_points.c?

• residual_all_points.c

  • The variables "f0_of_xx0" and so on are used in computing the residual, but they are defined with "__attribute__ ((unused))" in the headers in "rfm_files".

• L2_norm_of_gf.c

  • What is the radius "r"? Could it be possible (and would it be a problem) for the extent of r to be under integration_radius, so that the integration is not performed over the whole the radius?

• Files in the folder conformally_flat_BBH/rfm_files/

  • What is the purpose of the empty files with names ending in “read2”?

• Unused files: Should they be left out of the thorn?

  • apply_bcs_curvilinear_radiation.c
  • set_Cparameters-nopointer.h
  • set_Cparameters-SIMD.h
  • xx_to_Cart.c

‌

[Roland Haas]

@Cheng-Hsin Cheng, @Guiseppe Ficarra could you confirm, for the record, that the review concluded?

[Leonardo Werneck]

Hi all. Here are the updates made to NRPyEllipticET based on the reviewers' comments:

• A brief README file has been added;
• uuGF's description has been updated in interface.ccl;
• info_output_freq's comment has been updated in param.ccl;
• N2's “forbidden value” has been removed, as it was unused;
• initialize_offdiagonal_metric_to_zero was indeed an old parameter and has been removed;
• All instances of fprintf(stderr,...); exit(); have been replaced by either CCTK_*ERROR or CCTK_*WARN;
• We have added a new parameter called residual_integration_radius, with a proper description of what it does;
• Some of the algorithms are hard to describe using brief code comments. We recommend looking at the NRPy+ Jupyter notebooks for the complete description of such algorithms;
• Thank you for pointing out that CFL_FACTOR was missing in find_timestep.c; this was actually a bug fix;
• Unfortunately the code in wavespeed_gf_all_points.c depends on dt (the result of find_timestep.c) at every grid point and thus the two cannot be decoupled;
• The reason why the rfm parameters have the unused attribute is similar to the reason why parameters have the attribute, i.e., not every function that uses the rfm includes will use all of the included parameters, which could lead to warnings;
• All instances of __attribute__((unused)) have been replaced with CCTK_ATTRIBUTE_UNUSED;
• I have removed the empty files in conformally_flat_BBH/rfm_files/ manually and am trying to figure out a way of automating this;
• apply_bcs_curvilinear_radiation.c and set_Cparameters-SIMD.h have been removed for now;
• xx_to_Cart.c and set_Cparameters-nopointer.h have not been removed. Although currently unused by the thorn, they provide useful functionality that new users might want to use.

A general comment: the Jupyter notebook that generates NRPyEllipticET was very outdated. I am currently pushing on getting it up to speed, but at the moment we do not have a notebook that regenerates the latest version of the thorn.

[Cheng-Hsin Cheng]

@Roland Haas yes, I would like to confirm that we concluded the review.

[Leonardo Werneck]

@Roland Haas I have added a test directory to NRPyEllipticET.

[Roland Haas]

NRPyEllipticET has been included in the ET as of 2022_11 Kowaleski: https://www.einsteintoolkit.org/about/releases/ET_2022_11_announcement.html