Compactify the Scalar SF - #7469
Conversation
| // 2 r_u - r_u^2/r to compactify the wavezone. | ||
| // Elliptic PDEs will be re-written as a function of | ||
| // inv_r (1/r) instead of r. | ||
| inv_r = (2.0 * r_u - r) / (r_u * r_u); |
There was a problem hiding this comment.
This equation was derived as
x-coordinate r in
There was a problem hiding this comment.
Thanks, this comment was very helpful :D Please do expand the code comment a bit to make this clearer, and/or define sigma = get<0>(x) (see also comment above).
There was a problem hiding this comment.
I added more detail in the explanation. Also defined sigma, directly from the x coordinate.
| * advantages: first, the horizon is placed at a finite radius rather than | ||
| * $r_\star \rightarrow -\infty$; second, the flux component normal to the | ||
| * boundary vanishes at the horizon ($r=r_{\plus}$), which reduces to | ||
| * boundary vanishes at the horizon ($r=r_+$), which reduces to |
There was a problem hiding this comment.
Initially I had r_{\plus}, which appears as "undefined" in the documentation, so this has been changed. (irrelevant to this PR)
| * where $\Delta\phi(r) = \frac{a}{r_+ - r_-} | ||
| * \ln(\frac{r-r_+}{r-r_-})$. |
There was a problem hiding this comment.
Initially I had r_{\plus}, which appears as "undefined" in the documentation, so this has been changed. (irrelevant to this PR)
| * \begin{align} | ||
| * \Psi_m^P &= \frac{r}{2 \pi \sin(\theta)^{|m|}} | ||
| * e^{i m \left( \varphi - \phi\right)} \Phi_m^\mathrm{Wardell} \\ | ||
| * e^{-i m \left( \varphi - \phi\right)} \Phi_m^\mathrm{Wardell} \\ |
There was a problem hiding this comment.
Cillian pointed out that I had typo in this equation. This has been fixed (irrelevant to this PR)
| MModeNumber: *m_mode_number | ||
| HyperboloidalSlicing: True | ||
| Order: 2 | ||
| Upper: None |
There was a problem hiding this comment.
Outer BC was changed as we no longer need to impose anything
| } | ||
| } | ||
|
|
||
| SPECTRE_TEST_CASE( |
There was a problem hiding this comment.
This test checks that, elliptic coefficient (in v or t region) modified by Jacobian returns to the elliptic coefficient in u region. Might be useful if we change elliptic coefficients in the future.
(But I am not sure if this test is necessary)
|
@nilsvu If you could take a look when you have time, that would be great, thanks! |
24e6d91 to
a2938f3
Compare
|
For elliptic coefficients, I re-wrote them as a function of inv_r (=1/r) in the entire region (not just u-region) to simplify the expressions. This is to avoid defining elliptic coefficients differently whether they are in u-region or else where. |
| // 2 r_u - r_u^2/r to compactify the wavezone. | ||
| // Elliptic PDEs will be re-written as a function of | ||
| // inv_r (1/r) instead of r. | ||
| inv_r = (2.0 * r_u - r) / (r_u * r_u); |
There was a problem hiding this comment.
Thanks, this comment was very helpful :D Please do expand the code comment a bit to make this clearer, and/or define sigma = get<0>(x) (see also comment above).
There was a problem hiding this comment.
I believe this code will be wrong in the u region. Is this a problem? It's probably better not to perform BL-to-tortoise transforms where sigma, inv_r etc. here properly, rather than below.
There was a problem hiding this comment.
I removed r_star completely, as it was never used when penetrating_horizon is true (this is why the problem was never raised, even though the definition was wrong)
| get<0>(alpha) = delta / r_sq_plus_a_sq; | ||
| get<1>(alpha) = 1.0 / r_sq_plus_a_sq; | ||
| get<0>(alpha) = (delta_over_r_sq / one_plus_a_sq_inv_r_sq) * dsigma_dr; | ||
| get<1>(alpha) = (square(inv_r) / one_plus_a_sq_inv_r_sq) / dsigma_dr; |
There was a problem hiding this comment.
I think it's better to be a bit more careful here with dsigma_dr -> 0. I think it would be reasonable to have in_u_region branches that take care of cancellations (where needed), rather than force everything into the 1/r form.
There was a problem hiding this comment.
Yes I agree. thanks for catching this issue. When outside of u-region, I kept the original expression (without 1/r reformulation) and inside u-region, division by dsigma_dr is now taken care of, by some cancellation so that we no longer have 1/0 in u-region.
78db2fa to
d680ba8
Compare
| @@ -104,8 +104,6 @@ CircularOrbit::variables( | |||
| // NOLINTNEXTLINE | |||
| r.set_data_ref(const_cast<DataVector*>(&r_star_or_r)); | |||
There was a problem hiding this comment.
This line is incorrect in the u-region where get<0>(x) is actually sigma. This means that many of the computations below are incorrect (like delta etc). I think you're careful about only using correct computations in if (in_u_region) branches, but this is still confusing and easy to introduce bugs. Can you define sigma here properly and make sure all computations below are correct, or (if this makes the code too complicated or leads to unnecessary computations) add comments to clarify exactly which quantities are computed correctly in which region?
There was a problem hiding this comment.
You are right. Though they were never used in u-region, it was confusing and misleading. I re-defined r and sigma when penetrating_horizon is true, so that computations below (delta, etc) are always true.
Proposed changes
CircularOrbit.hppUpgrade instructions
Code review checklist
make docto generate the documentation locally intoBUILD_DIR/docs/html.Then open
index.html.code review guide.
bugfixornew featureif appropriate."Co-Authored-By: " as the last line of the
commit, e.g. "Co-Authored-By: Claude Sonnet 4.6 noreply@anthropic.com",
"Co-Authored-by: Codex noreply@openai.com", or
"Co-Authored-By: GitHub Copilot CLI noreply@microsoft.com".
Further comments