[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_Input.py

@@ -16,9 +16,9 @@ import numpy
 File_directory   = "GW150914"                    ## output file directory
 Output_directory = "binary_output"               ## binary data file directory
                                                  ## The file directory name should not be too long
-MPI_processes    = 2                             ## number of mpi processes used in the simulation
+MPI_processes    = 64                             ## number of mpi processes used in the simulation
 
-GPU_Calculation  = "yes"                          ## Use GPU or not 
+GPU_Calculation  = "no"                          ## Use GPU or not 
                                                  ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
 CPU_Part         = 1.0
 GPU_Part         = 0.0
@@ -31,7 +31,7 @@ GPU_Part         = 0.0
 ## Setting the physical system and numerical method
 
 Symmetry                 = "equatorial-symmetry"   ## Symmetry of System: choose equatorial-symmetry、no-symmetry、octant-symmetry
-Equation_Class           = "BSSN-EScalar"                  ## Evolution Equation: choose "BSSN", "BSSN-EScalar", "BSSN-EM", "Z4C"
+Equation_Class           = "BSSN"                  ## Evolution Equation: choose "BSSN", "BSSN-EScalar", "BSSN-EM", "Z4C" 
                                                    ## If "BSSN-EScalar" is chosen, it is necessary to set other parameters below
 Initial_Data_Method      = "Ansorg-TwoPuncture"    ## initial data method: choose "Ansorg-TwoPuncture", "Lousto-Analytical", "Cao-Analytical", "KerrSchild-Analytical"
 Time_Evolution_Method    = "runge-kutta-45"        ## time evolution method: choose "runge-kutta-45"

AMSS_NCKU_Program.py

@@ -58,36 +58,31 @@ File_directory = os.path.join(input_data.File_directory)
 
 ## If the specified output directory exists, ask the user whether to continue
 if os.path.exists(File_directory):
-    auto_overwrite = str(getattr(input_data, "Auto_Overwrite_Output", "yes")).strip().lower()
+    print( " Output dictionary has been existed !!!  "                                                              )
-    if auto_overwrite in ("1", "yes", "y", "true", "on", "continue"):
+    print( " If you want to overwrite the existing file directory, please input 'continue' in the terminal !! "     ) 
-        print( " Output dictionary has been existed; Auto_Overwrite_Output=yes, continue the calculation. " )
+    print( " If you want to retain the existing file directory, please input 'stop' in the terminal to stop the "   ) 
-        print(                                                                                              )
+    print( " simulation. Then you can reset the output dictionary in the input script file AMSS_NCKU_Input.py !!! " )
-    else:
+    print(                                                                                                          )
-        print( " Output dictionary has been existed !!!  "                                                              )
+    ## Prompt whether to overwrite the existing directory
-        print( " If you want to overwrite the existing file directory, please input 'continue' in the terminal !! "     ) 
+    while True:
-        print( " If you want to retain the existing file directory, please input 'stop' in the terminal to stop the "   ) 
+        try:
-        print( " simulation. Then you can reset the output dictionary in the input script file AMSS_NCKU_Input.py !!! " )
+            inputvalue = input()
-        print(                                                                                                          )
+            ## If the user agrees to overwrite, proceed and remove the existing directory
-        ## Prompt whether to overwrite the existing directory
+            if ( inputvalue == "continue" ):
-        while True:
+                print( " Continue the calculation !!! " )
-            try:
+                print(                                  )
-                inputvalue = input()
+                break  
-                ## If the user agrees to overwrite, proceed and remove the existing directory
+            ## If the user chooses not to overwrite, exit and keep the existing directory
-                if ( inputvalue == "continue" ):
+            elif ( inputvalue == "stop" ):
-                    print( " Continue the calculation !!! " )
+                print( " Stop the calculation !!! "    )
-                    print(                                  )
+                sys.exit() 
-                    break  
+            ## If the user input is invalid, prompt again
-                ## If the user chooses not to overwrite, exit and keep the existing directory
+            else:
-                elif ( inputvalue == "stop" ):
-                    print( " Stop the calculation !!! "    )
-                    sys.exit() 
-                ## If the user input is invalid, prompt again
-                else:
-                    print( " Please input your choice !!! "                   )
-                    print( " Input 'continue' or 'stop' in the terminal !!! " )
-            except ValueError:
                 print( " Please input your choice !!! "                   )
                 print( " Input 'continue' or 'stop' in the terminal !!! " )
+        except ValueError:
+            print( " Please input your choice !!! "                   )
+            print( " Input 'continue' or 'stop' in the terminal !!! " )
         
 ## Remove the existing output directory if present
 shutil.rmtree(File_directory, ignore_errors=True)
@@ -131,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(                                                                                   )
@@ -263,7 +253,7 @@ print()
 if (input_data.GPU_Calculation == "no"):
     ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
 elif (input_data.GPU_Calculation == "yes"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE_CUDA")
+    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
 
 if not os.path.exists( ABE_file ):
     print(                                                                                                  )
@@ -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_source/BSSN/bssnEM_class.C

@@ -18,9 +18,6 @@ using namespace std;
 #include "Parallel.h"
 #include "bssnEM_class.h"
 #include "bssn_rhs.h"
-#if USE_CUDA_BSSN
-#include "bssn_rhs_cuda.h"
-#endif
 #include "empart.h"
 #include "initial_puncture.h"
 #include "initial_maxwell.h"
@@ -39,106 +36,6 @@ using namespace std;
 
 //================================================================================================
 
-#if USE_CUDA_BSSN
-namespace {
-
-bool fill_bssn_cuda_views_prefix(Block *cg, MyList<var> *vars,
-                                 double **host_views,
-                                 double *propspeeds = nullptr,
-                                 double *soa_flat = nullptr)
-{
-  int idx = 0;
-  while (vars && idx < BSSN_CUDA_STATE_COUNT)
-  {
-    host_views[idx] = cg->fgfs[vars->data->sgfn];
-    if (propspeeds)
-      propspeeds[idx] = vars->data->propspeed;
-    if (soa_flat)
-    {
-      soa_flat[3 * idx + 0] = vars->data->SoA[0];
-      soa_flat[3 * idx + 1] = vars->data->SoA[1];
-      soa_flat[3 * idx + 2] = vars->data->SoA[2];
-    }
-    vars = vars->next;
-    ++idx;
-  }
-  return idx == BSSN_CUDA_STATE_COUNT;
-}
-
-void skip_bssn_cuda_prefix(MyList<var> *&a, MyList<var> *&b, MyList<var> *&c)
-{
-  for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c; ++i)
-  {
-    a = a->next;
-    b = b->next;
-    c = c->next;
-  }
-}
-
-void skip_bssn_cuda_prefix(MyList<var> *&a, MyList<var> *&b,
-                           MyList<var> *&c, MyList<var> *&d)
-{
-  for (int i = 0; i < BSSN_CUDA_STATE_COUNT && a && b && c && d; ++i)
-  {
-    a = a->next;
-    b = b->next;
-    c = c->next;
-    d = d->next;
-  }
-}
-
-int run_bssn_em_cuda_substep(Block *cg,
-                             MyList<var> *state_in_list,
-                             MyList<var> *state_out_list,
-                             Patch *patch,
-                             double &dT_lev,
-                             double &TRK4,
-                             int &iter_count,
-                             int &Symmetry,
-                             int lev,
-                             double &ndeps,
-                             int &co,
-                             double &chitiny,
-                             var *rho, var *Sx, var *Sy, var *Sz,
-                             var *Sxx, var *Sxy, var *Sxz,
-                             var *Syy, var *Syz, var *Szz)
-{
-  double *state_in[BSSN_CUDA_STATE_COUNT];
-  double *state_out[BSSN_CUDA_STATE_COUNT];
-  double *matter[BSSN_CUDA_MATTER_COUNT] = {
-      cg->fgfs[rho->sgfn], cg->fgfs[Sx->sgfn], cg->fgfs[Sy->sgfn], cg->fgfs[Sz->sgfn],
-      cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn],
-      cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn]};
-  double propspeed[BSSN_CUDA_STATE_COUNT];
-  double soa_flat[3 * BSSN_CUDA_STATE_COUNT];
-  if (!fill_bssn_cuda_views_prefix(cg, state_in_list, state_in, propspeed, soa_flat) ||
-      !fill_bssn_cuda_views_prefix(cg, state_out_list, state_out))
-    return 1;
-
-  int apply_bam_bc = 0;
-#if (SommerType == 0)
-#ifndef WithShell
-  apply_bam_bc = (lev == 0) ? 1 : 0;
-#endif
-#endif
-  int use_zero_matter = 0;
-  int keep_resident_state = 0;
-  int apply_enforce_ga = 0;
-  return bssn_cuda_rk4_substep(cg,
-                               cg->shape, cg->X[0], cg->X[1], cg->X[2],
-                               state_in, state_out, matter,
-                               propspeed, soa_flat, patch->bbox,
-                               dT_lev, TRK4, iter_count, apply_bam_bc,
-                               Symmetry, lev, ndeps, co,
-                               use_zero_matter,
-                               keep_resident_state, apply_enforce_ga, chitiny);
-}
-
-}
-#endif
-
-//================================================================================================
-
 // Define bssnEM_class
 
 // It inherits some members and methods from the parent class bssn_class and modifies others.
