Combine Shape map calls

[nikwit]

Currently, the shape map computes the coordinates, frame velocity and jacobians separately which duplicates a lot of work such as computing the spherical harmonic expansion.

This PR adds a combined call to these three quantities. I found that this speeds up simulations between 16 % and 27% (when combined with #6223)

depends on #6223

[knelli2]

Thank you for investigating this!

[knelli2]

Instead of checking if this is a shape map explicitly (and thus adding an strange dependency for the general coord map on the time dependent maps), add a meta function that checks if the Map object has a function named coords_frame_velocity_jacobian that can be called with the given arguments.

Then this easily allows us to edit the other time dependent maps if we want to to add this feature to those as well

[nilsdeppe]

Note: I think we could get away with a forward declaration of the map to avoid the coupling. Just providing another option.

[nilsdeppe]

LGTM, no need for my approval to merge.

[nilsdeppe]

3, 3?

[nilsdeppe]

Note: I think we could get away with a forward declaration of the map to avoid the coupling. Just providing another option.

[nikwit]

@knelli2 I have factored out large chunks of the jacobian, I do not think this is worth doing for the operator() or frame velocity as it is just two statements. I also avoid duplicate calls to the transition_function in the combined call which I had overlooked before.

[knelli2]

Looks good! You can squash this last change in with everything else

[knelli2]

Also fix the compilation errors from gcc-9 and the clang-tidy issues. The one about using l as a variable you can ignore or add it to our .clang-tidy file because it's pretty silly

[nikwit]

Thanks for your reviews. I changed everything but the long link domain::CoordinateMaps::ShapeMapTransitionFunctions::ShapeMapTransitionFunction::original_radius_over_radius in the documentation Shape.hpp:162 is over 80 characters and not resolved if I break it up over two lines. Does someone know how this can be resolved?

[knelli2]

You can squash directly. As for the long link, our commit hooks check for linked that contain either \link or \endlink and ignore those, so just make sure one of those is on the same line as the link.

[nikwit]

all done and squashed in, thank you for the review!

[knelli2]

Thanks for optimizing this!

[knelli2]

This test failure is concerning me since I (and @nilsvu) have never seen it before and I want to make sure it's not related to the shape map changes. The test failure happened twice in a row with the same error

[nikwit]

oh, that is weird. So it only fails with clang-13 in Release mode. I don't know the executable at all. Would this executable call the coords_framevel_jac function of CoordinateMap? If so, do you know which compile time and runtime path it would take in that function? @nilsvu

[nikwit]

Looks like reducing the LMax reduced the weird test failure of clang-13 in the Xcts executable @knelli2

[nikwit]

@knelli2 @nilsvu @nilsdeppe I looked into the test failure a bit more. It seems that the test failure is unrelated to combining the shape map calls. I can "fix" the failure by simply adding or removing some temporaries in the shape map that change the evaluation order of operations (all additions and multiplications). One change that causes the test to fail is e.g. here, where I just introduce a temporary that should definitely be safe: nikwit@6ecb730

Also, simply inlining the jacobian_helper function causes the test to work.

I am therefore convinced that this PR simply changes some floating point precision that propagates through the elliptic solver and changes the residuals that are computed (which are still convergent). IMO, we should simply raise the tolerance of that test.

[knelli2]

@nilsvu I think it is up to you what to do here since this only affects the elliptic solver.

Also @nikwit this has a conflict now and needs a rebase

[nilsvu]

I ran this branch on my machine and the test passes. It looks like specifically clang-13 Release fails, and specifically the first GMRES iteration after AMR iteration is different. The shape map gets evaluated on the new refined grid points, so at least it's plausible that an issue in the shape map can cause this failure. I worries me enough that I wouldn't just reduce the Lmax to "fix" it without checking that BBH evolutions aren't affected. Suggestions:

• How about adding a regression test to GH/KerrSchild.yaml like the one in Xcts/KerrSchild.yaml, to check that the first couple of time steps remain identical? This should use the same Kerr spin and shape map Lmax that caused issues in XCTS.
• ylm::Spherepack does some fancy caching in a static thread_local memory pool, maybe that's causing issues somehow? Just speculating. I guess you could try making it a normal member variable by removing the static thread_local just to see if that affects the issue.
• You could also run a BBH on clang-13 Release to make sure the results are the same before and after this PR (to numerical roundoff error, not just truncation error).

It's a very mysterious issue, and I just want to ensure that it doesn't affect our BBH simulations in mysterious ways 😅

[nikwit]

@nilsvu I think my analysis above shows that the issue is related simply to the order in which expressions are evaluated in the shape map (nikwit@6ecb730 fixes it!). The test appears to blow up a floating point difference to a difference of 1e-9 somehow. I do not see how this can be related to ylm::Spherepack.

My suspicion is that the test is not a fair test. Can you describe where the residuals that are compared in the test come from? It looks like the test uses 4 grid points per element and no h-refinement. How can this lead to a residual of 1e-11?

[nilsvu]

Can you describe where the residuals that are compared in the test come from? It looks like the test uses 4 grid points per element and no h-refinement. How can this lead to a residual of 1e-11?

The refinement of the domain and the DG scheme defines a discretized problem $A(u_q)_p = b_p$ over the grid points (it's a matrix equation if the problem is linear). The residual is the $L_2$ norm $|| A(u_q)_p - b_p||_2$ over the grid points. This discretized problem can be solved for the values $u_p$ down to machine precision (or to 1e-11 in this case), independent of resolution. If the problem is linear it's just a matrix inversion, which can be done to machine precision. This accuracy is only the numerical solution to the discretized problem, it doesn't mean the truncation error is so small as well. The truncation error is limited by how well the discretized problem represents the true problem, and is ~1e-3 in this case (see the "Discretization error" check in KerrSchild.yaml, which compares the solution $u_p$ to the analytic solution at every grid point).

Now I don't see how the numerical algorithm that solves the discretized problem should blow up roundoff error in the shape map to truncation error level differences, though I might be missing something of course.

[nikwit]

alright, I see.

At this point I would argue that since changing the order in which DataVectors are evaluated fixes the test, this must mean that any issue that may or may not be in the code is unrelated to this PR. If anything, the failing test is a symptom of an issue already on develop (though I do not believe there is one).

[nilsvu]

I have now seen two other cases of this test failing on CI, so I agree that it is not directly related to the changes in this PR and I'm happy with this PR getting merged. I still have no idea what causes the issue though :(

[knelli2]

@nikwit Just needs another rebase to fix the conflict and then I'll merge this. @nilsvu Thanks for looking into this!

[nikwit]

rebased. Thank you all for your due diligence!