Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add volume integrands for XCTS asymptotic quantities. #6276

Merged
merged 1 commit into from
Sep 23, 2024

Conversation

iago-mendes
Copy link
Member

Proposed changes

This PR adds volume integrands for the calculation of ADM mass, ADM linear momentum, and center of mass in the XCTS system.

The primary goal of using infinite volume integrals instead of infinite surface integrals is to improve numerical accuracy. We haven't seen significant accuracy improvements yet, and we have found some analytic test cases in which the volume integrals do not give the expected results (which could be a bug in the method or an assumption being broken somewhere in these cases). Then, we'll keep using the only surface integrals in the ID parameter control for now.

Upgrade instructions

Code review checklist

  • The code is documented and the documentation renders correctly. Run
    make doc to generate the documentation locally into BUILD_DIR/docs/html.
    Then open index.html.
  • The code follows the stylistic and code quality guidelines listed in the
    code review guide.
  • The PR lists upgrade instructions and is labeled bugfix or
    new feature if appropriate.

Further comments

Copy link
Member

@nilsvu nilsvu left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't reviewed the tests yet. Please push changes as fixup commits.

* \bar\gamma^{jk} \partial_i \bar\Gamma^i_{jk}
* + \bar\Gamma^i_{jk} \partial_i \bar\gamma^{jk}
* - \bar\gamma^{ij} \partial_i \bar\Gamma_j
* - \bar\Gamma_j \partial_i \bar\gamma^{ij}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is missing 2 terms: $+\gamma^{jk} \Gamma_l \Gamma^l_{jk} - \gamma^{lj} \Gamma_l \Gamma_j$ where $\Gamma_l = \Gamma^k_{kl}$.
I don't think these cancel since $\gamma^{jk}\Gamma^l_{jk} \neq \gamma^{lj}\Gamma^k_{jk}$ (or do they somehow?)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand it right, these extra terms would only show up if we take the covariant derivative $\bar D_i$ when applying Gauss' law instead of a partial derivative $\partial_i$ (this is what I'm currently doing). I remember we briefly talked about these two options (referencing Saul's recent DG paper), and we decided that $\partial_i$ was the way to go.

That said, after redoing the calculations, I found an extra term involving the conformal factor first derivative:
$\partial_i ( - 8 \bar\gamma^{ij} \partial_j \psi ) = - 8 \bar\gamma^{ij} \partial_i \partial_j \psi - 8 \partial_i \bar\gamma^{ij} \partial_j \psi = - 8 \bar D^2 \psi - 8 \partial_i \bar\gamma^{ij} \partial_j \psi $.
Again, this is only if we use $\partial_i$ in Gauss' Law. Looking at Baumgarte-Shapiro's equation (3.146), they seem to use $\bar D_i$ in Gauss' Law.