@@ -347,8 +244,6 @@ void bssnEM_class::Initialize()
     CheckPoint->readcheck_sh(SH, myrank);
 #endif
 
-  Initialize_Level_Runtime();
-
   double h = GH->PatL[0]->data->blb->data->getdX(0);
   for (int i = 1; i < dim; i++)
     h = Mymin(h, GH->PatL[0]->data->blb->data->getdX(i));
@@ -958,7 +853,6 @@ void bssnEM_class::Step(int lev, int YN)
                      cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
 #endif
 
-        bool used_gpu_substep = false;
         if (
             f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2],
                                  cg->fgfs[phi0->sgfn],
@@ -980,16 +874,7 @@ void bssnEM_class::Step(int lev, int YN)
                                  cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], 
                                  cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn],
                                  Symmetry, lev, ndeps) ||
-#if USE_CUDA_BSSN
+            f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-            ((used_gpu_substep =
-                  (run_bssn_em_cuda_substep(cg, StateList, SynchList_pre, Pp->data,
-                                            dT_lev, TRK4, iter_count, Symmetry, lev,
-                                            ndeps, pre, chitiny,
-                                            rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0))
-                 ? 0
-                 : 1) ||
-#endif
-            (!used_gpu_substep && f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
                                cg->fgfs[phi0->sgfn], cg->fgfs[trK0->sgfn],
                                cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], 
                                cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
@@ -1022,7 +907,7 @@ void bssnEM_class::Step(int lev, int YN)
                                cg->fgfs[Cons_Ham->sgfn],
                                cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
                                cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                               Symmetry, lev, ndeps, pre)))
