Add ADM Mass / Spin and Extrapolation to Infinity #121

TaigoFr · 2020-05-06T19:24:25Z

Addition of a class ADMMass that calculates the ADM mass and spin of a spacetime, together with a class ADMMassExtraction that extracts the corresponding variables.

This partially solves #90 , and will have to be changed once diagnostics are added. #90 also includes the Kretschmann scalar and such variables, but that's in my opinion a feature to be added separately.

This PR also includes Extrapolation to infinity of first and second order. First order is in Alcubierre. Second order I wrote some notes about it, as I didn't find it elsewhere and made the computation myself. Second order (requires 3 radii to cancel leading terms of order 1/r AND 1/r^2) greatly improves the result, so it must be correct.

I added this extraction to the Kerr example, as it's a neat place for it to be. For a Kerr BH with Mass=1.0 and Spin=0.3 I managed to get in my laptop a quick extraction that with 2nd order extraction to infinity results in M=1.0000609955 and a=0.29994697676, which looks pretty good to me. I used N=64, L=128, extraction radii at r = (30. 40. 50.) at levels (1 1 0). The simpson method of @mirenradia's SurfaceExtraction also greatly helps this. Other methods (except for boole rule) or first order extraction get slightly worse results, but still converge.

ADM extraction has some ambiguity in terms of what factors of chi should be added. In principle it doesn't matter. I tested that indeed it doesn't matter much, once extrapolation is done. So I added them trying to follow standard calculations, but it is open for review.

KAClough · 2020-05-10T12:40:33Z

This is also very nice! Especially nice to have the extrapolation of the quantities to infinity output automatically.

A few minor comments which I will put in the code.

Examples/KerrBH/KerrBHLevel.cpp

KAClough · 2020-05-10T12:52:08Z

Examples/KerrBH/KerrBHLevel.cpp

+    if (m_p.activate_extraction == 1)
+    {
+        // Populate the Weyl Scalar values on the grid
+        fillAllGhosts();


Comment should say Madm not Weyl.

I have always thought that filling the Weyl/Madm data on every level post timestep is quite inefficient. This is especially true here when we only actually use it every coarse or second to coarse timestep.

Perhaps we could do something like
if(m_time % extraction_dt == 0) {calculate diagnostics}
I think the modulo thing is usually non trivial as you have to check to numerical precision. Having a small function in GRAMRLevel that returns true if m_time is a multiple of dt for a given level would be quite useful generally, and could then be used here.

This can be filed as an issue if you don't want to do it now :-)

A quick alternative is to just move it to inside the following if(m_level == m_p.extraction_params.min_extraction_level()) no? Like that it is only calculated on the moments we want output. I'm not sure this is the same as what you wanted, let me know.
EDIT: ignore this silly comment. This does not work.

Careful @TaigoFr, you want the BoxLoops::loop called on all levels, not just the extraction level.

I realized that and edited my comment. I solved it in the last commit. It's a solution that is stable under numerical perturbations.

Nice fix :-)
Can we also add it to BinaryBHLevel::specificPostTimeStep()? It probably doesn't save us much, but it annoys me that we calculate Weyl on every level timestep.

Examples/KerrBH/params_cheap.txt

Source/CCZ4/ADMMass.hpp

Source/utils/ADMMassExtraction.hpp

KAClough · 2020-05-14T09:58:21Z

Subject to the small change I suggested for amending the BBH example I'm happy with this, @mirenradia are you?

mirenradia · 2020-06-17T16:24:21Z

I will make a full review at some point but my main comment at the moment is that I think the code for ADM angular momentum should be separated out from the code for the ADM mass (and their respective extraction classes) and ChomboParameters::variable_names_to_enum shouldn't be used to check for this (especially in a compute class where performance is important).

mirenradia · 2020-07-05T11:28:22Z

Now that #119 is merged, we should make these diagnostic variables.

KAClough · 2020-08-04T16:29:20Z