I'm not sure what approach we should take here, but both of them require updates to these integrands.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need the covariant Laplacian of $\psi$ to replace it with the Hamiltonian constraint, so I would go with the covariant derivatives and the Christoffel terms I wrote above. As you say, either way you have to account for the covariant derivatives with Christoffel terms.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, that makes sense. This raises another point. If we're using $\bar D_i$ when applying Gauss' law, then can we still assume a flat space in this transformation?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For Gauss' law you transform $\oint n_i \sqrt{\gamma} F^i = \int \partial_i (\sqrt{\gamma} F^i) = \int \sqrt{\gamma} \nabla_i F^i$ (see Eqs. (7-10) in https://arxiv.org/abs/physics/9802027v1). So you get the covariant derivative from moving the metric determinant through the partial derivative.
So, you integrate with the curved area/volume element, get a covariant derivative from Gauss' law, and expand it into Christoffel terms where needed.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is helpful, thanks! The test that I had in Test_AdmMass was using flat area/volume elements and started failing after these changes. Once I updated this test to use curved area/volume elements, I started getting the correct values again. So, things seem consistent.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok great! Both approaches work, we just have to be consistent. Please note in the docs for the integrand function if it needs to be integrated with the flat or conformal area/volume element.

src/PointwiseFunctions/Xcts/AdmMass.hpp Outdated Show resolved Hide resolved
src/PointwiseFunctions/Xcts/AdmMass.hpp Outdated Show resolved Hide resolved
src/PointwiseFunctions/Xcts/AdmMass.cpp Outdated Show resolved Hide resolved
src/PointwiseFunctions/Xcts/CenterOfMass.hpp Outdated Show resolved Hide resolved
src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp Outdated Show resolved Hide resolved
src/PointwiseFunctions/Xcts/CenterOfMass.cpp Outdated Show resolved Hide resolved
src/PointwiseFunctions/Xcts/CenterOfMass.hpp Outdated Show resolved Hide resolved
Copy link
Member

@nilsvu nilsvu left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can squash these two commits (including this change)

src/PointwiseFunctions/Xcts/CenterOfMass.cpp Outdated Show resolved Hide resolved
@nilsvu
Copy link
Member

nilsvu commented Sep 16, 2024

Looks good, you can squash the fixup.

@iago-mendes iago-mendes force-pushed the volume-integrands branch 2 times, most recently from 0b09e1f to 016f6ed Compare September 16, 2024 21:52
@nilsvu
Copy link
Member

nilsvu commented Sep 16, 2024

New fixups look good too, you can squash

Copy link
Member

@nilsvu nilsvu left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good. You can squash everything, including this.

@@ -195,6 +196,15 @@ void test_local_adm_integrals(const double& distance,
deriv_conformal_metric, inv_conformal_metric);
const auto conformal_christoffel_contracted = tenex::evaluate<ti::i>(
conformal_christoffel_second_kind(ti::J, ti::i, ti::j));
const auto deriv_conformal_christoffel_second_kind = partial_derivative(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At the moment you don't use this yet, right? Same for Ricci scalar/tensor and longitudinal shift below. Maybe hold off with adding these until you use them in this test.

}
}
}

// Check result
auto custom_approx = Approx::custom().epsilon(10. / distance).scale(1.0);
auto custom_approx = Approx::custom().epsilon(TOLERANCE).scale(1.0);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you keep this at O(1/distance)?

Copy link
Member Author

@iago-mendes iago-mendes Sep 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't work for this test because the CoM integral "diverges" as we increase R. I believe this happens due to round-off errors when canceling larger and larger terms. I describe a bit of this at the top of this file.


// Set up domain
const size_t h_refinement = 1;
const size_t p_refinement = 11;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is very high refinement. My understanding is that the accuracy of the integral is currently not limited by resolution, so can you lower this?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was the lowest resolution that I was able to find in order to keep the tolerance at $10^{-3}$. Figure 3.2.1 in these convergence plots shows that the volume integral residual gets slightly better for higher resolution.

I think there must be some change that we can make in the calculation of the CoM integrals to get around these round-off errors, but this is probably something for a future PR.

nilsvu
nilsvu previously approved these changes Sep 20, 2024
@nilsvu
Copy link
Member

nilsvu commented Sep 20, 2024

ADM mass test times out, probably just relax the timeout

@iago-mendes
Copy link
Member Author

ADM mass test times out, probably just relax the timeout

How do I relax the timeout?

@nilsvu
Copy link
Member

nilsvu commented Sep 21, 2024

Add // [[TimeOut, 10]] above the SPECTRE_TEST_CASE

@iago-mendes
Copy link
Member Author

It seems like it's still timing out for some compilers. Should I relax it even further?

@nilsvu
Copy link
Member

nilsvu commented Sep 23, 2024

Yes just set a reasonable timeout, and/or speed up the test

@iago-mendes
Copy link
Member Author

Yes just set a reasonable timeout, and/or speed up the test

Done. The ADM mass test now succeeds after <20s for the compilers that were failing before (gcc-9, gcc-11). Now the fails are from PartialDerivs, so it's unrelated to my changes.

@nilsvu nilsvu merged commit 20b3671 into sxs-collaboration:develop Sep 23, 2024
21 of 23 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants