[TEST]UPSTREAM: Pick some source changes from 48080d0a97

* Sync new folder structure

.gitignore

@@ -1,6 +1,6 @@
 __pycache__
 GW150914
-GW150914*
+GW150914-origin
 docs
 *.tmp
-.codex
+

.idea/vcs.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="VcsDirectoryMappings">
+    <mapping directory="" vcs="Git" />
+  </component>
+</project>

AMSS_NCKU_Program.py

@@ -126,11 +126,6 @@ setup.generate_AMSSNCKU_input()
 #inputvalue = input()  ## Wait for user input (press Enter) to proceed
 #print()
 
-setup.print_puncture_information()
-
-
-##################################################################
-
 ## Generate AMSS-NCKU program input files based on the configured parameters
 
 print(                                                                                   )
@@ -177,9 +172,6 @@ print( " AMSS-NCKU macro file macrodef.h has been generated. " )
 generate_macrodef.generate_macrodef_fh()
 print( " AMSS-NCKU macro file macrodef.fh has been generated. " )
 
-generate_macrodef.generate_build_config()
-print( " AMSS-NCKU build config AMSS_NCKU_build.mk has been generated. " )
-
 
 ##################################################################
 
@@ -222,11 +214,9 @@ shutil.copytree(AMSS_NCKU_source_path, AMSS_NCKU_source_copy)
 
 macrodef_h_path  = os.path.join(File_directory, "macrodef.h") 
 macrodef_fh_path = os.path.join(File_directory, "macrodef.fh") 
-build_config_path = os.path.join(File_directory, "AMSS_NCKU_build.mk")
 
 shutil.copy2(macrodef_h_path,  AMSS_NCKU_source_copy)
 shutil.copy2(macrodef_fh_path, AMSS_NCKU_source_copy)
-shutil.copy2(build_config_path, AMSS_NCKU_source_copy)
 
 # Notes on copying files:
 # shutil.copy2 preserves file metadata such as modification time.
@@ -317,7 +307,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
     
     import generate_TwoPuncture_input
     
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
+    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input(numerical_grid.puncture_data)
     
     print(                                                                                              )
     print( " The input parfile for the TwoPunctureABE executable has been generated. " )
@@ -359,7 +349,7 @@ if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ):
 
 import renew_puncture_parameter
     
-renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
+renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory, numerical_grid.puncture_data)
 
 
 ## Generated AMSS-NCKU input filename

AMSS_NCKU_Program_Plot.py

@@ -1,100 +0,0 @@
-##################################################################
-##
-## AMSS-NCKU Plot-Only Restart Script
-## Author: Xiaoqu / Claude
-## 2026/05/12
-##
-## This script checks for existing output data from AMSS_NCKU_Program.py.
-## If data exists, it skips all computation and goes directly to plotting,
-## saving time when plotting was interrupted.
-## If no data is found, it exits with a message.
-##
-##################################################################
-
-## Guard against re-execution by multiprocessing child processes.
-if __name__ != '__main__':
-    import sys as _sys
-    _sys.exit(0)
-
-
-import os
-import sys
-
-import AMSS_NCKU_Input as input_data
-
-##################################################################
-
-## Construct paths from input configuration
-File_directory         = os.path.join(input_data.File_directory)
-output_directory       = os.path.join(File_directory, "AMSS_NCKU_output")
-binary_results_directory = os.path.join(output_directory, input_data.Output_directory)
-figure_directory       = os.path.join(File_directory, "figure")
-
-##################################################################
-
-## Check whether the required output data files exist
-
-required_files = [
-    os.path.join(binary_results_directory, "bssn_BH.dat"),
-    os.path.join(binary_results_directory, "bssn_ADMQs.dat"),
-    os.path.join(binary_results_directory, "bssn_psi4.dat"),
-    os.path.join(binary_results_directory, "bssn_constraint.dat"),
-]
-
-missing_files = [f for f in required_files if not os.path.exists(f)]
-
-if missing_files:
-    print(" No existing AMSS_NCKU_Program.py output data found. ")
-    print(" The following required files are missing: ")
-    for f in missing_files:
-        print(f"   {f}")
-    print()
-    print(" Please run AMSS_NCKU_Program.py first to generate the simulation data. ")
-    print(" Exiting. ")
-    sys.exit(1)
-
-print(" Found existing AMSS_NCKU_Program.py output data. "          )
-print(" Skipping all computation and going directly to plotting. " )
-print()
-
-## Ensure the figure directory exists (it should, but be safe)
-os.makedirs(figure_directory, exist_ok=True)
-
-##################################################################
-
-## Plot the AMSS-NCKU program results
-
-import plot_xiaoqu
-import plot_GW_strain_amplitude_xiaoqu
-from parallel_plot_helper import run_plot_tasks_parallel
-
-plot_tasks = []
-
-## Plot black hole trajectory
-plot_tasks.append((plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory)))
-plot_tasks.append((plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory)))
-
-## Plot black hole separation vs. time
-plot_tasks.append((plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory)))
-
-## Plot gravitational waveforms (psi4 and strain amplitude)
-for i in range(input_data.Detector_Number):
-    plot_tasks.append((plot_xiaoqu.generate_gravitational_wave_psi4_plot, (binary_results_directory, figure_directory, i)))
-    plot_tasks.append((plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i)))
-
-## Plot ADM mass evolution
-for i in range(input_data.Detector_Number):
-    plot_tasks.append((plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i)))
-
-## Plot Hamiltonian constraint violation over time
-for i in range(input_data.grid_level):
-    plot_tasks.append((plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i)))
-
-run_plot_tasks_parallel(plot_tasks)
-
-## Plot stored binary data (runs serially, not in the parallel pool)
-plot_xiaoqu.generate_binary_data_plot(binary_results_directory, figure_directory)
-
-print()
-print(" Plotting completed successfully. ")
-print()

AMSS_NCKU_source/BSSN/bssnEM_class.C

@@ -258,8 +258,6 @@ void bssnEM_class::Initialize()
     PhysTime = StartTime;
     Setup_Black_Hole_position();
   }
-
-  setup_transfer_caches();
 }
 
 //================================================================================================

AMSS_NCKU_source/BSSN/bssn_class.C

@@ -299,28 +299,6 @@ bssn_class::bssn_class(double Couranti, double StartTimei, double TotalTimei,
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
 
-  // Derived classes override Initialize(), so ownership-sensitive members must
-  // be in a known state before any specialized setup path runs.
-  GH = 0;
-  SH = 0;
-  PhysTime = 0.0;
-  BH_num = 0;
-  BH_num_input = 0;
-  Porg0 = 0;
-  Porgbr = 0;
-  Porg = 0;
-  Porg1 = 0;
-  Porg_rhs = 0;
-  Mass = 0;
-  Pmom = 0;
-  Spin = 0;
-  sync_cache_pre = 0;
-  sync_cache_cor = 0;
-  sync_cache_rp_coarse = 0;
-  sync_cache_rp_fine = 0;
-  sync_cache_restrict = 0;
-  sync_cache_outbd = 0;
-
   // setup Monitors
   {
     stringstream a_stream;
@@ -1008,7 +986,13 @@ void bssn_class::Initialize()
     Setup_Black_Hole_position();
   }
 
-  setup_transfer_caches();
+  // Initialize sync caches (per-level, for predictor and corrector)
+  sync_cache_pre = new Parallel::SyncCache[GH->levels];
+  sync_cache_cor = new Parallel::SyncCache[GH->levels];
+  sync_cache_rp_coarse = new Parallel::SyncCache[GH->levels];
+  sync_cache_rp_fine = new Parallel::SyncCache[GH->levels];
+  sync_cache_restrict = new Parallel::SyncCache[GH->levels];
+  sync_cache_outbd = new Parallel::SyncCache[GH->levels];
 }
 
 //================================================================================================
@@ -1263,7 +1247,30 @@ bssn_class::~bssn_class()
 #endif
 
   // Destroy sync caches before GH
-  destroy_transfer_caches();
+  if (sync_cache_pre)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_pre[i].destroy();
+    delete[] sync_cache_pre;
+  }
+  if (sync_cache_cor)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_cor[i].destroy();
+    delete[] sync_cache_cor;
+  }
+  if (sync_cache_rp_coarse)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_rp_coarse[i].destroy();
+    delete[] sync_cache_rp_coarse;
+  }
+  if (sync_cache_rp_fine)
+  {
+    for (int i = 0; i < GH->levels; i++)
+      sync_cache_rp_fine[i].destroy();
+    delete[] sync_cache_rp_fine;
+  }
 
   delete GH;
 #ifdef WithShell
@@ -2482,7 +2489,7 @@ void bssn_class::Evolve(int Steps)
     GH->Regrid(Symmetry, BH_num, Porgbr, Porg0,
                SynchList_cor, OldStateList, StateList, SynchList_pre,
                fgt(PhysTime - dT_mon, StartTime, dT_mon / 2), ErrorMonitor);
-    invalidate_transfer_caches();
+    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(); }
     STEP_TIMER_ADD(TB_REGRID, timer_regrid);
 #endif
 
@@ -2723,7 +2730,7 @@ void bssn_class::RecursiveStep(int lev)
   {
   if (ConstraintRefreshLevels)
     ConstraintRefreshLevels[lev] = 1;
-  invalidate_transfer_caches();
+  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(); }
   }
   STEP_TIMER_ADD(TB_REGRID, timer_regrid_onelevel);
 #endif
@@ -2904,7 +2911,7 @@ void bssn_class::ParallelStep()
   if (GH->Regrid_Onelevel(GH->mylev, Symmetry, BH_num, Porgbr, Porg0,
                       SynchList_cor, OldStateList, StateList, SynchList_pre,
                       fgt(PhysTime - dT_lev, StartTime, dT_lev / 2), ErrorMonitor))
-    invalidate_transfer_caches();
+  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
 }
 
@@ -3068,10 +3075,10 @@ void bssn_class::ParallelStep()
       if (lev + 1 >= GH->movls)
       {
         //	       GH->Regrid_Onelevel_aux(lev,Symmetry,BH_num,Porgbr,Porg0,
-          if (GH->Regrid_Onelevel(lev + 1, Symmetry, BH_num, Porgbr, Porg0,
+        if (GH->Regrid_Onelevel(lev + 1, Symmetry, BH_num, Porgbr, Porg0,
-                              SynchList_cor, OldStateList, StateList, SynchList_pre,
+                            SynchList_cor, OldStateList, StateList, SynchList_pre,
-                              fgt(PhysTime - dT_levp1, StartTime, dT_levp1 / 2), ErrorMonitor))
+                            fgt(PhysTime - dT_levp1, StartTime, dT_levp1 / 2), ErrorMonitor))
-            invalidate_transfer_caches();
+        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(); }
 
         //               a_stream.clear();
         //               a_stream.str("");
@@ -3086,7 +3093,7 @@ void bssn_class::ParallelStep()
       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))
-        invalidate_transfer_caches();
+      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(); }
 
       //               a_stream.clear();
       //               a_stream.str("");
@@ -3105,7 +3112,7 @@ void bssn_class::ParallelStep()
           if (GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
                               SynchList_cor, OldStateList, StateList, SynchList_pre,
                               fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
-            invalidate_transfer_caches();
+          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(); }
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -3121,7 +3128,7 @@ void bssn_class::ParallelStep()
           if (GH->Regrid_Onelevel(lev - 1, Symmetry, BH_num, Porgbr, Porg0,
                               SynchList_cor, OldStateList, StateList, SynchList_pre,
                               fgt(PhysTime - dT_lev, StartTime, dT_levm1 / 2), ErrorMonitor))
-            invalidate_transfer_caches();
+          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(); }
 
           //               a_stream.clear();
           //               a_stream.str("");
@@ -3652,7 +3659,7 @@ void bssn_class::Step(int lev, int YN)
 
   STEP_TIMER_DECL(timer_predictor_sync);
   Parallel::AsyncSyncState async_pre;
-  sync_predictor_start(lev, SynchList_pre, async_pre);
+  Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
 
 #ifdef WithShell
   if (lev == 0)
@@ -3671,7 +3678,7 @@ void bssn_class::Step(int lev, int YN)
     }
   }
 #endif
-  sync_predictor_finish(lev, async_pre, SynchList_pre);
+  Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry);
 
 #ifdef WithShell
   // Complete non-blocking error reduction and check
@@ -4017,7 +4024,7 @@ void bssn_class::Step(int lev, int YN)
 
     STEP_TIMER_DECL(timer_corrector_sync);
     Parallel::AsyncSyncState async_cor;
-    sync_corrector_start(lev, SynchList_cor, async_cor);
+    Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
 
 #ifdef WithShell
     if (lev == 0)
@@ -4036,7 +4043,7 @@ void bssn_class::Step(int lev, int YN)
       }
     }
 #endif
-    sync_corrector_finish(lev, async_cor, SynchList_cor);
+    Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry);
 
 #ifdef WithShell
     // Complete non-blocking error reduction and check
@@ -4525,7 +4532,7 @@ void bssn_class::Step(int lev, int YN)
 #endif
 
   Parallel::AsyncSyncState async_pre;
-  sync_predictor_start(lev, SynchList_pre, async_pre);
+  Parallel::Sync_start(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev], async_pre);
 
 #ifdef WithShell
   if (lev == 0)
@@ -4544,7 +4551,7 @@ void bssn_class::Step(int lev, int YN)
     }
   }
 #endif
-  sync_predictor_finish(lev, async_pre, SynchList_pre);
+  Parallel::Sync_finish(sync_cache_pre[lev], async_pre, SynchList_pre, Symmetry);
 
 #ifdef WithShell
   // Complete non-blocking error reduction and check
@@ -4873,7 +4880,7 @@ void bssn_class::Step(int lev, int YN)
 #endif
 
     Parallel::AsyncSyncState async_cor;
-    sync_corrector_start(lev, SynchList_cor, async_cor);
+    Parallel::Sync_start(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev], async_cor);
 
 #ifdef WithShell
     if (lev == 0)
@@ -4892,7 +4899,7 @@ void bssn_class::Step(int lev, int YN)
       }
     }
 #endif
-    sync_corrector_finish(lev, async_cor, SynchList_cor);
+    Parallel::Sync_finish(sync_cache_cor[lev], async_cor, SynchList_cor, Symmetry);
 
 #ifdef WithShell
     // Complete non-blocking error reduction and check
@@ -5284,7 +5291,7 @@ void bssn_class::Step(int lev, int YN)
 
   //   misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor sync");
 
-  sync_evolution(lev, SynchList_pre, sync_cache_pre);
+  Parallel::Sync_cached(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev]);
 
   // Complete non-blocking error reduction and check
   MPI_Wait(&err_req, MPI_STATUS_IGNORE);
@@ -5485,7 +5492,7 @@ void bssn_class::Step(int lev, int YN)
 
     //    misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector sync");
 
-    sync_evolution(lev, SynchList_cor, sync_cache_cor);
+    Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev]);
 
     //    misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after Corrector sync");
 
@@ -6074,92 +6081,6 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //
 // SynchList_cor  old -----------
 {
-#if (ABEtype == 1 || ABEtype == 2)
-#if (PSTR == 1 || PSTR == 2)
-//  stringstream a_stream;
-//  a_stream.setf(ios::left);
-#endif
-
-  if (lev > 0)
-  {
-    MyList<Patch> *Pp, *Ppc;
-    if (lev > trfls && YN == 0)
-    {
-      Pp = GH->PatL[lev - 1];
-      while (Pp)
-      {
-        if (BB)
-          Parallel::prepare_inter_time_level(Pp->data, SL, OL, corL,
-                                             SynchList_pre, 0);
-        else
-          Parallel::prepare_inter_time_level(Pp->data, SL, OL,
-                                             SynchList_pre, 0);
-        Pp = Pp->next;
-      }
-
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], SynchList_pre, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, GH->bdsul[lev], Symmetry);
-#endif
-    }
-    else
-    {
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], SL, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SL, SL, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SL, SL, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->bdsul[lev], Symmetry);
-#endif
-    }
-
-    Parallel::Sync(GH->PatL[lev], SL, Symmetry);
-  }
-  return;
-#endif
   STEP_TIMER_DECL(timer_restrict_prolong);
 #if (PSTR == 1 || PSTR == 2)
 //  stringstream a_stream;
@@ -6202,7 +6123,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 #endif
 
 #if (RPB == 0)
-      restrict_evolution(lev, SL, SynchList_pre);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SynchList_pre,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
@@ -6215,7 +6136,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
-      sync_evolution(lev - 1, SynchList_pre, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -6226,7 +6147,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SynchList_pre, SL);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #endif
@@ -6253,7 +6174,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 #endif
 
 #if (RPB == 0)
-      restrict_evolution(lev, SL, SL);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SL,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
@@ -6266,7 +6187,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
-      sync_evolution(lev - 1, SL, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -6277,7 +6198,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SL, SL);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #endif
@@ -6294,7 +6215,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 #endif
     }
 
-    sync_evolution(lev, SL, sync_cache_rp_fine);
+    Parallel::Sync_cached(GH->PatL[lev], SL, Symmetry, sync_cache_rp_fine[lev]);
 
 #if (PSTR == 1 || PSTR == 2)
 //    a_stream.clear();
@@ -6323,91 +6244,6 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
 //
 // SynchList_cor  old -----------
 {
-#if (ABEtype == 1 || ABEtype == 2)
-  if (lev >= GH->levels - 1)
-    return;
-  lev = lev + 1;
-
-  if (lev > 0)
-  {
-    MyList<Patch> *Pp, *Ppc;
-    if (lev > trfls && YN == 0)
-    {
-      Pp = GH->PatL[lev - 1];
-      while (Pp)
-      {
-        if (BB)
-          Parallel::prepare_inter_time_level(Pp->data, SL, OL, corL,
-                                             SynchList_pre, 0);
-        else
-          Parallel::prepare_inter_time_level(Pp->data, SL, OL,
-                                             SynchList_pre, 0);
-        Pp = Pp->next;
-      }
-
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], SynchList_pre, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SL, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, GH->bdsul[lev], Symmetry);
-#endif
-    }
-    else
-    {
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], SL, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SL, SL, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SL, SL, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->bdsul[lev], Symmetry);
-#endif
-    }
-
-    Parallel::Sync(GH->PatL[lev], SL, Symmetry);
-  }
-  return;
-#endif
   STEP_TIMER_DECL(timer_restrict_prolong);
   //  misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"starting RestrictProlong_aux");
 
@@ -6433,17 +6269,17 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       }
 
 #if (RPB == 0)
-      restrict_evolution(lev, SL, SynchList_pre);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SynchList_pre,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
-      sync_evolution(lev - 1, SynchList_pre, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SynchList_pre, SL);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SL, Symmetry);
 #endif
@@ -6455,17 +6291,17 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
     else // no time refinement levels and for all same time levels
     {
 #if (RPB == 0)
-      restrict_evolution(lev, SL, SL);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SL,SL,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
 #endif
 
-      sync_evolution(lev - 1, SL, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SL, SL);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, Symmetry);
 #endif
@@ -6475,11 +6311,7 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
 #endif
     }
 
-#if (ABEtype == 1 || ABEtype == 2)
+    Parallel::Sync_cached(GH->PatL[lev], SL, Symmetry, sync_cache_rp_fine[lev]);
-    Parallel::Sync(GH->PatL[lev], SL, Symmetry);
-#else
-    sync_evolution(lev, SL, sync_cache_rp_fine);
-#endif
   }
   STEP_TIMER_ADD(TB_RESTRICT_PROLONG, timer_restrict_prolong);
 }
@@ -6492,93 +6324,8 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
 
 void bssn_class::RestrictProlong(int lev, int YN, bool BB)
 {
-  double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
-#if (ABEtype == 1 || ABEtype == 2)
-  if (lev > 0)
-  {
-    MyList<Patch> *Pp, *Ppc;
-    if (lev > trfls && YN == 0)
-    {
-      if (myrank == 0)
-        cout << "/=: " << GH->Lt[lev - 1] << "," << GH->Lt[lev] + dT_lev << endl;
-      Pp = GH->PatL[lev - 1];
-      while (Pp)
-      {
-        if (BB)
-          Parallel::prepare_inter_time_level(Pp->data, StateList, OldStateList, SynchList_cor,
-                                             SynchList_pre, 0);
-        else
-          Parallel::prepare_inter_time_level(Pp->data, StateList, OldStateList,
-                                             SynchList_pre, 0);
-        Pp = Pp->next;
-      }
-
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], SynchList_pre, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, SynchList_pre, SynchList_cor, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, GH->bdsul[lev], Symmetry);
-#endif
-    }
-    else
-    {
-      if (myrank == 0)
-        cout << "===: " << GH->Lt[lev - 1] << "," << GH->Lt[lev] + dT_lev << endl;
-#if (RPB == 0)
-      Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, Symmetry);
-#elif (RPB == 1)
-      Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, GH->rsul[lev], Symmetry);
-#endif
-
-      Parallel::Sync(GH->PatL[lev - 1], StateList, Symmetry);
-
-#if (RPB == 0)
-      Ppc = GH->PatL[lev - 1];
-      while (Ppc)
-      {
-        Pp = GH->PatL[lev];
-        while (Pp)
-        {
-#if (MIXOUTB == 0)
-          Parallel::OutBdLow2Hi(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
-#elif (MIXOUTB == 1)
-          Parallel::OutBdLow2Himix(Ppc->data, Pp->data, StateList, SynchList_cor, Symmetry);
-#endif
-          Pp = Pp->next;
-        }
-        Ppc = Ppc->next;
-      }
-#elif (RPB == 1)
-      Parallel::OutBdLow2Hi_bam(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, GH->bdsul[lev], Symmetry);
-#endif
-    }
-
-    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
-  }
-  return;
-#endif
   STEP_TIMER_DECL(timer_restrict_prolong);
+  double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
   // we assume  for fine
   // SynchList_cor 1   -----------
   //
@@ -6611,17 +6358,17 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       }
 
 #if (RPB == 0)
-      restrict_evolution(lev, SynchList_cor, SynchList_pre);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_cor,SynchList_pre,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
-      sync_evolution(lev - 1, SynchList_pre, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SynchList_pre, SynchList_cor);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #endif
@@ -6635,17 +6382,17 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       if (myrank == 0)
         cout << "===: " << GH->Lt[lev - 1] << "," << GH->Lt[lev] + dT_lev << endl;
 #if (RPB == 0)
-      restrict_evolution(lev, SynchList_cor, StateList);
+      Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, Symmetry, sync_cache_restrict[lev]);
 #elif (RPB == 1)
       //       Parallel::Restrict_bam(GH->PatL[lev-1],GH->PatL[lev],SynchList_cor,StateList,Symmetry);
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, GH->rsul[lev], Symmetry);
 #endif
 
-      sync_evolution(lev - 1, StateList, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, StateList, SynchList_cor);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #endif
@@ -6655,7 +6402,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
 #endif
     }
 
-    sync_evolution(lev, SynchList_cor, sync_cache_rp_fine);
+    Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_rp_fine[lev]);
   }
   STEP_TIMER_ADD(TB_RESTRICT_PROLONG, timer_restrict_prolong);
 }
@@ -6687,7 +6434,7 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, SynchList_pre, SynchList_cor);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], SynchList_pre, SynchList_cor, Symmetry);
 #endif
@@ -6700,7 +6447,7 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
     {
 #if (RPB == 0)
 #if (MIXOUTB == 0)
-      outbdlow2hi_evolution(lev, StateList, SynchList_cor);
+      Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry, sync_cache_outbd[lev]);
 #elif (MIXOUTB == 1)
       Parallel::OutBdLow2Himix(GH->PatL[lev - 1], GH->PatL[lev], StateList, SynchList_cor, Symmetry);
 #endif
@@ -6719,10 +6466,10 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
 #else
       Parallel::Restrict_after(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, Symmetry);
 #endif
-      sync_evolution(lev - 1, StateList, sync_cache_rp_coarse);
+      Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
     }
 
-    sync_evolution(lev, SynchList_cor, sync_cache_rp_fine);
+    Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_rp_fine[lev]);
   }
 }
 #undef MIXOUTB
@@ -7452,169 +7199,6 @@ void bssn_class::compute_Porg_rhs(double **BH_PS, double **BH_RHS, var *forx, va
     }
   }
 }
-
-bool bssn_class::use_transfer_cache() const
-{
-#if BSSN_USE_TRANSFER_CACHE
-  return true;
-#else
-  return false;
-#endif
-}
-
-void bssn_class::setup_transfer_caches()
-{
-  sync_cache_pre = 0;
-  sync_cache_cor = 0;
-  sync_cache_rp_coarse = 0;
-  sync_cache_rp_fine = 0;
-  sync_cache_restrict = 0;
-  sync_cache_outbd = 0;
-
-  if (!use_transfer_cache() || !GH)
-    return;
-
-  sync_cache_pre = new Parallel::SyncCache[GH->levels];
-  sync_cache_cor = new Parallel::SyncCache[GH->levels];
-  sync_cache_rp_coarse = new Parallel::SyncCache[GH->levels];
-  sync_cache_rp_fine = new Parallel::SyncCache[GH->levels];
-  sync_cache_restrict = new Parallel::SyncCache[GH->levels];
-  sync_cache_outbd = new Parallel::SyncCache[GH->levels];
-}
-
-void bssn_class::invalidate_transfer_caches()
-{
-  if (!use_transfer_cache() || !GH || !sync_cache_pre || !sync_cache_cor ||
-      !sync_cache_rp_coarse || !sync_cache_rp_fine || !sync_cache_restrict || !sync_cache_outbd)
-    return;
-
-  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();
-  }
-}
-
-void bssn_class::destroy_transfer_caches()
-{
-  if (sync_cache_pre)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_pre[i].destroy();
-    delete[] sync_cache_pre;
-    sync_cache_pre = 0;
-  }
-  if (sync_cache_cor)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_cor[i].destroy();
-    delete[] sync_cache_cor;
-    sync_cache_cor = 0;
-  }
-  if (sync_cache_rp_coarse)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_rp_coarse[i].destroy();
-    delete[] sync_cache_rp_coarse;
-    sync_cache_rp_coarse = 0;
-  }
-  if (sync_cache_rp_fine)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_rp_fine[i].destroy();
-    delete[] sync_cache_rp_fine;
-    sync_cache_rp_fine = 0;
-  }
-  if (sync_cache_restrict)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_restrict[i].destroy();
-    delete[] sync_cache_restrict;
-    sync_cache_restrict = 0;
-  }
-  if (sync_cache_outbd)
-  {
-    if (use_transfer_cache() && GH)
-      for (int i = 0; i < GH->levels; i++)
-        sync_cache_outbd[i].destroy();
-    delete[] sync_cache_outbd;
-    sync_cache_outbd = 0;
-  }
-}
-
-void bssn_class::sync_predictor_start(int lev, MyList<var> *VarList, Parallel::AsyncSyncState &async_state)
-{
-  if (use_transfer_cache())
-    Parallel::Sync_start(GH->PatL[lev], VarList, Symmetry, sync_cache_pre[lev], async_state);
-  else
-    Parallel::Sync(GH->PatL[lev], VarList, Symmetry);
-}
-
-void bssn_class::sync_predictor_finish(int lev, Parallel::AsyncSyncState &async_state, MyList<var> *VarList)
-{
-  if (use_transfer_cache())
-    Parallel::Sync_finish(sync_cache_pre[lev], async_state, VarList, Symmetry);
-}
-
-void bssn_class::sync_corrector_start(int lev, MyList<var> *VarList, Parallel::AsyncSyncState &async_state)
-{
-  if (use_transfer_cache())
-    Parallel::Sync_start(GH->PatL[lev], VarList, Symmetry, sync_cache_cor[lev], async_state);
-  else
-    Parallel::Sync(GH->PatL[lev], VarList, Symmetry);
-}
-
-void bssn_class::sync_corrector_finish(int lev, Parallel::AsyncSyncState &async_state, MyList<var> *VarList)
-{
-  if (use_transfer_cache())
-    Parallel::Sync_finish(sync_cache_cor[lev], async_state, VarList, Symmetry);
-}
-
-void bssn_class::sync_evolution(int lev, MyList<var> *VarList, Parallel::SyncCache *cache_array)
-{
-  if (use_transfer_cache() && cache_array)
-    Parallel::Sync_cached(GH->PatL[lev], VarList, Symmetry, cache_array[lev]);
-  else
-    Parallel::Sync(GH->PatL[lev], VarList, Symmetry);
-}
-
-void bssn_class::restrict_evolution(int lev, MyList<var> *src_var_list, MyList<var> *dst_var_list)
-{
-  if (use_transfer_cache())
-    Parallel::Restrict_cached(GH->PatL[lev - 1], GH->PatL[lev], src_var_list, dst_var_list, Symmetry, sync_cache_restrict[lev]);
-  else
-    Parallel::Restrict(GH->PatL[lev - 1], GH->PatL[lev], src_var_list, dst_var_list, Symmetry);
-}
-
-void bssn_class::outbdlow2hi_evolution(int lev, MyList<var> *src_var_list, MyList<var> *dst_var_list)
-{
-  if (use_transfer_cache())
-  {
-    Parallel::OutBdLow2Hi_cached(GH->PatL[lev - 1], GH->PatL[lev], src_var_list, dst_var_list, Symmetry, sync_cache_outbd[lev]);
-    return;
-  }
-
-  MyList<Patch> *Ppc = GH->PatL[lev - 1];
-  while (Ppc)
-  {
-    MyList<Patch> *Pp = GH->PatL[lev];
-    while (Pp)
-    {
-      Parallel::OutBdLow2Hi(Ppc->data, Pp->data, src_var_list, dst_var_list, Symmetry);
-      Pp = Pp->next;
-    }
-    Ppc = Ppc->next;
-  }
-}
 #endif
 
 //================================================================================================

AMSS_NCKU_source/BSSN/bssn_class.h