I am wondering pedantically about the "ADMMass" name. It also includes the angular momentum, so I think that as Miren suggests this could be separated out, or it could become ADMQuantities or ADMMetrics or some other name which makes it sound more general. I would also like to add the ADM Momentum in some specified direction (doesn't have to be now, but maybe this influences how we solve this issue).

Apart from this (and potential merge conflicts with the other PRs) I think it looks great.

TaigoFr · 2021-01-12T19:13:59Z

I changed the name from ADMMass to ADMQuantities. I think it's best to not separate the mass and momentum calculations, as they involve the common calculation of the normal dS. Since now the mass and momentum diagnostic components are specified in the constructor and defaulted to "don't calculate" (a_c_Madm = -1, int a_c_Jadm = -1), it's hopefully good enough.
Updated to the most recent code as well.

mirenradia · 2021-01-16T17:13:56Z

Can it be made possible for the user to select (via a parameter) which radii to use for the extrapolation?

Add Richardson Extrapolation to Infinity

In BinaryBH Example, calculate Weyl only on the minimum extraction level full cycles Allow the user to extract only the ADM mass without adding the spin as a variable

…r storage (more like the style of Miren's NewConstraints class)

# Conflicts: # Examples/KerrBH/params.txt # Examples/KerrBH/params_cheap.txt # Source/AMRInterpolator/SurfaceExtraction.hpp # Source/GRChomboCore/SimulationParametersBase.hpp

…ation)

mirenradia

Thanks for adding this code and sorry for the long delay in reviewing.

Would it be possible to add the Richardson extrapolation in a new write_xxx function e.g. SurfaceExtraction::write_integrals_and_extrapolate rather than hacking it into the existing function?

mirenradia · 2021-09-09T15:51:54Z

Examples/KerrBH/KerrBHLevel.cpp

+{
+    CH_TIME("KerrBHLevel::specificPostTimeStep");
+    // Do the extraction on the min extraction level
+    if (m_p.activate_extraction == 1)


Pretty sure you corrected this elsewhere.

Suggested change

if (m_p.activate_extraction == 1)

if (m_p.activate_extraction)

mirenradia · 2021-09-09T15:55:17Z

Examples/KerrBH/SimulationParameters.hpp

@@ -29,6 +29,8 @@ class SimulationParameters : public SimulationParametersBase
        pp.load("kerr_mass", kerr_params.mass);
        pp.load("kerr_spin", kerr_params.spin);
        pp.load("kerr_center", kerr_params.center, center);
+
+        pp.load("activate_extraction", activate_extraction, false);


Might as well use the one in SimulationParametersBase

mirenradia · 2021-09-09T15:55:27Z

Examples/KerrBH/SimulationParameters.hpp

@@ -51,6 +53,7 @@ class SimulationParameters : public SimulationParametersBase
    }

    KerrBH::params_t kerr_params;
+    bool activate_extraction;


mirenradia · 2021-09-09T15:56:06Z

Examples/KerrBH/params_cheap.txt

@@ -118,7 +119,7 @@ max_steps = 4
 # Spatial derivative order (only affects CCZ4 RHS)
 max_spatial_derivative_order = 4 # can be 4 or 6

-nan_check = 1
+# nan_check = 1


mirenradia · 2021-09-09T15:57:34Z

Source/AMRInterpolator/SurfaceExtraction.hpp

        int min_extraction_level()
        {
            return *(std::min_element(extraction_levels.begin(),
                                      extraction_levels.end()));
        }
+
+        //! deletes invalid or repeated
+        void validate_extrapolation_radii();


Can we do this via parameter checking instead? This will make it more transparent to the user.

mirenradia · 2021-09-09T16:12:03Z

Examples/BinaryBH/params.txt

+extraction_extrapolation_order = 3
+extraction_extrapolation_radii = 0 1 2


Can't extrapolate Weyl4 as time is not retarded.

mirenradia · 2021-09-09T16:17:08Z

Examples/KerrBH/params.txt

+#        2 2
+
+# integral_file_prefix = "Weyl4_mode_"
+


Can you add the extrapolation parameters here (even if just commented out)?

mirenradia · 2021-09-09T16:19:11Z

Source/GRChomboCore/SimulationParametersBase.hpp

+                pp.load("extraction_extrapolation_radii",
+                        extraction_params.radii_idxs_for_extrapolation,
+                        extrapolation_order);


Can we provide the actual radii rather than the radii indices? This can then be checked with parameter checking so that you're not requesting extrapolation from a radii you're not actually extracting on.

mirenradia · 2021-09-10T11:33:35Z

Source/CCZ4/ADMQuantities.hpp

+                Madm += dS_L[i] / (16. * M_PI * m_G_Newton) /
+                        (vars.chi * sqrt(vars.chi)) * h_UU[j][k] * h_UU[i][l] *
+                        (vars.chi * (d1.h[l][k][j] - d1.h[j][k][l]) -
+                         (vars.h[l][k] * d1.chi[j] - vars.h[j][k] * d1.chi[l]));


Maybe add a reference for these formulae as a comment e.g. Baumgarte and Shapiro p85 Eq. (3.128)?

mirenradia · 2021-09-10T11:35:44Z

Source/CCZ4/ADMQuantities.hpp

+                Jadm += -dS_L[i] / (8. * M_PI * m_G_Newton) *
+                        epsilon[m_dir][j][k] * x[j] * vars.K *
+                        TensorAlgebra::delta(i, k);
+
+                FOR(l, m)
+                {
+                    Jadm += dS_L[i] / (8. * M_PI * m_G_Newton) *
+                            epsilon[m_dir][j][k] * x[j] * h_UU[i][l] *
+                            h_UU[k][m] * vars.chi *
+                            (vars.A[l][m] + vars.K * vars.h[l][m] / 3.);


Add a reference e.g. Baumgarte and Shapiro p95 Eq. (3.191).

TaigoFr requested review from KAClough and mirenradia May 6, 2020 19:24

mirenradia added the feature New feature label May 7, 2020

KAClough reviewed May 10, 2020

View reviewed changes

KAClough marked this pull request as draft October 23, 2020 13:27

KAClough mentioned this pull request Jan 5, 2021

Add class to call specificPostTimeStep at t=0 #142

Merged

TaigoFr force-pushed the feature/adm_mass branch 3 times, most recently from f468bfa to 4a60c40 Compare January 19, 2021 17:27

Base automatically changed from master to main January 21, 2021 16:42

TaigoFr added 7 commits March 17, 2021 17:42

Add ADM Mass and Angular Momentum Extraction

e733757

Add Richardson Extrapolation to Infinity

Resolve Katy's comments for ADM Mass PR

6e165f7

Extraction improvements for PR

23e2bfb

In BinaryBH Example, calculate Weyl only on the minimum extraction level full cycles Allow the user to extract only the ADM mass without adding the spin as a variable

Add arguments for user to specify which DiagnosticVariables to use fo…

229908c

…r storage (more like the style of Miren's NewConstraints class)

Change ADMMass name to ADMQuantities

b77e7dd

Add extrapolation radii as parameter

8d2f682

Fixes after rebase

613b4b9

TaigoFr force-pushed the feature/adm_mass branch from 4a60c40 to 613b4b9 Compare March 17, 2021 17:53

TaigoFr marked this pull request as ready for review March 17, 2021 18:21

TaigoFr added 4 commits March 26, 2021 18:33

Merge branch 'main' into feature/adm_mass

6a80629

# Conflicts: # Examples/KerrBH/params.txt # Examples/KerrBH/params_cheap.txt # Source/AMRInterpolator/SurfaceExtraction.hpp # Source/GRChomboCore/SimulationParametersBase.hpp

Add ADMQuantitiesTest (<1% error on M and J on level 0 after extrapol…

03f8455

…ation)

Merge branch 'main' into feature/adm_mass

26a87a5

Merge remote-tracking branch 'origin/main' into feature/adm_mass

ecc83d2

mirenradia requested changes Sep 10, 2021

View reviewed changes

KAClough assigned the-florist Dec 18, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add ADM Mass / Spin and Extrapolation to Infinity #121

Add ADM Mass / Spin and Extrapolation to Infinity #121

TaigoFr commented May 6, 2020

KAClough commented May 10, 2020

KAClough May 10, 2020

TaigoFr May 13, 2020 •

edited

Loading

mirenradia May 13, 2020

TaigoFr May 13, 2020

KAClough May 14, 2020

KAClough commented May 14, 2020

mirenradia commented Jun 17, 2020

mirenradia commented Jul 5, 2020

KAClough commented Aug 4, 2020

TaigoFr commented Jan 12, 2021

mirenradia commented Jan 16, 2021

mirenradia left a comment

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 9, 2021

mirenradia Sep 10, 2021

mirenradia Sep 10, 2021

	if (m_p.activate_extraction == 1)
	if (m_p.activate_extraction)

		extraction_extrapolation_order = 3
		extraction_extrapolation_radii = 0 1 2

Add ADM Mass / Spin and Extrapolation to Infinity #121

Are you sure you want to change the base?

Add ADM Mass / Spin and Extrapolation to Infinity #121

Conversation

TaigoFr commented May 6, 2020

KAClough commented May 10, 2020

Choose a reason for hiding this comment

TaigoFr May 13, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KAClough commented May 14, 2020

mirenradia commented Jun 17, 2020

mirenradia commented Jul 5, 2020

KAClough commented Aug 4, 2020

TaigoFr commented Jan 12, 2021

mirenradia commented Jan 16, 2021

mirenradia left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TaigoFr May 13, 2020 •

edited

Loading