+                               Symmetry, lev, ndeps, pre))
         {
           cout << "find NaN in domain: (" 
                << cg->bbox[0] << ":" << cg->bbox[3] << "," 
@@ -1035,10 +920,6 @@ void bssnEM_class::Step(int lev, int YN)
         {
           MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; 
           // we do not check the correspondence here
-#if USE_CUDA_BSSN
-          if (used_gpu_substep)
-            skip_bssn_cuda_prefix(varl0, varl, varlrhs);
-#endif
           
           while (varl0)
           {
@@ -1340,7 +1221,7 @@ void bssnEM_class::Step(int lev, int YN)
   }
 #endif
 
-  Parallel::Sync_cached(GH->PatL[lev], SynchList_pre, Symmetry, sync_cache_pre[lev]);
+  Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
 
 #ifdef WithShell
   if (lev == 0)
@@ -1428,7 +1309,6 @@ void bssnEM_class::Step(int lev, int YN)
                          cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
 #endif
 
-          bool used_gpu_substep = false;
           if (
               f_compute_rhs_empart(cg->shape, cg->X[0], cg->X[1], cg->X[2],
                                    cg->fgfs[phi->sgfn],
@@ -1450,16 +1330,7 @@ void bssnEM_class::Step(int lev, int YN)
                                    cg->fgfs[Sxx->sgfn], cg->fgfs[Sxy->sgfn], cg->fgfs[Sxz->sgfn], 
                                    cg->fgfs[Syy->sgfn], cg->fgfs[Syz->sgfn], cg->fgfs[Szz->sgfn],
                                    Symmetry, lev, ndeps) ||
-#if USE_CUDA_BSSN
+              f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
-              ((used_gpu_substep =
-                    (run_bssn_em_cuda_substep(cg, SynchList_pre, SynchList_cor, Pp->data,
-                                              dT_lev, TRK4, iter_count, Symmetry, lev,
-                                              ndeps, cor, chitiny,
-                                              rho, Sx, Sy, Sz, Sxx, Sxy, Sxz, Syy, Syz, Szz) == 0))
-                   ? 0
-                   : 1) ||
-#endif
-              (!used_gpu_substep && f_compute_rhs_bssn(cg->shape, TRK4, cg->X[0], cg->X[1], cg->X[2],
                                  cg->fgfs[phi->sgfn], cg->fgfs[trK->sgfn],
                                  cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], 
                                  cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
@@ -1491,7 +1362,7 @@ void bssnEM_class::Step(int lev, int YN)
                                  cg->fgfs[Cons_Ham->sgfn],
                                  cg->fgfs[Cons_Px->sgfn], cg->fgfs[Cons_Py->sgfn], cg->fgfs[Cons_Pz->sgfn],
                                  cg->fgfs[Cons_Gx->sgfn], cg->fgfs[Cons_Gy->sgfn], cg->fgfs[Cons_Gz->sgfn],
-                                 Symmetry, lev, ndeps, cor)))
+                                 Symmetry, lev, ndeps, cor))
           {
             cout << "find NaN in domain: (" 
                  << cg->bbox[0] << ":" << cg->bbox[3] << "," 
@@ -1503,10 +1374,6 @@ void bssnEM_class::Step(int lev, int YN)
           {
             MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; 
             // we do not check the correspondence here
-#if USE_CUDA_BSSN
-            if (used_gpu_substep)
-              skip_bssn_cuda_prefix(varl0, varl, varl1, varlrhs);
-#endif
             
             while (varl0)
             {
@@ -1816,7 +1683,7 @@ void bssnEM_class::Step(int lev, int YN)
     }
 #endif
 
-    Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_cor[lev]);
+    Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
 
 #ifdef WithShell
     if (lev == 0)

AMSS_NCKU_source/BSSN/bssn_class.h

@@ -48,7 +48,6 @@ public:
        double StartTime, TotalTime;
        double AnasTime, DumpTime, d2DumpTime, CheckTime;
        double LastAnas, LastConsOut;
-       bool cuda_level0_constraint_cache_valid;
        int *ConstraintRefreshLevels;
        double Courant;
        double numepss, numepsb, numepsh;
@@ -184,9 +183,6 @@ public:
        virtual void Constraint_Out();
        virtual void Compute_Constraint();
 
-protected:
-       void Initialize_Level_Runtime();
-
 #ifdef With_AHF
 protected:
        MyList<var> *AHList, *AHDList, *GaugeList;

AMSS_NCKU_source/BSSN/bssn_rhs.h

@@ -6,7 +6,6 @@
 #define f_compute_rhs_bssn compute_rhs_bssn
 #define f_compute_rhs_bssn_ss compute_rhs_bssn_ss
 #define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar
-#define f_compute_rhs_bssn_escalar_matter compute_rhs_bssn_escalar_matter
 #define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss
 #define f_compute_rhs_Z4c compute_rhs_z4c
 #define f_compute_rhs_Z4cnot compute_rhs_z4cnot
@@ -17,7 +16,6 @@
 #define f_compute_rhs_bssn COMPUTE_RHS_BSSN
 #define f_compute_rhs_bssn_ss COMPUTE_RHS_BSSN_SS
 #define f_compute_rhs_bssn_escalar COMPUTE_RHS_BSSN_ESCALAR
-#define f_compute_rhs_bssn_escalar_matter COMPUTE_RHS_BSSN_ESCALAR_MATTER
 #define f_compute_rhs_bssn_escalar_ss COMPUTE_RHS_BSSN_ESCALAR_SS
 #define f_compute_rhs_Z4c COMPUTE_RHS_Z4C
 #define f_compute_rhs_Z4cnot COMPUTE_RHS_Z4CNOT
@@ -28,7 +26,6 @@
 #define f_compute_rhs_bssn compute_rhs_bssn_
 #define f_compute_rhs_bssn_ss compute_rhs_bssn_ss_
 #define f_compute_rhs_bssn_escalar compute_rhs_bssn_escalar_
-#define f_compute_rhs_bssn_escalar_matter compute_rhs_bssn_escalar_matter_
 #define f_compute_rhs_bssn_escalar_ss compute_rhs_bssn_escalar_ss_
 #define f_compute_rhs_Z4c compute_rhs_z4c_
 #define f_compute_rhs_Z4cnot compute_rhs_z4cnot_
@@ -99,20 +96,6 @@ extern "C"
                                   int &, int &, double &, int &, int &);
 }
 
