Port puncture tracker

[mirenradia]

This depends on #2.

[mirenradia]

We have decided to use AMReX's particle capabilities to implement this. I have been working on this this week and what I have done is on the enhancement/puncture_tracker branch. Unfortunately as at 134ff8e, there is a problem that leads to a failed assertion. AMReX seems to think the puncture particles are out of range. I need to check their positions carefully.

[WeiqunZhang]

Is there an inputs file I can try?

One issue might be particle id must be positive. Could you try the following changes?

diff --git a/Source/BlackHoles/PunctureTracker.cpp b/Source/BlackHoles/PunctureTracker.cpp
index 357d4cf..c61d990 100644
--- a/Source/BlackHoles/PunctureTracker.cpp
+++ b/Source/BlackHoles/PunctureTracker.cpp
@@ -96,15 +96,13 @@ void PunctureTracker::set_initial_punctures()
     // use a vector for the write out
     punctures_file.write_time_data_line(get_puncture_vector());
 
+    if (amrex::ParallelDescriptor::MyProc() != 0) return;
+
     // It doesn't matter where we put the puncture particles initially.
     // They will be redistributed later
     const int base_level = 0;
-    for (amrex::MFIter mfi = MakeMFIter(base_level); mfi.isValid(); ++mfi)
     {
-        auto &particle_tile = DefineAndReturnParticleTile(
-            base_level, mfi.index(), mfi.LocalTileIndex());
-        if (mfi.index() != 0 || mfi.LocalTileIndex() != 0)
-            continue;
+        auto &particle_tile = DefineAndReturnParticleTile(base_level, 0, 0);
 
         particle_tile.resize(m_num_punctures);
 
@@ -118,7 +116,8 @@ void PunctureTracker::set_initial_punctures()
                 auto &puncture_particle = particle_tile_data[ipuncture];
                 puncture_particle.pos(idir) =
                     m_puncture_coords[ipuncture][idir];
-                puncture_particle.id() = ipuncture;
+                puncture_particle.id() = ipuncture + 1; // or add an integer data to particels
+                puncture_particle.cpu() = 0;
             }
         }
     }
@@ -148,7 +147,7 @@ void PunctureTracker::redistribute()
             for (int ipunc = 0; ipunc < num_punc_tile; ipunc++)
             {
                 ParticleType &p = punc_particles[ipunc];
-                int ipuncture   = p.id();
+                int ipuncture   = p.id()-1;
                 int linear_idx =
                     m_num_punctures * amrex::ParallelDescriptor::MyProc() +
                     ipuncture;
@@ -277,7 +276,7 @@ void PunctureTracker::execute_tracking(double a_time, double a_restart_time,
             {
                 ParticleType &p = punc_particles[ipunc];
 
-                m_puncture_coords[p.id()] = p.pos();
+                m_puncture_coords[p.id()-1] = p.pos();
             }
         }
     } // ilevel

[WeiqunZhang]

Another patch

diff --git a/Source/BlackHoles/PunctureTracker.cpp b/Source/BlackHoles/PunctureTracker.cpp
index c61d990..fd3288f 100644
--- a/Source/BlackHoles/PunctureTracker.cpp
+++ b/Source/BlackHoles/PunctureTracker.cpp
@@ -132,7 +132,7 @@ void PunctureTracker::redistribute()
     amrex::Vector<int> proc_has_punctures(
         m_num_punctures * amrex::ParallelDescriptor::NProcs(), 0);
 
-    for (int ilevel = 0; ilevel < m_gr_amr->finestLevel(); ilevel++)
+    for (int ilevel = 0; ilevel <= m_gr_amr->finestLevel(); ilevel++)
     {
         if (this->NumberOfParticlesAtLevel(ilevel) == 0L)
         {
@@ -197,7 +197,7 @@ void PunctureTracker::execute_tracking(double a_time, double a_restart_time,
     // Redistribute punctures to the correct grid
     redistribute();
 
-    for (int ilevel = 0; ilevel < m_gr_amr->finestLevel(); ilevel++)
+    for (int ilevel = 0; ilevel <= m_gr_amr->finestLevel(); ilevel++)
     {
         if (this->NumberOfParticlesAtLevel(ilevel) == 0L)
         {

[mirenradia]

Thanks for the patches @WeiqunZhang. I had to change to calling cic_interpolate()(b8c712e). Any idea why my original code with the direct invocation of linear_interpolate_to_particle() doesn't work?

Is there an inputs file I can try?

I've added the necessary parameters to test_params.txt in d0ab14e. However, the resolution is too low with these parameters for the tracking of the punctures to be meaningful. I've been using this as a quick test to check the code doesn't crash but to check if the results are meaningful, I've added

track_punctures = 1
puncture_tracking_level = 7

to the end of params_profile.txt.

One question I have for you is in the case of reflective BCs where the punctures lie right on the edge (face) of the problem domain. Due to finite FP precision, they may drift very slightly outside the problem domain (but still very much within the first layer of ghost cells on the finest level). How do you think we should handle this case? I think I've seen that some of the AMReX particle infrastructure marks particles as invalid if they leave the problem domain.

[WeiqunZhang]

Re: cic_interpolate vs. linear_interpolate_to_particle, I think it's a flaw in amrex. https://github.com/AMReX-Codes/amrex/blob/c49d35e9d45cca08b012c61f1224b459b5264f8c/Src/Particle/AMReX_TracerParticle_mod_K.H#L137 This does not seem to make sense.

for (int comp = start_comp; comp < ncomp; ++comp) {

Since it's called ncomp (number of components), people would expect it to be

for (int comp = start_comp; comp < start_comp+ncomp; ++comp) {

[WeiqunZhang]

Re: reflective BC, we have something called roundoff_lo and roundoff_hi in Geometry. It's guaranteed that a number in the range of [roundoff_lo, roundoff_hi] is considered inside the domain and int(std::floor((x - ProbLo)*InvCellSize)) will return a number in [0,ncell). So maybe you try to make sure all particles are within [roundoff_lo, roundoff_hi].

Currentl roundoff_lo and roundoff_hi are private members. I can add public functions to get them.

[WeiqunZhang]

AMReX-Codes/amrex#4097

[WeiqunZhang]

AMReX-Codes/amrex#4098