Add ADM Mass / Spin and Extrapolation to Infinity

Examples/BinaryBH/params.txt

@@ -180,6 +180,9 @@ modes = 2 0 # l m for spherical harmonics
         4 3
         4 4
 
+extraction_extrapolation_order = 3
+extraction_extrapolation_radii = 0 1 2
+
 # integral_file_prefix = "Weyl4_mode_"
 
 # write_extraction = 0

Examples/KerrBH/DiagnosticVariables.hpp

@@ -15,6 +15,9 @@ enum
     c_Mom2,
     c_Mom3,
 
+    c_Madm,
+    c_Jadm,
+
     NUM_DIAGNOSTIC_VARS
 };
 
@@ -23,7 +26,9 @@ namespace DiagnosticVariables
 static const std::array<std::string, NUM_DIAGNOSTIC_VARS> variable_names = {
     "Ham",
 
-    "Mom1", "Mom2", "Mom3"};
+    "Mom1",  "Mom2", "Mom3",
+
+    "M_adm", "J_adm"};
 }
 
 #endif /* DIAGNOSTICVARIABLES_HPP */

Examples/KerrBH/KerrBHLevel.cpp

@@ -6,7 +6,7 @@
 #include "KerrBHLevel.hpp"
 #include "BoxLoops.hpp"
 #include "CCZ4RHS.hpp"
-#include "ChiTaggingCriterion.hpp"
+#include "ChiExtractionTaggingCriterion.hpp"
 #include "ComputePack.hpp"
 #include "KerrBHLevel.hpp"
 #include "NanCheck.hpp"
@@ -20,6 +20,9 @@
 #include "GammaCalculator.hpp"
 #include "KerrBH.hpp"
 
+#include "ADMQuantities.hpp"
+#include "ADMQuantitiesExtraction.hpp"
+
 void KerrBHLevel::specificAdvance()
 {
     // Enforce the trace free A_ij condition and positive chi and alpha
@@ -38,9 +41,9 @@ void KerrBHLevel::initialData()
     if (m_verbosity)
         pout() << "KerrBHLevel::initialData " << m_level << endl;
 
-    // First set everything to zero then calculate initial data  Get the Kerr
-    // solution in the variables, then calculate the \tilde\Gamma^i numerically
-    // as these are non zero and not calculated in the Kerr ICs
+    // First set everything to zero then calculate initial data
+    // Get the Kerr solution in the variables, then calculate the \tilde\Gamma^i
+    // numerically as these are non zero and not calculated in the Kerr ICs
     BoxLoops::loop(
         make_compute_pack(SetValue(0.), KerrBH(m_p.kerr_params, m_dx)),
         m_state_new, m_state_new, INCLUDE_GHOST_CELLS);
@@ -97,5 +100,39 @@ void KerrBHLevel::preTagCells()
 void KerrBHLevel::computeTaggingCriterion(FArrayBox &tagging_criterion,
                                           const FArrayBox &current_state)
 {
-    BoxLoops::loop(ChiTaggingCriterion(m_dx), current_state, tagging_criterion);
+    BoxLoops::loop(ChiExtractionTaggingCriterion(m_dx, m_level, m_p.max_level,
+                                                 m_p.extraction_params,
+                                                 m_p.activate_extraction),
+                   current_state, tagging_criterion);
 }
+
+void KerrBHLevel::specificPostTimeStep()
+{
+    CH_TIME("KerrBHLevel::specificPostTimeStep");
+    // Do the extraction on the min extraction level
+    if (m_p.activate_extraction == 1)
+    {
+        int min_level = m_p.extraction_params.min_extraction_level();
+        bool calculate_adm = at_level_timestep_multiple(min_level);
+        if (calculate_adm)
+        {
+            // Populate the ADM Mass and Spin values on the grid
+            fillAllGhosts();
+            BoxLoops::loop(ADMQuantities(m_p.extraction_params.center, m_dx,
+                                         c_Madm, c_Jadm),
+                           m_state_new, m_state_diagnostics,
+                           EXCLUDE_GHOST_CELLS);
+
+            if (m_level == min_level)
+            {
+                CH_TIME("ADMExtraction");
+                // Now refresh the interpolator and do the interpolation
+                m_gr_amr.m_interpolator->refresh();
+                ADMQuantitiesExtraction my_extraction(
+                    m_p.extraction_params, m_dt, m_time, m_restart_time, c_Madm,
+                    c_Jadm);
+                my_extraction.execute_query(m_gr_amr.m_interpolator);
+            }
+        }
+    }
+}

Examples/KerrBH/KerrBHLevel.hpp

@@ -39,6 +39,9 @@ class KerrBHLevel : public GRAMRLevel
     /// Things to do before tagging cells (i.e. filling ghosts)
     virtual void preTagCells() override;
 
+    // to do post each time step on every level
+    virtual void specificPostTimeStep() override;
+
     virtual void
     computeTaggingCriterion(FArrayBox &tagging_criterion,
                             const FArrayBox &current_state) override;

Examples/KerrBH/Main_KerrBH.cpp

@@ -40,6 +40,16 @@ int runGRChombo(int argc, char *argv[])
     DefaultLevelFactory<KerrBHLevel> kerr_bh_level_fact(gr_amr, sim_params);
     setupAMRObject(gr_amr, kerr_bh_level_fact);
 
+    // Set up interpolator:
+    // call this after amr object setup so grids known
+    // and need it to stay in scope throughout run
+    // Note: 'interpolator' needs to be in scope when gr_amr.run() is called,
+    // otherwise pointer is lost
+    AMRInterpolator<Lagrange<4>> interpolator(
+        gr_amr, sim_params.origin, sim_params.dx, sim_params.boundary_params,
+        sim_params.verbosity);
+    gr_amr.set_interpolator(&interpolator);
+
     using Clock = std::chrono::steady_clock;
     using Minutes = std::chrono::duration<double, std::ratio<60, 1>>;
 

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);
     }
 
     void check_params()
