Add Lean_public and Proca thorns to the ET

[Miguel Zilhão]

As discussed during the European ET meeting in Mallorca and during the 2018-01-22 ET telecon, we have a couple of Cactus thorns that we'd like to make available to the general ET community. These include: an updated evolution code that we used to evolve Proca fields (https://arxiv.org/abs/1505.00797), the corresponding analysis and initial data thorns, as well as a general metric evolution thorn (based on Uli Sperhake's Lean). We have been cleaning up the codes, and they are now available in following (public) bitbucket repositories:

https://bitbucket.org/canuda/lean_public https://bitbucket.org/canuda/proca

Included are some testsuites as well as some basic README files. We are currently working on a paper that would also serve as documentation. We would be very grateful if these could be considered for inclusion in the 2018_08 release.

Sincerely,

Helvi Witek Miguel Zilhão

Keyword: Lean_public

[anonymous]

I am reviewing the "LeanBSSNMoL" code for submission to the Einstein Toolkit. Executive summary: The code is beautifully written, but the thorn lacks test cases, which are necessary.

Downloading and building the code worked fine.

The code is written in Fortran 90. It has a clear structure and is very readable. It consists of essentially one large loop per function, which is efficient and compiles quickly, as opposed to using Fortran 90's array syntax that unfortunately leads to bad slow-downs with some compilers.

The code has the usual idiosyncracies where it expects that "CCTK_REAL" is the same as "double precision", and "CCTK_INT" is the same as "integer". These hold usually in practice, and the code is fine as is, although the main point of introducing "CCTK_REAL" and "CCTK_INT" was to allow these types to be distinct.

I am reading the code, and I have the following remarks:

The file "zero_rhs.F90" (and maybe others) do not include "cctk_Arguments.h". How do they even compile? The missing include statement should be added. I wonder whether there is a bug in the Einstein Toolkit in that including "cctk.h" already includes "cctk_Arguments.h"?

I have, for obvious reasons, not thoroughly checked the implementations of finite differences and BSSN equations in the file "bssn_constraints.F90" and "calc_bssn_rhs.F90". I assume that the authors checked their code and am willing to trust them. Personally, I would have defined macros for the lengthy finite differencing expressions.

I ran all test cases successfully. There are no test cases, so this is not a strong statement. Test cases need to be added for a few simple cases, e.g. for some of the Apples with Apples examples (nonlinear gauge wave, Gowdy spacetim) as well as a single, preferably rotating black hole.

I have run the example parameter file on 1 node with 40 cores of the Niagara system of Compute Canada, using 10 MPI processes with 4 OpenMP threads each. This seems to work nicely, albeit progress was so slow that after a few hour only 16 iterations had completed, and horizons were found only once. A quick glance shows that the output looks correct.

For future reference, example parameter files should contain comments near the top describing approximately how much memory and CPU time they use. This makes it much simpler to run them, and prevents people from accidentally running them on their laptop if there is no chance of success. Put differently, the driver should have a notion of the amount of memory available on a node (or laptop), and should refuse to allocate significantly more than that. Simfactory already has the necessary information.

Additional remarks, as suggestions for future improvements to the authors:

The file "evolve_utils.F90" uses "omp parallel do collapse(3)". It might be faster to use only "collapse(2)" here, and instead use "omp simd" for the innermost loop. The same applies to all OpenMP loops.

Some functions have very generic names, such as "remove_trA" or "reset_detmetric". This risks name collisions throughout the executable. I suggest to add a prefix to these functions, even if only a short one such as "lbm_". Alternatively, these functions can be put into a module with a unique name such as "LeanBSSN_evolve_utils", which would have the added benefit of checking the function's signatures. (Scheduled functions cannot be put into modules.)

Other, minor remarks:

The thorn list entry contains the comment "Lean public thorns". This gives the impression as if there was a proper set of Lean thorns, and a few of them were made public. I'd prefer to call them simply "Lean thorns", as these thorns need to be interesting and useful by themselves, and should see some development by themselves. I object to the notion that Einstein Toolkit thorns are merely mirror images of thorns that truly live elsewhere, hidden, and which from time to time receive a data dump as upgrade. Instead, the thorns in the Einstein Toolkit should be the main development versions, and where private copies exist, they should contribute via pull request to the main development versions.

The shift condition mentioned in the README of LeanBSSNMoL should recieve a proper pointer to a particular set of equations (e.g. (S2)?).

I see comments such as "$Header:$" in some files. These are now-unused left-overs from CVS, and can be removed.