-extern "C"
-{
-        int f_compute_rhs_bssn_escalar_matter(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 *,                                                                        // Sphi, Spi rhs
-                                              double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, // stress-energy
-                                              int &, int &, double &);
-}
-
 extern "C"
 {
         int f_compute_rhs_bssn_escalar(int *, double &, double *, double *, double *,                                                      // ex,T,X,Y,Z

AMSS_NCKU_source/BSSN/bssn_rhs_c.C

@@ -1098,12 +1098,12 @@ int f_compute_rhs_bssn(int *ex, double &T,
             betaz_rhs[i] = FF * dtSfz[i];
 
             reta[i] =
-                gupxx[i] * chix[i] * chix[i]
+                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
-                + gupyy[i] * chiy[i] * chiy[i]
+                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
-                + gupzz[i] * chiz[i] * chiz[i]
+                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
-                + TWO * ( gupxy[i] * chix[i] * chiy[i]
+                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
-                        + gupxz[i] * chix[i] * chiz[i]
+                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
-                        + gupyz[i] * chiy[i] * chiz[i] );
+                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
 
             #if (GAUGE == 2)
             reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
@@ -1116,12 +1116,12 @@ int f_compute_rhs_bssn(int *ex, double &T,
             dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
             #elif (GAUGE == 4 || GAUGE == 5)
             reta[i] =
-                gupxx[i] * chix[i] * chix[i]
+                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
-                + gupyy[i] * chiy[i] * chiy[i]
+                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
-                + gupzz[i] * chiz[i] * chiz[i]
+                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
-                + TWO * ( gupxy[i] * chix[i] * chiy[i]
+                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
-                        + gupxz[i] * chix[i] * chiz[i]
+                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
-                        + gupyz[i] * chiy[i] * chiz[i] );
+                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
 
             #if (GAUGE == 4)
             reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );

AMSS_NCKU_source/BSSN_GPU/bssn_gpu.h

@@ -0,0 +1,73 @@
+
+#ifndef BSSN_GPU_H_
+#define BSSN_GPU_H_
+#include "bssn_macro.h"
+#include "macrodef.fh"
+
+#define DEVICE_ID 0
+// #define DEVICE_ID_BY_MPI_RANK
+#define GRID_DIM 256
+#define BLOCK_DIM 128
+
+#define _FH2_(i, j, k) fh[(i) + (j) * _1D_SIZE[2] + (k) * _2D_SIZE[2]]
+#define _FH3_(i, j, k) fh[(i) + (j) * _1D_SIZE[3] + (k) * _2D_SIZE[3]]
+#define pow2(x) ((x) * (x))
+#define TimeBetween(a, b) ((b.tv_sec - a.tv_sec) + (b.tv_usec - a.tv_usec) / 1000000.0f)
+#define M_ metac.
+#define Mh_ meta->
+#define Ms_ metassc.
+#define Msh_ metass->
+
+// #define TIMING
+
+#define RHS_SS_PARA int calledby, int mpi_rank, int *ex, double &T, double *crho, double *sigma, double *R, double *X, double *Y, double *Z, double *drhodx, double *drhody, double *drhodz, double *dsigmadx, double *dsigmady, double *dsigmadz, double *dRdx, double *dRdy, double *dRdz, double *drhodxx, double *drhodxy, double *drhodxz, double *drhodyy, double *drhodyz, double *drhodzz, double *dsigmadxx, double *dsigmadxy, double *dsigmadxz, double *dsigmadyy, double *dsigmadyz, double *dsigmadzz, double *dRdxx, double *dRdxy, double *dRdxz, double *dRdyy, double *dRdyz, double *dRdzz, double *chi, double *trK, double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz, double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz, double *Gamx, double *Gamy, double *Gamz, double *Lap, double *betax, double *betay, double *betaz, double *dtSfx, double *dtSfy, double *dtSfz, double *chi_rhs, double *trK_rhs, double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs, double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs, double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs, double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs, double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs, double *rho, double *Sx, double *Sy, double *Sz, double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz, double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz, double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz, double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz, double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz, double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res, double *Gmx_Res, double *Gmy_Res, double *Gmz_Res, int &Symmetry, int &Lev, double &eps, int &sst, int &co
+
+/**  main function */
+int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,
+            double *X, double *Y, double *Z,
+
+            double *chi, double *trK,
+
+            double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
+
+            double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
+
+            double *Gamx, double *Gamy, double *Gamz,
+
+            double *Lap, double *betax, double *betay, double *betaz,
+
+            double *dtSfx, double *dtSfy, double *dtSfz,
+
+            double *chi_rhs, double *trK_rhs,
+
+            double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
+
+            double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
+
+            double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
+
+            double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
+
+            double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
+
+            double *rho, double *Sx, double *Sy, double *Sz, double *Sxx,
+            double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
+
+            double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
+
+            double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
+
+            double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
+
+            double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
+
+            double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res,
+            double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
+            int &Symmetry, int &Lev, double &eps, int &co);
+
+int gpu_rhs_ss(RHS_SS_PARA);
+
+/** Init GPU side data in GPUMeta. */
+// void init_fluid_meta_gpu(GPUMeta *gpu_meta);
+
+#endif

AMSS_NCKU_source/BSSN_GPU/bssn_gpu_class.h

@@ -0,0 +1,210 @@
+
+#ifndef BSSN_GPU_CLASS_H
+#define BSSN_GPU_CLASS_H
+
+#ifdef newc
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cstdlib>
+#include <string>
+#include <cmath>
+using namespace std;
+#else
+#include <iostream.h>
+#include <iomanip.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#endif
+
+#include <mpi.h>
+
+#include "macrodef.h"
+#include "cgh.h"
+#include "ShellPatch.h"
+#include "misc.h"
+#include "var.h"
+#include "MyList.h"
+#include "monitor.h"
+#include "surface_integral.h"
+#include "checkpoint.h"
+
+// added by yangquan
+#include "bssn_macro.h"
+
+extern void setpbh(int iBHN, double **iPBH, double *iMass, int rBHN);
+
+class bssn_class
+{
+public:
+       // added by yangquan
+       //----------------------
+       int gpu_num_mynode;
+       int cpu_core_num_mynode;
+       int mpi_process_num_mynode;
+       int my_sequence_mynode;
+       int mynode_id;
+       int use_gpu;
+
+       virtual void Step_GPU(int lev, int YN);
+       virtual void Get_runtime_envirment();
+       // virtual void Step_OPENMP(int lev,int YN);
+       //----------------------
+
+       int ngfs;
+       int nprocs, myrank;
+       cgh *GH;
+       ShellPatch *SH;
+       double PhysTime;
+
+       int checkrun;
+       char checkfilename[50];
+       int Steps;
+       double StartTime, TotalTime;
+       double AnasTime, DumpTime, d2DumpTime, CheckTime;
+       double LastAnas, LastConsOut;
+       double Courant;
+       double numepss, numepsb, numepsh;
+       int Symmetry;
+       int maxl, decn;
+       double maxrex, drex;
+       int trfls, a_lev;
+
+       double dT;
+       double chitiny;
+
+       double **Porg0, **Porgbr, **Porg, **Porg1, **Porg_rhs;
+       int BH_num, BH_num_input;
+       double *Mass, *Pmom, *Spin;
+       double ADMMass;
+
+       var *phio, *trKo;
+       var *gxxo, *gxyo, *gxzo, *gyyo, *gyzo, *gzzo;
+       var *Axxo, *Axyo, *Axzo, *Ayyo, *Ayzo, *Azzo;
+       var *Gmxo, *Gmyo, *Gmzo;
+       var *Lapo, *Sfxo, *Sfyo, *Sfzo;
+       var *dtSfxo, *dtSfyo, *dtSfzo;
+
+       var *phi0, *trK0;
+       var *gxx0, *gxy0, *gxz0, *gyy0, *gyz0, *gzz0;
+       var *Axx0, *Axy0, *Axz0, *Ayy0, *Ayz0, *Azz0;
+       var *Gmx0, *Gmy0, *Gmz0;
+       var *Lap0, *Sfx0, *Sfy0, *Sfz0;
+       var *dtSfx0, *dtSfy0, *dtSfz0;
+
+       var *phi, *trK;
+       var *gxx, *gxy, *gxz, *gyy, *gyz, *gzz;
+       var *Axx, *Axy, *Axz, *Ayy, *Ayz, *Azz;
+       var *Gmx, *Gmy, *Gmz;
+       var *Lap, *Sfx, *Sfy, *Sfz;
+       var *dtSfx, *dtSfy, *dtSfz;
+
+       var *phi1, *trK1;
+       var *gxx1, *gxy1, *gxz1, *gyy1, *gyz1, *gzz1;
+       var *Axx1, *Axy1, *Axz1, *Ayy1, *Ayz1, *Azz1;
+       var *Gmx1, *Gmy1, *Gmz1;
+       var *Lap1, *Sfx1, *Sfy1, *Sfz1;
+       var *dtSfx1, *dtSfy1, *dtSfz1;
+
+       var *phi_rhs, *trK_rhs;
+       var *gxx_rhs, *gxy_rhs, *gxz_rhs, *gyy_rhs, *gyz_rhs, *gzz_rhs;
+       var *Axx_rhs, *Axy_rhs, *Axz_rhs, *Ayy_rhs, *Ayz_rhs, *Azz_rhs;
+       var *Gmx_rhs, *Gmy_rhs, *Gmz_rhs;
+       var *Lap_rhs, *Sfx_rhs, *Sfy_rhs, *Sfz_rhs;
+       var *dtSfx_rhs, *dtSfy_rhs, *dtSfz_rhs;
+
+       var *rho, *Sx, *Sy, *Sz, *Sxx, *Sxy, *Sxz, *Syy, *Syz, *Szz;
+
+       var *Gamxxx, *Gamxxy, *Gamxxz, *Gamxyy, *Gamxyz, *Gamxzz;
+       var *Gamyxx, *Gamyxy, *Gamyxz, *Gamyyy, *Gamyyz, *Gamyzz;
+       var *Gamzxx, *Gamzxy, *Gamzxz, *Gamzyy, *Gamzyz, *Gamzzz;
+
+       var *Rxx, *Rxy, *Rxz, *Ryy, *Ryz, *Rzz;
+
+       var *Rpsi4, *Ipsi4;
+       var *t1Rpsi4, *t1Ipsi4, *t2Rpsi4, *t2Ipsi4;
+
+       var *Cons_Ham, *Cons_Px, *Cons_Py, *Cons_Pz, *Cons_Gx, *Cons_Gy, *Cons_Gz;
+
+#ifdef Point_Psi4
+       var *phix, *phiy, *phiz;
+       var *trKx, *trKy, *trKz;
+       var *Axxx, *Axxy, *Axxz;
+       var *Axyx, *Axyy, *Axyz;
+       var *Axzx, *Axzy, *Axzz;
+       var *Ayyx, *Ayyy, *Ayyz;
+       var *Ayzx, *Ayzy, *Ayzz;
+       var *Azzx, *Azzy, *Azzz;
+#endif
+       // FIXME: uc = StateList, up = OldStateList, upp = SynchList_cor; so never touch these three data
+       MyList<var> *StateList, *SynchList_pre, *SynchList_cor, *RHSList;
+       MyList<var> *OldStateList, *DumpList;
+       MyList<var> *ConstraintList;
+
+       monitor *ErrorMonitor, *Psi4Monitor, *BHMonitor, *MAPMonitor;
+       monitor *ConVMonitor;
+       surface_integral *Waveshell;
+       checkpoint *CheckPoint;
+
+public:
+       bssn_class(double Couranti, double StartTimei, double TotalTimei, double DumpTimei, double d2DumpTimei, double CheckTimei, double AnasTimei,
+                  int Symmetryi, int checkruni, char *checkfilenamei, double numepssi, double numepsbi, double numepshi,
+                  int a_levi, int maxli, int decni, double maxrexi, double drexi);
+       ~bssn_class();
+
+       void Evolve(int Steps);
+       void RecursiveStep(int lev);
+#if (PSTR == 1)
+       void ParallelStep();
+       void SHStep();
+#endif
+       void RestrictProlong(int lev, int YN, bool BB, MyList<var> *SL, MyList<var> *OL, MyList<var> *corL);
+       void RestrictProlong_aux(int lev, int YN, bool BB, MyList<var> *SL, MyList<var> *OL, MyList<var> *corL);
+       void RestrictProlong(int lev, int YN, bool BB);
+       void ProlongRestrict(int lev, int YN, bool BB);
+       void Setup_Black_Hole_position();
+       void compute_Porg_rhs(double **BH_PS, double **BH_RHS, var *forx, var *fory, var *forz, int lev);
+       bool read_Pablo_file(int *ext, double *datain, char *filename);
+       void write_Pablo_file(int *ext, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax,
+                             char *filename);
+       void AnalysisStuff(int lev, double dT_lev);
+       void Setup_KerrSchild();
+       void Enforce_algcon(int lev, int fg);
+
+       void testRestrict();
+       void testOutBd();
+
+       virtual void Setup_Initial_Data_Lousto();
+       virtual void Setup_Initial_Data_Cao();
+       virtual void Initialize();
+       virtual void Read_Ansorg();
+       virtual void Read_Pablo() {};
+       virtual void Compute_Psi4(int lev);
+       virtual void Step(int lev, int YN);
+       virtual void Interp_Constraint(bool infg);
+       virtual void Constraint_Out();
+       virtual void Compute_Constraint();
+
+#ifdef With_AHF
+protected:
+       MyList<var> *AHList, *AHDList, *GaugeList;
+       int AHfindevery;
+       double AHdumptime;
+       int *lastahdumpid, HN_num; // number of possible horizons
+       int *findeveryl;
+       double *xc, *yc, *zc, *xr, *yr, *zr;
+       bool *trigger;
+       double *dTT;
+       int *dumpid;
+
+public:
+       void AH_Prepare_derivatives();
+       bool AH_Interp_Points(MyList<var> *VarList,
+                             int NN, double **XX,
+                             double *Shellf, int Symmetryi);
+       void AH_Step_Find(int lev, double dT_lev);
+#endif
+};
+#endif /* BSSN_GPU_CLASS_H */