@@ -51,6 +53,7 @@ class SimulationParameters : public SimulationParametersBase
     }
 
     KerrBH::params_t kerr_params;
+    bool activate_extraction;
 };
 
 #endif /* SIMULATIONPARAMETERS_HPP_ */

Examples/KerrBH/params.txt

@@ -95,6 +95,7 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
                          0 1 2 3          #Theta and Gamma
                          0 1 2 3 1 2 3    #lapse shift and B
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
+                         0 0              #ADM M and J
 
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
@@ -144,4 +145,26 @@ sigma = 0.3
 # min_lapse = 1.e-4
 
 #################################################
+# Extraction parameters
+
+# extraction_center = 256 256 256 # defaults to center
+activate_extraction = 1
+num_extraction_radii = 3
+extraction_radii = 30. 40. 50.
+extraction_levels = 2 1 0
+num_points_phi = 50
+num_points_theta = 51
+# num_modes = 3
+# modes = 2 0 # l m for spherical harmonics
+#        2 1
+#        2 2
+
+# integral_file_prefix = "Weyl4_mode_"
+
+# write_extraction = 0
+# extraction_subpath = "data/extraction" # directory for 'write_extraction = 1'
+# extraction_file_prefix = "Weyl4_extraction_"
+
+#################################################
+
 

Examples/KerrBH/params_cheap.txt

@@ -95,6 +95,7 @@ vars_parity            = 0 0 4 6 0 5 0    #chi and hij
                          0 1 2 3          #Theta and Gamma
                          0 1 2 3 1 2 3    #lapse shift and B
 vars_parity_diagnostic = 0 1 2 3          #Ham and Mom
+                         0 0              #ADM M and J
 
 # if sommerfeld boundaries selected, must select
 # non zero asymptotic values
@@ -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
 
 # Lapse evolution
 lapse_advec_coeff = 1.0
@@ -144,3 +145,29 @@ sigma = 0.3
 # min_lapse = 1.e-4
 
 #################################################
+# Extraction parameters
+
+# extraction_center = 256 256 256 # defaults to center
+activate_extraction = 1
+num_extraction_radii = 3
+extraction_radii = 30. 40. 50.
+extraction_levels = 1 1 0
+num_points_phi = 50
+num_points_theta = 51
+# num_modes = 3
+# modes = 2 0 # l m for spherical harmonics
+#        2 1
+#        2 2
+
+extraction_extrapolation_order = 3
+extraction_extrapolation_radii = 0 1 2
+
+# integral_file_prefix = "Weyl4_mode_"
+
+# write_extraction = 0
+# extraction_subpath = "data/extraction" # directory for 'write_extraction = 1'
+# extraction_file_prefix = "Weyl4_extraction_"
+
+#################################################
+
+

Source/AMRInterpolator/SurfaceExtraction.hpp

@@ -54,19 +54,25 @@ template <class SurfaceGeometry> class SurfaceExtraction
         std::string data_path, integral_file_prefix;
         std::string extraction_path, extraction_file_prefix;
 
+        std::vector<int>
+            radii_idxs_for_extrapolation; // 2 for 2nd order, 3 for 3rd order
+
         int min_extraction_level()
         {
             return *(std::min_element(extraction_levels.begin(),
                                       extraction_levels.end()));
         }
+
+        //! deletes invalid or repeated
+        void validate_extrapolation_radii();
     };
 
     using vars_t = std::tuple<int, VariableType, Derivative>;
 
   protected:
     const SurfaceGeometry m_geom; //!< the geometry class which knows about
                                   //!< the particular surface
-    const params_t m_params;
+    params_t m_params;
     std::vector<std::tuple<int, VariableType, Derivative>>
         m_vars; //!< the vector of pairs of
     //!< variables and derivatives to extract
@@ -189,8 +195,12 @@ template <class SurfaceGeometry> class SurfaceExtraction
     //! convenience caller for write_integrals in the case of just integral per
     //! surface
     void write_integral(const std::string &a_filename,
-                        const std::vector<double> a_integrals,
-                        const std::string a_label = "") const;
+                        const std::vector<double> &a_integrals,
+                        const std::string &a_label = "") const;
+
+    // returns empty vector if no extrapolation is done
+    std::vector<double> richardson_extrapolation(
+        const std::vector<std::vector<double>> &integrals) const;
 };
 
 #include "SurfaceExtraction.impl.hpp"