@@ -33,14 +33,6 @@ using namespace std;
 
 extern void setpbh(int iBHN, double **iPBH, double *iMass, int rBHN);
 
-#ifndef BSSN_USE_TRANSFER_CACHE
-#define BSSN_USE_TRANSFER_CACHE 1
-#endif
-
-#ifndef BSSN_USE_ESCALAR_C_KERNEL
-#define BSSN_USE_ESCALAR_C_KERNEL 1
-#endif
-
 class bssn_class
 {
 public:
@@ -179,17 +171,6 @@ public:
        void testOutBd();
        
        bool check_Stdin_Abort(); 
-       bool use_transfer_cache() const;
-       void setup_transfer_caches();
-       void invalidate_transfer_caches();
-       void destroy_transfer_caches();
-       void sync_predictor_start(int lev, MyList<var> *VarList, Parallel::AsyncSyncState &async_state);
-       void sync_predictor_finish(int lev, Parallel::AsyncSyncState &async_state, MyList<var> *VarList);
-       void sync_corrector_start(int lev, MyList<var> *VarList, Parallel::AsyncSyncState &async_state);
-       void sync_corrector_finish(int lev, Parallel::AsyncSyncState &async_state, MyList<var> *VarList);
-       void sync_evolution(int lev, MyList<var> *VarList, Parallel::SyncCache *cache_array = 0);
-       void restrict_evolution(int lev, MyList<var> *src_var_list, MyList<var> *dst_var_list);
-       void outbdlow2hi_evolution(int lev, MyList<var> *src_var_list, MyList<var> *dst_var_list);
 
        virtual void Setup_Initial_Data_Cao();
        virtual void Setup_Initial_Data_Lousto();

AMSS_NCKU_source/BSSN/bssn_rhs.h

@@ -67,27 +67,6 @@ extern "C"
                                int &, int &, double &, int &);
 }
 
-int f_compute_rhs_bssn_escalar_c(int *, double &, double *, double *, double *,                                                      // ex,T,X,Y,Z
-                                 double *, double *,                                                                                 // chi, trK
-                                 double *, double *, double *, double *, double *, double *,                                         // gij
-                                 double *, double *, double *, double *, double *, double *,                                         // Aij
-                                 double *, double *, double *,                                                                       // Gam
-                                 double *, double *, double *, double *, double *, double *, double *,                               // Gauge
-                                 double *, double *,                                                                                 // Sphi, Spi
-                                 double *, double *,                                                                                 // chi, trK
-                                 double *, double *, double *, double *, double *, double *,                                         // gij
-                                 double *, double *, double *, double *, double *, double *,                                         // Aij
-                                 double *, double *, double *,                                                                       // Gam
-                                 double *, double *, double *, double *, double *, double *, double *,                               // Gauge
-                                 double *, double *,                                                                                 // Sphi, Spi
-                                 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, // stress-energy
-                                 double *, double *, double *, double *, double *, double *,                                         // Christoffel
-                                 double *, double *, double *, double *, double *, double *,                                         // Christoffel
-                                 double *, double *, double *, double *, double *, double *,                                         // Christoffel
-                                 double *, double *, double *, double *, double *, double *,                                         // Ricci
-                                 double *, double *, double *, double *, double *, double *, double *,                               // constraint violation
-                                 int &, int &, double &, int &);
-
 extern "C"
 {
         int f_compute_rhs_bssn_ss(int *, double &, double *, double *, double *,                                                      // ex,T,rho,sigma,R
@@ -262,31 +241,4 @@ extern "C"
                                      double *);
 } // FR_cons
 
-// BSSN-EM C kernel (replaces empart.f90 + bssn_rhs.f90 for BSSN+Maxwell)
-int f_compute_rhs_bssn_em_c(int *, double &, double *, double *, double *,
-                             double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *,
-                             double *, double *,
-                             double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *,
-                             double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             double *, double *, double *, double *, double *, double *,
-                             int &, int &, double &, int &);
-
 #endif /* BSSN_H */

AMSS_NCKU_source/BSSN/bssn_rhs_c.C

@@ -1075,10 +1075,6 @@ int f_compute_rhs_bssn(int *ex, double &T,
         }
         #endif
 
-        #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
-        fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
-        #endif
-
         for (int i = 0; i < all; i += 1) {
             #if (GAUGE == 0)
             betax_rhs[i] = FF * dtSfx[i];
@@ -1164,17 +1160,11 @@ int f_compute_rhs_bssn(int *ex, double &T,
         lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps);
         lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
         lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps);
-        #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
         lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps);
-        #endif
         lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
-        #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
         lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps);
-        #endif
         lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps);
-        #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
         lopsided_kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
-        #endif
         lopsided_kodis(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,eps);
         lopsided_kodis(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,eps);
         lopsided_kodis(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,eps);

AMSS_NCKU_source/BSSN/lopsided_c.C

