First draft unifying WeylExtraction and particles

The following logic is enforced in this commit:
- Particle interpolator is stored in BHAMR (as well as particles are
  populated once)
- Particle interpolator is also now templated over the number of
  components
- We ask for Redistribute() of particles after regrid
- GRAMR contains a very short function to post-process the MultiFabs of derived
  variables into the correct format for interpolation
- GRAMR contains some virtual function to be overriden by BHAMR in order
  to set up query and populate the particles only once
- Handling for parity of derived vars is implemented. Parity handling
  for both state and derived is to be tested!
- Weyl4 is computed in the post-processing step for now; to be addressed later but
  is needed for WeylExtraction testing

Author: tamaraevst

Examples/BinaryBH/BinaryBHLevel.cpp

@@ -12,11 +12,11 @@
 #include "PunctureTagger.hpp"
 #include "PunctureTracker.hpp"
 // xxxxx #include "SixthOrderDerivatives.hpp"
-#include "CustomExtraction.hpp"
+// #include "CustomExtraction.hpp"
 #include "TraceARemoval.hpp"
 #include "TwoPuncturesInitialData.hpp"
 #include "Weyl4.hpp"
-// #include "WeylExtraction.hpp"
+#include "WeylExtractionParticle.hpp"
 
 BHAMR<BinaryBHLevel::num_punctures> *BinaryBHLevel::get_bhamr_ptr()
 {
@@ -252,6 +252,25 @@ void BinaryBHLevel::tag_cells(amrex::TagBoxArray &a_tag_box_array,
     amrex::Gpu::streamSynchronize();
 }
 
+// This is called by amrex on every level after regridding. After regrid we
+// should redistribute the particles again if BoxArrays or DistributionMappings
+// have changed, for example?
+void BinaryBHLevel::specific_post_regrid(int a_lbase, int a_new_finest)
+{
+    amrex::Print() << "BinaryBHLevel::specific_post_regrid() on level "
+                   << Level() << "\n";
+    if (get_gramr_ptr()->cumTime() > 0.0)
+    {
+        if (auto *bh = get_bhamr_ptr())
+        {
+            if (bh->m_weyl_interpolator)
+            {
+                bh->m_weyl_interpolator->force_redistribute(true);
+            }
+        }
+    }
+}
+
 void BinaryBHLevel::specific_post_init()
 {
     BL_PROFILE("BinaryBHLevel::specific_post_init()");
@@ -295,32 +314,91 @@ void BinaryBHLevel::specific_post_checkpoint(const std::string &a_chk_dir,
 
 void BinaryBHLevel::specificPostTimeStep()
 {
-
-    // Custon extraction
+    // std::cout << "BinaryBHLevel::specificPostTimeStep() on level " << Level()
+    //           << std::endl;
     bool first_step = (parent->levelSteps(0) == 0);
 
-    if (Level() == 1)
+    if (Level() == 0)
     {
-        // set the interpolator
-        ParticleInterpolators<1> interpolator(simParams().boundary_params,
-                                              c_chi);
-        interpolator.set_gramr_ptr(get_gramr_ptr());
-
-        // set up the query and execute it
-        std::array<double, AMREX_SPACEDIM> extraction_origin = {
-            0.0, simParams().L / 2, -4.0};
-
-        double m_time       = get_state_data(State_Type).curTime();
-        double m_dt         = get_gramr_ptr()->dtLevel(Level());
-        double restart_time = get_gramr_ptr()->get_restart_time();
-
-        // a random chi lineout
-        CustomExtraction chi_extraction(c_chi, 1, 15, simParams().L / 2.,
-                                        extraction_origin, m_dt, m_time,
-                                        restart_time, first_step);
-        chi_extraction.execute_query(interpolator, "chi_lineout");
+        // Weyl extraction
+        const int finest = get_gramr_ptr()->finestLevel();
+        amrex::Vector<std::unique_ptr<amrex::MultiFab>> out_weyl(finest + 1);
+
+        // get all the derived info on Weyl
+        auto &lst       = amrex::AmrLevel::get_derive_lst();
+        const auto *rec = lst.get("Weyl4");
+        int state_index = -1, scomp = -1, ncomp_in = -1;
+        rec->getRange(
+            0, state_index, scomp,
+            ncomp_in); // here scomp indicated the starting comp in the state
+                       // (needed to calculate the derived var); similarly
+                       // ncomp_in is the overall number of state comps needed.
+                       // Currently Weyl4 fills all CCZ4 vars, which may be
+                       // problematic.
+
+        // std::cout << "Number of components in Weyl query " << ncomp_in <<
+        // std::endl; std::cout << "Starting component " << scomp << std::endl;
+        // std::cout << "Number of components in Weyl derive " <<
+        // rec->numDerive() << std::endl;
+
+        const int ngrow = 2;
+        const int nout  = rec->numDerive();
+
+        for (int lev = 0; lev <= finest; ++lev)
+        {
+            auto &L           = get_gramr_ptr()->getLevel(lev);
+            const auto &geom  = L.Geom();
+            const auto &state = L.get_new_data(state_index);
+
+            out_weyl[lev] = std::make_unique<amrex::MultiFab>(
+                state.boxArray(), state.DistributionMap(), nout, ngrow);
+
+            // fill ghosts
+            amrex::MultiFab src(state.boxArray(), state.DistributionMap(),
+                                ncomp_in, ngrow);
+            amrex::Copy(src, state, /*scomp=*/scomp, /*dcomp=*/0,
+                        /*ncomp=*/ncomp_in, /*nghost=*/ngrow);
+            src.FillBoundary(geom.periodicity());
+
+            Weyl4::compute_mf(*out_weyl[lev], /*dcomp=*/0, /*ncomp=*/nout, src,
+                              geom, /*time=*/0.0, /*bcrec=*/nullptr,
+                              /*level=*/lev);
+        }
+
+        amrex::Vector<const amrex::MultiFab *> fields;
+        get_gramr_ptr()->convert_derived_multifabs(out_weyl, fields);
+
+        // Perform extraction
+        double m_time         = get_state_data(State_Type).curTime();
+        double m_dt           = get_gramr_ptr()->dtLevel(Level());
+        double m_restart_time = get_gramr_ptr()->get_restart_time();
+
+        WeylExtractionParticle my_extraction(simParams().extraction_params,
+                                             m_dt, m_time, first_step,
+                                             m_restart_time);
+        my_extraction.execute_query(get_bhamr_ptr()->m_weyl_interpolator,
+                                    fields);
     }
 
+    // Custom extraction
+    // if (Level() == 1)
+    // {
+    //     // set up the query and execute it
+    //     std::array<double, AMREX_SPACEDIM> extraction_origin = {
+    //         0.0, simParams().L / 2, -4.0};
+
+    //     double m_time       = get_state_data(State_Type).curTime();
+    //     double m_dt         = get_gramr_ptr()->dtLevel(Level());
+    //     double restart_time = get_gramr_ptr()->get_restart_time();
+
+    //     // a random chi lineout
+    //     CustomExtraction chi_extraction(c_chi, 1, 15, simParams().L / 2.,
+    //                                     extraction_origin, m_dt, m_time,
+    //                                     restart_time, first_step);
+    //     chi_extraction.execute_query(*get_bhamr_ptr()->m_weyl_interpolator,
+    //                                  "chi_lineout");
+    // }
+
     // do puncture tracking on requested level
     if (simParams().puncture_tracking_enabled &&
         Level() == simParams().puncture_tracking_level)

Examples/BinaryBH/BinaryBHLevel.hpp

@@ -50,6 +50,9 @@ class BinaryBHLevel : public GRAMRLevel
     void tag_cells(amrex::TagBoxArray &a_tag_box_array,
                    amrex::Real a_regrid_threshold) final;
 
+    /// Things to do after regrid
+    void specific_post_regrid(int a_lbase, int a_new_finest) override;
+
     //! Things to do after a restart
     void specific_post_restart() override;
 

Examples/BinaryBH/Main_BinaryBH.cpp

@@ -15,6 +15,8 @@
 // Problem specific includes:
 #include "BinaryBHLevel.hpp"
 
+#include "ParticleInterpolators.hpp"
+
 // System includes
 #include <chrono>
 #include <iostream>
@@ -35,7 +37,6 @@ int runGRTeclyn(int /*argc*/, char * /*argv*/[])
     }
 
     GRAMR::set_simulation_parameters(sim_params);
-
     DefaultLevelFactory<BinaryBHLevel> bh_level_bld;
 
 #ifdef USE_TWOPUNCTURES
@@ -49,6 +50,12 @@ int runGRTeclyn(int /*argc*/, char * /*argv*/[])
 
     bh_amr.init(0., sim_params.stop_time);
 
+    // Hand the interpolation pointer to BHAMR
+    ParticleInterpolators<2> weyl_interpolator(
+        sim_params.boundary_params,
+        0); // the derived component usually starts from 0
+    bh_amr.set_interpolator(&weyl_interpolator);
+
     while (
         (bh_amr.okToContinue() != 0) &&
         (bh_amr.levelSteps(0) < sim_params.max_steps ||

Examples/BinaryBH/params_new.txt

@@ -0,0 +1,122 @@
+# GRAMReX BinaryBH parameter file for GitHub Action regression test
+
+#################################################
+# Filesystem parameters
+
+verbosity = 1
+
+output_path = .
+# amr.check_file = chk
+# amr.plot_file = plt
+# amr.restart = chk00000
+# amr.file_name_digits = 5
+
+amr.check_int = 0
+checkpoint_interval = 0
+amr.plot_int = 8
+plot_interval = 8
+amr.derive_plot_vars = ALL
+
+#################################################
+# BH Initial Data parameters
+
+massA = 0.5
+massB = 0.5
+
+offsetA = -10 0 0
+offsetB =  10 0 0
+
+momentumA = 0.0 -0.1 0.0
+momentumB = 0.0  0.1 0.0
+
+#################################################
+# Grid parameters
+
+N_full = 128
+L_full = 64
+
+max_level = 1 # There are (max_level+1) levels
+
+regrid_interval = -1 -1 # no point regridding with no AMR
+
+max_grid_size = 16
+block_factor = 16
+
+# tag_buffer_size = 3
+# num_ghosts = 3
+# center = 8 8 8 # defaults to center of the grid
+
+#################################################
+# Boundary Conditions parameters
+
+# Periodic directions - 0 = false, 1 = true
+isPeriodic = 0 0 0
+# if not periodic, then specify the boundary type
+# 0 = static, 1 = sommerfeld, 2 = reflective
+# (see BoundaryConditions.hpp for details)
+hi_boundary = 1 1 1
+lo_boundary = 1 1 1
+
+# if sommerfeld boundaries selected, must select
+# non zero asymptotic values
+nonzero_asymptotic_vars = chi h11 h22 h33 lapse
+nonzero_asymptotic_values = 1.0 1.0 1.0 1.0 1.0
+
+# if you are using extrapolating BC:
+# extrapolation_order = 1
+
+# extrapolating_vars =
+
+#################################################
+# Evolution parameters
+
+# dt will be dx*dt_multiplier on each grid level
+dt_multiplier = 0.25
+# stop_time = 1.0
+max_steps = 3
+
+# max_spatial_derivative_order = 4 # only 4 implemented at the moment
+
+# nan_check = true
+
+# Lapse evolution
+lapse_advec_coeff = 1.0
+lapse_coeff = 2.0
+lapse_power = 1.0
+
+# Shift evolution
+shift_advec_coeff = 1.0
+shift_Gamma_coeff = 0.75
+eta = 1.0
+
+# CCZ4 parameters
+formulation = 0 # 1 for BSSN, 0 for CCZ4
+kappa1 = 0.1
+kappa2 = 0.
+kappa3 = 1.
+covariantZ4 = 1 # 0: keep kappa1; 1 [default]: replace kappa1 -> kappa1/lapse
+
+# coefficient for KO numerical dissipation
+sigma = 0.5
+
+# min_chi = 1.e-4
+# min_lapse = 1.e-4
+
+#################################################
+# Extraction parameters
+
+# extraction_center = 4 4 0 # defaults to center
+activate_extraction = true
+num_extraction_radii = 1
+extraction_radii = 30.0
+extraction_levels = 0
+num_points_phi = 24
+num_points_theta = 37
+num_modes = 3
+modes = 2 0  2 1  2 2
+
+#################################################
+# Puncture Tracking parameters
+
+puncture_tracking.enabled = 1
+puncture_tracking.level = 0

Source/BlackHoles/BHAMR.hpp

@@ -7,6 +7,8 @@
 #define BHAMR_HPP_
 
 #include "GRAMR.hpp"
+#include "InterpolationQueryParticle.hpp"
+#include "ParticleInterpolators.hpp"
 #include "PunctureTracker.hpp"
 
 #include <AMReX_ParmParse.H>
@@ -21,8 +23,18 @@ template <int num_punctures> class BHAMR : public GRAMR
 {
   private:
     PunctureTracker<num_punctures> m_puncture_tracker;
+    InterpolationQueryParticle *m_query =
+        nullptr; // query used for interpolation
+    bool m_query_populated =
+        false; // flag to identify whether the query has been populated: for
+               // particles that will be fixed we want to do this only once
+               // here!
 
   public:
+
+    ParticleInterpolators<2> *m_weyl_interpolator =
+        nullptr; // weyl interpolator
+
     BHAMR(amrex::LevelBld *a_levelbld) : GRAMR(a_levelbld)
     {
         amrex::ParmParse puncture_tracking_pp("puncture_tracking");
@@ -39,6 +51,42 @@ template <int num_punctures> class BHAMR : public GRAMR
     {
         return m_puncture_tracker;
     }
+
+    // set weyl interpolator
+    void set_interpolator(ParticleInterpolators<2> *a_interpolator)
+    {
+        AMREX_ASSERT(a_interpolator != nullptr);
+
+        m_weyl_interpolator = a_interpolator;
+        m_weyl_interpolator->set_gramr_ptr(this);
+    }
+
+    // set query
+    void set_query(InterpolationQueryParticle &q) override
+    {
+        if (m_query != &q)
+        { // only if the query object changed
+            m_query           = &q;
+            m_query_populated = false;
+        }
+    }
+
+    // populate only once
+    void ensure_query_populated() override
+    {
+        AMREX_ASSERT(m_query);
+        if (!m_query_populated)
+        {
+            m_weyl_interpolator->populate_from_query(*m_query);
+            m_query_populated = true;
+        }
+    }
+
+    // access to a cached query
+    InterpolationQueryParticle *query() override
+    {
+        return m_query ? &*m_query : nullptr;
+    }
 };
 
 #endif /* BHAMR_HPP_ */