Optimize constraint refresh after regrid

AMSS_NCKU_source/bssn_class.C

@@ -65,6 +65,7 @@ bssn_class::bssn_class(double Couranti, double StartTimei, double TotalTimei,
                          xc(0), yc(0), zc(0), xr(0), yr(0), zr(0), trigger(0), dTT(0), dumpid(0),
 #endif
                          a_lev(a_levi), maxl(maxli), decn(decni), maxrex(maxrexi), drex(drexi),
+                         ConstraintRefreshLevels(0),
                          CheckPoint(0)
                          // CheckPoint(0)
 {
@@ -702,6 +703,9 @@ void bssn_class::Initialize()
     }
   }
   GH = new cgh(0, ngfs, Symmetry, pname, checkrun, ErrorMonitor);
+  ConstraintRefreshLevels = new int[GH->levels];
+  for (int il = 0; il < GH->levels; il++)
+    ConstraintRefreshLevels[il] = 0;
   if (checkrun)
     CheckPoint->readcheck_cgh(PhysTime, GH, myrank, nprocs, Symmetry);
   else
@@ -791,6 +795,8 @@ bssn_class::~bssn_class()
   DumpList->clearList();
   ConstraintList->clearList();
 
+  delete[] ConstraintRefreshLevels;
+
   delete phio;
   delete trKo;
   delete gxxo;
@@ -2441,7 +2447,11 @@ void bssn_class::RecursiveStep(int lev)
   if (GH->Regrid_Onelevel(lev, Symmetry, BH_num, Porgbr, Porg0,
                       SynchList_cor, OldStateList, StateList, SynchList_pre,
                       fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
+  {
+  if (ConstraintRefreshLevels)
+    ConstraintRefreshLevels[lev] = 1;
   for (int il = 0; il < GH->levels; il++) { sync_cache_pre[il].invalidate(); sync_cache_cor[il].invalidate(); sync_cache_rp_coarse[il].invalidate(); sync_cache_rp_fine[il].invalidate(); sync_cache_restrict[il].invalidate(); sync_cache_outbd[il].invalidate(); }
+  }
 #endif
 }
 
@@ -7261,7 +7271,7 @@ void bssn_class::Constraint_Out()
     for (int lev = 0; lev < GH->levels; lev++)
     {
       // make sure the data consistent for higher levels
-      if (lev > 0) // if the constrait quantities can be reused from the step rhs calculation
+      if (lev > 0 && ConstraintRefreshLevels && ConstraintRefreshLevels[lev]) // only refresh levels whose grid layout changed after evolution
       {
         double TRK4 = PhysTime;
         double ndeps = numepsb;
@@ -7474,6 +7484,9 @@ void bssn_class::Constraint_Out()
     Interp_Constraint(false);
 
     LastConsOut = 0;
+    if (ConstraintRefreshLevels)
+      for (int lev = 0; lev < GH->levels; lev++)
+        ConstraintRefreshLevels[lev] = 0;
   }
 }
 

AMSS_NCKU_source/bssn_class.h

@@ -48,6 +48,7 @@ public:
        double StartTime, TotalTime;
        double AnasTime, DumpTime, d2DumpTime, CheckTime;
        double LastAnas, LastConsOut;
+       int *ConstraintRefreshLevels;
        double Courant;
        double numepss, numepsb, numepsh;
        int Symmetry;