@@ -0,0 +1,255 @@
+#include "tool.h"
+/*
+ * 你需要提供 symmetry_bd 的 C 版本（或 Fortran 绑到 C 的接口）。
+ * Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+ *
+ * 约定：
+ *   nghost = 3
+ *   ex[3]  = {ex1,ex2,ex3}
+ *   f      = 原始网格 (ex1*ex2*ex3)
+ *   fh     = 扩展网格 ((ex1+3)*(ex2+3)*(ex3+3))，对应 Fortran 的 (-2:ex1, ...)
+ *   SoA[3] = 输入参数
+ */
+void lopsided(const int ex[3],
+              const double *X, const double *Y, const double *Z,
+              const double *f, double *f_rhs,
+              const double *Sfx, const double *Sfy, const double *Sfz,
+              int Symmetry, const double SoA[3])
+{
+    const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
+    const double TWO = 2.0, F6 = 6.0, F18 = 18.0;
+    const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
+
+    const int NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2;
+    (void)OCTANT; // 这里和 Fortran 一样只是定义了不用也没关系
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    // 对应 Fortran: dX = X(2)-X(1)  （Fortran 1-based）
+    // C: X[1]-X[0]
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    // Fortran 里算了 d2dx/d2dy/d2dz 但本 subroutine 里没用到（保持一致也算出来）
+    const double d2dx  = ONE / TWO / dX;
+    const double d2dy  = ONE / TWO / dY;
+    const double d2dz  = ONE / TWO / dZ;
+    (void)d2dx; (void)d2dy; (void)d2dz;
+
+    // Fortran:
+    // imax = ex(1); jmax = ex(2); kmax = ex(3)
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    // Fortran:
+    // imin=jmin=kmin=1; 若满足对称条件则设为 -2
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
+
+    // 分配 fh：大小 (ex1+3)*(ex2+3)*(ex3+3)
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return; // 内存不足：直接返回（你也可以改成 abort/报错）
+
+    // Fortran: call symmetry_bd(3,ex,f,fh,SoA)
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    /*
+     * Fortran 主循环：
+     * do k=1,ex(3)-1
+     * do j=1,ex(2)-1
+     * do i=1,ex(1)-1
+     *
+     * 转成 C 0-based：
+     * k0 = 0..ex3-2, j0 = 0..ex2-2, i0 = 0..ex1-2
+     *
+     * 并且 Fortran 里的 i/j/k 在 fh 访问时，仍然是 Fortran 索引值：
+     * iF=i0+1, jF=j0+1, kF=k0+1
+     */
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                // ---------------- x direction ----------------
+                const double sfx = Sfx[p];
+                if (sfx > ZEO) {
+                    // Fortran: if(i+3 <= imax)
+                    // iF+3 <= ex1  <=> i0+4 <= ex1 <=> i0 <= ex1-4
+                    if (i0 <= ex1 - 4) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                    // elseif(i+2 <= imax)  <=> i0 <= ex1-3
+                    else if (i0 <= ex1 - 3) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i+1 <= imax)  <=> i0 <= ex1-2（循环里总成立）
+                    else if (i0 <= ex1 - 2) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                } else if (sfx < ZEO) {
+                    // Fortran: if(i-3 >= imin)
+                    // (iF-3) >= iminF  <=> (i0-2) >= iminF
+                    if ((i0 - 2) >= iminF) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                    // elseif(i-2 >= imin) <=> (i0-1) >= iminF
+                    else if ((i0 - 1) >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    }
+                    // elseif(i-1 >= imin) <=> i0 >= iminF
+                    else if (i0 >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                }
+
+                // ---------------- y direction ----------------
+                const double sfy = Sfy[p];
+                if (sfy > ZEO) {
+                    // jF+3 <= ex2 <=> j0+4 <= ex2 <=> j0 <= ex2-4
+                    if (j0 <= ex2 - 4) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    } else if (j0 <= ex2 - 3) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 <= ex2 - 2) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    }
+                } else if (sfy < ZEO) {
+                    if ((j0 - 2) >= jminF) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    } else if ((j0 - 1) >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    }
+                }
+
+                // ---------------- z direction ----------------
+                const double sfz = Sfz[p];
+                if (sfz > ZEO) {
+                    if (k0 <= ex3 - 4) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    } else if (k0 <= ex3 - 3) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 <= ex3 - 2) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    }
+                } else if (sfz < ZEO) {
+                    if ((k0 - 2) >= kminF) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    } else if ((k0 - 1) >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    }
+                }
+            }
+        }
+    }
+    free(fh);
+}
+
+
+
+
+

AMSS_NCKU_source/BSSN/lopsided_kodis_c.C