Calls to "CCTK_WARN(0, ...)" can be replaced by "CCTK_ERROR(..."), which might be more readable.

[Peter Diener]

Here is my review of the "Proca" suite of thorns that are proposed for inclusion in the Einstein Toolkit. Proca consists of ProcaBase, Proca_simpleID, TwoPunctures_KerrProca and ProcaEvolve and NPScalars_Proca. The code is nicely written, but mainly lack some tests and documentation.

Questions, comments and recommendations regarding Proca overall:

There are two tests in ProcaEvolve of which one uses Lean and the other ML_BSSN while otherwise appear to be the same. Whereas the Ai and Ei variables agree quite well, there seems to be significant differences between Aphi and Zeta. Is this expected? If so, why?

In ProcaEvolve there is another parameter file in the test directory (ML_BSSN_Aphi_mu0.4_c100.1_r06_S0.0_4th.par) that does not have any associated output. If this was meant to be a test, output should be provided. Otherwise it should be moved to the example parameter file directory.

NPScalars_Proca, Proca_simpleID and TwoPunctures_KerrProca have no tests of their own. Proca_simpleID is tested by the tests in ProcaEvolve, but NPScalars_Proca and TwoPunctures_KerrProca should definitely be covered by tests. You may also consider making a standalone test for Proca_simpleID.

None of the thorns in the Proca arrangement have any documentation of their own. As there is a paper describing the formulation, it is okay to refer to the paper, but there should be some description of the variables and parameters used in the code and how they relate to the physics and equations.

It would be nice if all he parameter files for the runs in the paper could be made available as examples of how to use the thorn. Currently there are only two example parameter files available and it's not immediately clear if they have been used for the paper. There could be a comment at the top of the parameterfiles indicating what section and/or figure in the paper they were used for.

Questions, comments and recommendations specifically to ProcaBase:

In parameter.ccl the comment to the value of the parameter mu says: "any non-negative number". However, the range allows for zero and in fact the default value is zero. Either the comment should be changed or the range and default value should be changed.

Are the integer power for the fall off rate for the outgoing boundary conditions for the Proca variables not known? Or can they be different for different physical systems? If they are known constants, should they be parameters at all? Also could the power be different for different components of the E and A vectors. If not, there's no need to make these vector parameters.

I propose to give the file (src/Basegrid.F90) with the symmetry registering routine a more descriptive name.

It would be nice to have a proper (but short) thorn documentation to relate the variables with the paper.

Questions, comments and recommendations specifically to Proca_simpleID:

The Simple Initial Data solution in the paper is a spherical symmetric solution. The parameters for Proca_simpleID on the other hand seem to allow for a superposition of 2 such solutions separated by a non-zero distance. As the parameter for the separation can not be zero (imposed by the parameter range) the initial data thorn will thus never produce initial data that is spherically symmetric around the origin of the cactus grid. Would it make sense for the range of par_b to include zero in order to allow for that. This is something worth explaining in the (missing) documentation.

Questions, comments and recommendations specifically to TwoPunctures_KerrProca:

It looks like the code implements initial data with a Y_11 angular dependence in addition to Y_00 and Y_10 mentioned in the paper. This should be in the documentation.

There's a comment in Equations.c that some of the equations were copied from Mathematica. Could that mathematica notebook be made available in the thorn as well?

In the future it might be worth considering some way of merging this back into the original TwoPunctures, so we don't have multiple thorns implementing the same spectral solver infrastructure.

There should be tests for all the different types of angular profiles for the initial data.

The discussion in the README should be converted into a proper thorn documentation.

Questions, comments and recommendations specifically to ProcaEvolve:

Need a thorn documentation to relate variables to the paper and to explain the use of parameters.

I don't know how much performance is a concern. If only a tiny amount of time is spent in Proca compared to BSSN (or other things) this is obviously not important. However, if the issues with the code, that prevents the compiler from optimizations, makes the code slow enough that the time spent in Proca is comparable to the time spent in BSSN it might be worth it to fundamentally rewrite the code. One of the problems, is the many nested loops over tensor indices. These prevents the compiler from vectorizing most of the RHS code (I checked this by looking at compiler vectorization output from the Intel compiler). As the loop lengths are just 3, the vectorizer deems that vectorization will be ineffecient and hence doesn't do it. On modern CPUs with vector lengths of 4 or 8 (in double precision), efficiently using the vector units is really important for good performance. However, this is not an issue (if there really is a problem with performance) that should prevent inclusion in the ET, but just something to think about for later.