@@ -0,0 +1,248 @@
+#include "tool.h"
+
+/*
+ * Combined advection (lopsided) + KO dissipation (kodis).
+ * Uses one shared symmetry_bd buffer per call.
+ */
+void lopsided_kodis(const int ex[3],
+                    const double *X, const double *Y, const double *Z,
+                    const double *f, double *f_rhs,
+                    const double *Sfx, const double *Sfy, const double *Sfz,
+                    int Symmetry, const double SoA[3], double eps)
+{
+    const double ZEO = 0.0, ONE = 1.0, F3 = 3.0;
+    const double F6 = 6.0, F18 = 18.0;
+    const double F12 = 12.0, F10 = 10.0, EIT = 8.0;
+    const double SIX = 6.0, FIT = 15.0, TWT = 20.0;
+    const double cof = 64.0; // 2^6
+
+    const int NO_SYMM = 0, EQ_SYMM = 1;
+
+    const int ex1 = ex[0], ex2 = ex[1], ex3 = ex[2];
+
+    const double dX = X[1] - X[0];
+    const double dY = Y[1] - Y[0];
+    const double dZ = Z[1] - Z[0];
+
+    const double d12dx = ONE / F12 / dX;
+    const double d12dy = ONE / F12 / dY;
+    const double d12dz = ONE / F12 / dZ;
+
+    const int imaxF = ex1;
+    const int jmaxF = ex2;
+    const int kmaxF = ex3;
+
+    int iminF = 1, jminF = 1, kminF = 1;
+    if (Symmetry > NO_SYMM && fabs(Z[0]) < dZ) kminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(X[0]) < dX) iminF = -2;
+    if (Symmetry > EQ_SYMM && fabs(Y[0]) < dY) jminF = -2;
+
+    // fh for Fortran-style domain (-2:ex1,-2:ex2,-2:ex3)
+    const size_t nx = (size_t)ex1 + 3;
+    const size_t ny = (size_t)ex2 + 3;
+    const size_t nz = (size_t)ex3 + 3;
+    const size_t fh_size = nx * ny * nz;
+
+    double *fh = (double*)malloc(fh_size * sizeof(double));
+    if (!fh) return;
+
+    symmetry_bd(3, ex, f, fh, SoA);
+
+    // Advection (same stencil logic as lopsided_c.C)
+    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
+        const int kF = k0 + 1;
+        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
+            const int jF = j0 + 1;
+            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
+                const int iF = i0 + 1;
+                const size_t p = idx_ex(i0, j0, k0, ex);
+
+                const double sfx = Sfx[p];
+                if (sfx > ZEO) {
+                    if (i0 <= ex1 - 4) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    } else if (i0 <= ex1 - 3) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    } else if (i0 <= ex1 - 2) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    }
+                } else if (sfx < ZEO) {
+                    if ((i0 - 2) >= iminF) {
+                        f_rhs[p] -= sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF - 3, jF, kF, ex)]);
+                    } else if ((i0 - 1) >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            ( fh[idx_fh_F(iF - 2, jF, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -      fh[idx_fh_F(iF + 2, jF, kF, ex)]);
+                    } else if (i0 >= iminF) {
+                        f_rhs[p] += sfx * d12dx *
+                            (-F3  * fh[idx_fh_F(iF - 1, jF, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF    , jF, kF, ex)]
+                             +F18 * fh[idx_fh_F(iF + 1, jF, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF + 2, jF, kF, ex)]
+                             +      fh[idx_fh_F(iF + 3, jF, kF, ex)]);
+                    }
+                }
+
+                const double sfy = Sfy[p];
+                if (sfy > ZEO) {
+                    if (j0 <= ex2 - 4) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    } else if (j0 <= ex2 - 3) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 <= ex2 - 2) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    }
+                } else if (sfy < ZEO) {
+                    if ((j0 - 2) >= jminF) {
+                        f_rhs[p] -= sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF - 3, kF, ex)]);
+                    } else if ((j0 - 1) >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            ( fh[idx_fh_F(iF, jF - 2, kF, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -      fh[idx_fh_F(iF, jF + 2, kF, ex)]);
+                    } else if (j0 >= jminF) {
+                        f_rhs[p] += sfy * d12dy *
+                            (-F3  * fh[idx_fh_F(iF, jF - 1, kF, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF    , kF, ex)]
+                             +F18 * fh[idx_fh_F(iF, jF + 1, kF, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF + 2, kF, ex)]
+                             +      fh[idx_fh_F(iF, jF + 3, kF, ex)]);
+                    }
+                }
+
+                const double sfz = Sfz[p];
+                if (sfz > ZEO) {
+                    if (k0 <= ex3 - 4) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    } else if (k0 <= ex3 - 3) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 <= ex3 - 2) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    }
+                } else if (sfz < ZEO) {
+                    if ((k0 - 2) >= kminF) {
+                        f_rhs[p] -= sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF - 3, ex)]);
+                    } else if ((k0 - 1) >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            ( fh[idx_fh_F(iF, jF, kF - 2, ex)]
+                             -EIT * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             +EIT * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -      fh[idx_fh_F(iF, jF, kF + 2, ex)]);
+                    } else if (k0 >= kminF) {
+                        f_rhs[p] += sfz * d12dz *
+                            (-F3  * fh[idx_fh_F(iF, jF, kF - 1, ex)]
+                             -F10 * fh[idx_fh_F(iF, jF, kF    , ex)]
+                             +F18 * fh[idx_fh_F(iF, jF, kF + 1, ex)]
+                             -F6  * fh[idx_fh_F(iF, jF, kF + 2, ex)]
+                             +      fh[idx_fh_F(iF, jF, kF + 3, ex)]);
+                    }
+                }
+            }
+        }
+    }
+
+    // KO dissipation (same domain restriction as kodiss_c.C)
+    if (eps > ZEO) {
+        const int i0_lo = (iminF + 2 > 0) ? iminF + 2 : 0;
+        const int j0_lo = (jminF + 2 > 0) ? jminF + 2 : 0;
+        const int k0_lo = (kminF + 2 > 0) ? kminF + 2 : 0;
+        const int i0_hi = imaxF - 4; // inclusive
+        const int j0_hi = jmaxF - 4;
+        const int k0_hi = kmaxF - 4;
+
+        if (!(i0_lo > i0_hi || j0_lo > j0_hi || k0_lo > k0_hi)) {
+            for (int k0 = k0_lo; k0 <= k0_hi; ++k0) {
+                const int kF = k0 + 1;
+                for (int j0 = j0_lo; j0 <= j0_hi; ++j0) {
+                    const int jF = j0 + 1;
+                    for (int i0 = i0_lo; i0 <= i0_hi; ++i0) {
+                        const int iF = i0 + 1;
+                        const size_t p = idx_ex(i0, j0, k0, ex);
+
+                        const double Dx_term =
+                            ((fh[idx_fh_F(iF - 3, jF, kF, ex)] + fh[idx_fh_F(iF + 3, jF, kF, ex)]) -
+                             SIX * (fh[idx_fh_F(iF - 2, jF, kF, ex)] + fh[idx_fh_F(iF + 2, jF, kF, ex)]) +
+                             FIT * (fh[idx_fh_F(iF - 1, jF, kF, ex)] + fh[idx_fh_F(iF + 1, jF, kF, ex)]) -
+                             TWT *  fh[idx_fh_F(iF,     jF, kF, ex)]) / dX;
+
+                        const double Dy_term =
+                            ((fh[idx_fh_F(iF, jF - 3, kF, ex)] + fh[idx_fh_F(iF, jF + 3, kF, ex)]) -
+                             SIX * (fh[idx_fh_F(iF, jF - 2, kF, ex)] + fh[idx_fh_F(iF, jF + 2, kF, ex)]) +
+                             FIT * (fh[idx_fh_F(iF, jF - 1, kF, ex)] + fh[idx_fh_F(iF, jF + 1, kF, ex)]) -
+                             TWT *  fh[idx_fh_F(iF, jF,     kF, ex)]) / dY;
+
+                        const double Dz_term =
+                            ((fh[idx_fh_F(iF, jF, kF - 3, ex)] + fh[idx_fh_F(iF, jF, kF + 3, ex)]) -
+                             SIX * (fh[idx_fh_F(iF, jF, kF - 2, ex)] + fh[idx_fh_F(iF, jF, kF + 2, ex)]) +
+                             FIT * (fh[idx_fh_F(iF, jF, kF - 1, ex)] + fh[idx_fh_F(iF, jF, kF + 1, ex)]) -
+                             TWT *  fh[idx_fh_F(iF, jF, kF,     ex)]) / dZ;
+
+                        f_rhs[p] += (eps / cof) * (Dx_term + Dy_term + Dz_term);
+                    }
+                }
+            }
+        }
+    }
+
+    free(fh);
+}