feat: port GPU code to CUDA 13 and enable GPU computation

Major changes:
   - Update makefile.inc for CUDA 13.1 with sm_89 architecture (RTX 4050)
   - Replace deprecated cudaThreadSynchronize() with cudaDeviceSynchronize()
   - Add CUDA_SAFE_CALL macro for CUDA 13 compatibility
   - Fix duplicate function definitions (compare_result_gpu, SHStep)
   - Fix syntax error in bssn_step_gpu.C
   - Enable GPU calculation in AMSS_NCKU_Input.py
   - Successfully build ABEGPU executable

.gitignore

@@ -1,5 +1,3 @@
 __pycache__
 GW150914
-GW150914*
+
-.codex
-docs/

AMSS_NCKU_Input.py

@@ -16,12 +16,12 @@ 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    = 64                             ## number of mpi processes used in the simulation
+MPI_processes    = 96                             ## number of mpi processes used in the simulation
 
-GPU_Calculation  = "no"                          ## Use GPU or not 
+GPU_Calculation  = "yes"                         ## Use GPU or not
-                                                 ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
+                                                 ## GPU support has been updated for CUDA 13
-CPU_Part         = 1.0
+CPU_Part         = 0.0
-GPU_Part         = 0.0
+GPU_Part         = 1.0
 
 #################################################
 

AMSS_NCKU_Program.py

@@ -9,19 +9,9 @@
 ##################################################################
 
 
-##################################################################
+##################################################################
-
+
-## Guard against re-execution by multiprocessing child processes.
+## Print program introduction
-## Without this, using 'spawn' or 'forkserver' context would cause every
-## worker to re-run the entire script.
-if __name__ != '__main__':
-    import sys as _sys
-    _sys.exit(0)
-
-
-##################################################################
-
-## Print program introduction
 
 import print_information
 
@@ -432,36 +422,31 @@ print( " Plotting the txt and binary results data from the AMSS-NCKU simulation
 print(                                                                          )
 
 
-import plot_xiaoqu
+import plot_xiaoqu
-import plot_GW_strain_amplitude_xiaoqu
+import plot_GW_strain_amplitude_xiaoqu
-from parallel_plot_helper import run_plot_tasks_parallel
+
-
+## Plot black hole trajectory
-plot_tasks = []
+plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
-
+plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
-## Plot black hole trajectory
+
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory) ) )
+## Plot black hole separation vs. time
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory) ) )
+plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
-
+
-## Plot black hole separation vs. time
+## Plot gravitational waveforms (psi4 and strain amplitude)
-plot_tasks.append( ( plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory) ) )
+for i in range(input_data.Detector_Number):
-
+    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
-## Plot gravitational waveforms (psi4 and strain amplitude)
+    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
-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 ADM mass evolution
-    plot_tasks.append( ( plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i) ) )
+for i in range(input_data.Detector_Number):
-
+    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
-## Plot ADM mass evolution
+
-for i in range(input_data.Detector_Number):
+## Plot Hamiltonian constraint violation over time
-    plot_tasks.append( ( plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i) ) )
+for i in range(input_data.grid_level):
-
+    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
-## Plot Hamiltonian constraint violation over time
+
-for i in range(input_data.grid_level):
+## Plot stored binary data
-    plot_tasks.append( ( plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i) ) )
+plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
-
-run_plot_tasks_parallel(plot_tasks)
-
-## Plot stored binary data
-plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
 
 print(                                                 )
 print( f" This Program Cost = {elapsed_time} Seconds " )

AMSS_NCKU_source/TwoPunctures.h

@@ -1,8 +1,7 @@
+
 #ifndef TWO_PUNCTURES_H
 #define TWO_PUNCTURES_H
 
-#include <omp.h>
-
 #define StencilSize 19
 #define N_PlaneRelax 1
 #define NRELAX 200
@@ -33,7 +32,7 @@ private:
        int npoints_A, npoints_B, npoints_phi;
 
        double target_M_plus, target_M_minus;
-
+       
        double admMass;
 
        double adm_tol;
@@ -43,18 +42,6 @@ private:
 
        int ntotal;
 
-       // ===== Precomputed spectral derivative matrices =====
-       double *D1_A, *D2_A;
-       double *D1_B, *D2_B;
-       double *DF1_phi, *DF2_phi;
-
-       // ===== Pre-allocated workspace for LineRelax (per-thread) =====
-       int max_threads;
-       double **ws_diag_be, **ws_e_be, **ws_f_be, **ws_b_be, **ws_x_be;
-       double **ws_l_be, **ws_u_be, **ws_d_be, **ws_y_be;
-       double **ws_diag_al, **ws_e_al, **ws_f_al, **ws_b_al, **ws_x_al;
-       double **ws_l_al, **ws_u_al, **ws_d_al, **ws_y_al;
-
        struct parameters
        {
               int nvar, n1, n2, n3;
@@ -71,28 +58,6 @@ public:
                     int Newtonmaxit);
        ~TwoPunctures();
 
-       // 02/07: New/modified methods
-       void allocate_workspace();
-       void free_workspace();
-       void precompute_derivative_matrices();
-       void build_cheb_deriv_matrices(int n, double *D1, double *D2);
-       void build_fourier_deriv_matrices(int N, double *DF1, double *DF2);
-       void Derivatives_AB3_MatMul(int nvar, int n1, int n2, int n3, derivs v);
-       void ThomasAlgorithm_ws(int N, double *b, double *a, double *c, double *x, double *q,
-                                double *l, double *u_ws, double *d, double *y);
-       void LineRelax_be_omp(double *dv,
-                         int const i, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols, int **cols,
-                         double **JFD, int tid);
-       void LineRelax_al_omp(double *dv,
-                         int const j, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols,
-                         int **cols, double **JFD, int tid);
-       void relax_omp(double *dv, int const nvar, int const n1, int const n2, int const n3,
-                  double const *rhs, int const *ncols, int **cols, double **JFD);
-
        void Solve();
        void set_initial_guess(derivs v);
        int index(int i, int j, int k, int l, int a, int b, int c, int d);
@@ -151,11 +116,23 @@ public:
        double BY_KKofxyz(double x, double y, double z);
        void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix);
        void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u);
+       void relax(double *dv, int const nvar, int const n1, int const n2, int const n3,
+                  double const *rhs, int const *ncols, int **cols, double **JFD);
+       void LineRelax_be(double *dv,
+                         int const i, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols, int **cols,
+                         double **JFD);
        void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
                          int n3, derivs dv, derivs u, double *values);
        void LinEquations(double A, double B, double X, double R,
                          double x, double r, double phi,
                          double y, double z, derivs dU, derivs U, double *values);
+       void LineRelax_al(double *dv,
+                         int const j, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols,
+                         int **cols, double **JFD);
        void ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q);
        void Save(char *fname);
        // provided by Vasileios Paschalidis (vpaschal@illinois.edu)
@@ -164,4 +141,4 @@ public:
        void SpecCoef(parameters par, int ivar, double *v, double *cf);
 };
 
-#endif /* TWO_PUNCTURES_H */
+#endif /* TWO_PUNCTURES_H */

AMSS_NCKU_source/bssn_gpu.cu

@@ -18,7 +18,7 @@ using namespace std;
 #include <fstream>
 #endif
 
-void compare_result_gpu(int ftag1,double * datac,int data_num){
+static void compare_result_gpu(int ftag1,double * datac,int data_num){
 	double * data = (double*)malloc(sizeof(double)*data_num);
 	cudaMemcpy(data, datac, data_num * sizeof(double), cudaMemcpyDeviceToHost);
 	compare_result(ftag1,data,data_num);
@@ -83,7 +83,7 @@ inline void sub_enforce_ga(int matrix_size){
 	double * trA = M_ chin1;
 	enforce_ga<<<GRID_DIM,BLOCK_DIM>>>(trA);
 	cudaMemset(trA,0,matrix_size * sizeof(double));
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 	
 	//cudaMemset(Mh_ gupxx,0,matrix_size * sizeof(double));
 	//trA gxx,gyy,gzz gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
@@ -273,13 +273,13 @@ __global__ void sub_symmetry_bd_partK(int ord,double * func, double * funcc,doub
 #endif //ifdef Vertex
 inline void sub_symmetry_bd(int ord,double * func, double * funcc,double * SoA){
 	sub_symmetry_bd_partF<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_symmetry_bd_partI<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc,SoA[0]);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_symmetry_bd_partJ<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc,SoA[1]);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_symmetry_bd_partK<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc,SoA[2]);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 }
 
 
@@ -378,9 +378,9 @@ inline void sub_fdderivs(double * f,double *fh,double *fxx,double *fxy,double *f
 	cudaMemset(fyy,0,_3D_SIZE[0] * sizeof(double));
 	cudaMemset(fyz,0,_3D_SIZE[0] * sizeof(double));
 	cudaMemset(fzz,0,_3D_SIZE[0] * sizeof(double));
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 	sub_fdderivs_part1<<<GRID_DIM,BLOCK_DIM>>>(f,fh,fxx,fxy,fxz,fyy,fyz,fzz);
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 }
 
 __global__ void sub_fderivs_part1(double * f,double * fh,double *fx,double *fy,double *fz  )
@@ -445,9 +445,9 @@ inline void sub_fderivs(double * f,double * fh,double *fx,double *fy,double *fz,
 	cudaMemset(fy,0,_3D_SIZE[0] * sizeof(double));
 	cudaMemset(fz,0,_3D_SIZE[0] * sizeof(double));
 	
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 	sub_fderivs_part1<<<GRID_DIM,BLOCK_DIM>>>(f,fh,fx,fy,fz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 }
 
 __global__ void computeRicci_part1(double * dst)
@@ -465,9 +465,9 @@ __global__ void computeRicci_part1(double * dst)
  inline void computeRicci(double * src,double* dst,double * SoA, Meta* meta)
 {
 	sub_fdderivs(src,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,SoA);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	computeRicci_part1<<<GRID_DIM,BLOCK_DIM>>>(dst);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 }/*Exception*/
 
@@ -524,9 +524,9 @@ __global__ void sub_kodis_part1(double *f,double *fh,double *f_rhs)
 inline void sub_kodis(double *f,double *fh,double *f_rhs,double *SoA)
 {
 	sub_symmetry_bd(3,f,fh,SoA);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_kodis_part1<<<GRID_DIM,BLOCK_DIM>>>(f,fh,f_rhs);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 }
  
 __global__ void  sub_lopsided_part1(double *f,double* fh,double *f_rhs,double *Sfx,double *Sfy,double *Sfz)
@@ -617,9 +617,9 @@ __global__ void  sub_lopsided_part1(double *f,double* fh,double *f_rhs,double *S
 
 inline void  sub_lopsided(double *f,double*fh,double *f_rhs,double *Sfx,double *Sfy,double *Sfz,double *SoA){
 	sub_symmetry_bd(3,f,fh,SoA);
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 	sub_lopsided_part1<<<GRID_DIM,BLOCK_DIM>>>(f,fh,f_rhs,Sfx,Sfy,Sfz);
-	cudaThreadSynchronize(); 
+	cudaDeviceSynchronize(); 
 }
 
 __global__ void compute_rhs_bssn_part1() 
@@ -2656,13 +2656,13 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 
 
 #ifdef TIMING1
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	gettimeofday(&tv2, NULL);
    	cout<<"TIME USED"<<TimeBetween(tv1, tv2)<<endl; 
 #endif	
 	//cout<<"GPU meta data ready.\n";
 	
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 
 //--------------test constant memory address & value--------------
 /*	double rank = mpi_rank;
@@ -2685,7 +2685,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 	//sub_enforce_ga(matrix_size);
 	//4.1-----compute rhs---------
 	compute_rhs_bssn_part1<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 
 	sub_fderivs(Mh_ betax,Mh_ fh,Mh_ betaxx,Mh_ betaxy,Mh_ betaxz,ass);
 	sub_fderivs(Mh_ betay,Mh_ fh,Mh_ betayx,Mh_ betayy,Mh_ betayz,sas);
@@ -2701,7 +2701,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 	sub_fderivs(Mh_ gyz,Mh_ fh,Mh_ gyzx,Mh_ gyzy,Mh_ gyzz, saa);
   	
   	compute_rhs_bssn_part2<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs(Mh_ betax,Mh_ fh,Mh_ gxxx,Mh_ gxyx,Mh_ gxzx,Mh_ gyyx,Mh_ gyzx,Mh_ gzzx,ass);
 	sub_fdderivs(Mh_ betay,Mh_ fh,Mh_ gxxy,Mh_ gxyy,Mh_ gxzy,Mh_ gyyy,Mh_ gyzy,Mh_ gzzy,sas);
@@ -2711,7 +2711,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 	sub_fderivs( Mh_ Gamz, Mh_ fh,Mh_ Gamzx, Mh_ Gamzy, Mh_ Gamzz,ssa);
 	
 	compute_rhs_bssn_part3<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	computeRicci(Mh_ dxx,Mh_ Rxx,sss, meta);
 	computeRicci(Mh_ dyy,Mh_ Ryy,sss, meta);
@@ -2720,20 +2720,20 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 	computeRicci(Mh_ gxz,Mh_ Rxz,asa, meta);
 	computeRicci(Mh_ gyz,Mh_ Ryz,saa, meta);
 	
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	compute_rhs_bssn_part4<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs(Mh_ chi,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,sss);
 	
 	compute_rhs_bssn_part5<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs(Mh_ Lap,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,sss);
 	
 	compute_rhs_bssn_part6<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
 	sub_fderivs(Mh_ chi,Mh_ fh, Mh_ dtSfx_rhs, Mh_ dtSfy_rhs, Mh_ dtSfz_rhs,sss);
@@ -2805,7 +2805,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 	
 	if(co == 0){
 		compute_rhs_bssn_part7<<<GRID_DIM,BLOCK_DIM>>>();
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 
 		sub_fderivs(Mh_ Axx,Mh_ fh,Mh_ gxxx,Mh_ gxxy,Mh_ gxxz,sss);
 		sub_fderivs(Mh_ Axy,Mh_ fh,Mh_ gxyx,Mh_ gxyy,Mh_ gxyz,aas);
@@ -2814,7 +2814,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 		sub_fderivs(Mh_ Ayz,Mh_ fh,Mh_ gyzx,Mh_ gyzy,Mh_ gyzz,saa);
 		sub_fderivs(Mh_ Azz,Mh_ fh,Mh_ gzzx,Mh_ gzzy,Mh_ gzzz,sss);
 		compute_rhs_bssn_part8<<<GRID_DIM,BLOCK_DIM>>>();
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 	}
 
 #if (ABV == 1)
@@ -2895,7 +2895,7 @@ int gpu_rhs(int calledby, int mpi_rank, int *ex, double &T,double *X, double *Y,
 //-------------------FOR GPU TEST----------------------
 //-----------------------------------------------------
 #ifdef TIMING
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	gettimeofday(&tv2, NULL);
    	cout<<"MPI rank is: "<<mpi_rank<<" GPU TIME is"<<TimeBetween(tv1, tv2)<<" (s)."<<endl; 
 #endif

AMSS_NCKU_source/bssn_gpu.h

@@ -4,6 +4,17 @@
 #include "bssn_macro.h"
 #include "macrodef.fh"
 
+// CUDA error checking macro for CUDA 13 compatibility
+#define CUDA_SAFE_CALL(call) \
+    do { \
+        cudaError_t err = call; \
+        if (err != cudaSuccess) { \
+            fprintf(stderr, "CUDA error in %s:%d: %s\n", __FILE__, __LINE__, \
+                    cudaGetErrorString(err)); \
+            exit(EXIT_FAILURE); \
+        } \
+    } while(0)
+
 #define DEVICE_ID 0
 // #define DEVICE_ID_BY_MPI_RANK
 #define GRID_DIM 256

AMSS_NCKU_source/bssn_gpu_rhs_ss.cu

@@ -20,7 +20,7 @@ using namespace std;
 
 __device__ volatile unsigned int global_count = 0;
 
-void compare_result_gpu(int ftag1,double * datac,int data_num){
+static void compare_result_gpu(int ftag1,double * datac,int data_num){
 	double * data = (double*)malloc(sizeof(double)*data_num);
 	cudaMemcpy(data, datac, data_num * sizeof(double), cudaMemcpyDeviceToHost);
 	compare_result(ftag1,data,data_num);
@@ -153,11 +153,11 @@ __global__ void sub_symmetry_bd_ss_partJ(int ord,double * func, double * funcc,d
 
 inline void sub_symmetry_bd_ss(int ord,double * func, double * funcc,double * SoA){
 	sub_symmetry_bd_ss_partF<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_symmetry_bd_ss_partI<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc,SoA[0]);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	sub_symmetry_bd_ss_partJ<<<GRID_DIM,BLOCK_DIM>>>(ord,func,funcc,SoA[1]);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 }
 
 __global__ void sub_fderivs_shc_part1(double *fx,double *fy,double *fz){
@@ -247,13 +247,13 @@ inline void sub_fderivs_shc(int& sst,double * f,double * fh,double *fx,double *f
 	//cudaMemset(Msh_ gy,0,h_3D_SIZE[0] * sizeof(double));
 	//cudaMemset(Msh_ gz,0,h_3D_SIZE[0] * sizeof(double));
 	sub_symmetry_bd_ss(2,f,fh,SoA1);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(0,fh,h_3D_SIZE[2]);	
 	sub_fderivs_sh<<<GRID_DIM,BLOCK_DIM>>>(fh,Msh_ gx,Msh_ gy,Msh_ gz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fderivs_shc_part1<<<GRID_DIM,BLOCK_DIM>>>(fx,fy,fz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(1,fx,h_3D_SIZE[0]);
 			//compare_result_gpu(2,fy,h_3D_SIZE[0]);
 			//compare_result_gpu(3,fz,h_3D_SIZE[0]);
@@ -451,17 +451,17 @@ inline void sub_fdderivs_shc(int& sst,double * f,double * fh,
 	
 	//fderivs_sh
 	sub_symmetry_bd_ss(2,f,fh,SoA1);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(1,fh,h_3D_SIZE[2]);	
 	sub_fderivs_sh<<<GRID_DIM,BLOCK_DIM>>>(fh,Msh_ gx,Msh_ gy,Msh_ gz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	//fdderivs_sh
 	sub_symmetry_bd_ss(2,f,fh,SoA1);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(21,fh,h_3D_SIZE[2]);
 	sub_fdderivs_sh<<<GRID_DIM,BLOCK_DIM>>>(fh,Msh_ gxx,Msh_ gxy,Msh_ gxz,Msh_ gyy,Msh_ gyz,Msh_ gzz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			/*compare_result_gpu(11,Msh_ gx,h_3D_SIZE[0]);
 			compare_result_gpu(12,Msh_ gy,h_3D_SIZE[0]);
 			compare_result_gpu(13,Msh_ gz,h_3D_SIZE[0]);
@@ -472,7 +472,7 @@ inline void sub_fdderivs_shc(int& sst,double * f,double * fh,
 			compare_result_gpu(5,Msh_ gyz,h_3D_SIZE[0]);
 			compare_result_gpu(6,Msh_ gzz,h_3D_SIZE[0]);*/
 	sub_fdderivs_shc_part1<<<GRID_DIM,BLOCK_DIM>>>(fxx,fxy,fxz,fyy,fyz,fzz);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			/*compare_result_gpu(1,fxx,h_3D_SIZE[0]);
 			compare_result_gpu(2,fxy,h_3D_SIZE[0]);
 			compare_result_gpu(3,fxz,h_3D_SIZE[0]);
@@ -496,9 +496,9 @@ __global__ void computeRicci_ss_part1(double * dst)
  inline void computeRicci_ss(int &sst,double * src,double* dst,double * SoA, Meta* meta)
 {
 	sub_fdderivs_shc(sst,src,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,SoA);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	computeRicci_ss_part1<<<GRID_DIM,BLOCK_DIM>>>(dst);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 }
 __global__ void sub_lopsided_ss_part1(double * dst)
@@ -516,9 +516,9 @@ __global__ void sub_lopsided_ss_part1(double * dst)
 inline void sub_lopsided_ss(int& sst,double *src,double* dst,double *SoA)
 {
 		sub_fderivs_shc(sst,src,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,SoA);
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 		sub_lopsided_ss_part1<<<GRID_DIM,BLOCK_DIM>>>(dst);
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 }
 
 __global__ void sub_kodis_sh_part1(double *f,double *fh,double *f_rhs)
@@ -590,11 +590,11 @@ inline void sub_kodis_ss(int &sst,double *f,double *fh,double *f_rhs,double *SoA
 	}
 			//compare_result_gpu(10,f,h_3D_SIZE[0]);
 	sub_symmetry_bd_ss(3,f,fh,SoA1);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(0,fh,h_3D_SIZE[3]);
 			
 	sub_kodis_sh_part1<<<GRID_DIM,BLOCK_DIM>>>(f,fh,f_rhs);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 			//compare_result_gpu(1,f_rhs,h_3D_SIZE[0]);
 }
 
@@ -2287,13 +2287,13 @@ int gpu_rhs_ss(RHS_SS_PARA)
 
 
 #ifdef TIMING1
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	gettimeofday(&tv2, NULL);
    	cout<<"TIME USED"<<TimeBetween(tv1, tv2)<<endl; 
 #endif	
 	//cout<<"GPU meta data ready.\n";
 	
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 
 
 //-------------get device info-------------------------------------
@@ -2306,7 +2306,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	//sub_enforce_ga(matrix_size);
 	//4.1-----compute rhs---------
 	compute_rhs_ss_part1<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 
 	sub_fderivs_shc(sst,Mh_ betax,Mh_ fh,Mh_ betaxx,Mh_ betaxy,Mh_ betaxz,ass);
 	sub_fderivs_shc(sst,Mh_ betay,Mh_ fh,Mh_ betayx,Mh_ betayy,Mh_ betayz,sas);
@@ -2322,7 +2322,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	sub_fderivs_shc(sst,Mh_ gyz,Mh_ fh,Mh_ gyzx,Mh_ gyzy,Mh_ gyzz, saa);
 	
 	compute_rhs_ss_part2<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs_shc(sst,Mh_ betax,Mh_ fh,Mh_ gxxx,Mh_ gxyx,Mh_ gxzx,Mh_ gyyx,Mh_ gyzx,Mh_ gzzx,ass);
 	sub_fdderivs_shc(sst,Mh_ betay,Mh_ fh,Mh_ gxxy,Mh_ gxyy,Mh_ gxzy,Mh_ gyyy,Mh_ gyzy,Mh_ gzzy,sas);
@@ -2332,7 +2332,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	sub_fderivs_shc( sst,Mh_ Gamz, Mh_ fh,Mh_ Gamzx, Mh_ Gamzy, Mh_ Gamzz,ssa);
 	
 	compute_rhs_ss_part3<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	computeRicci_ss(sst,Mh_ dxx,Mh_ Rxx,sss, meta);
 	computeRicci_ss(sst,Mh_ dyy,Mh_ Ryy,sss, meta);
@@ -2340,25 +2340,25 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	computeRicci_ss(sst,Mh_ gxy,Mh_ Rxy,aas, meta);
 	computeRicci_ss(sst,Mh_ gxz,Mh_ Rxz,asa, meta);
 	computeRicci_ss(sst,Mh_ gyz,Mh_ Ryz,saa, meta);
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	compute_rhs_ss_part4<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs_shc(sst,Mh_ chi,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,sss);
 	
-	//cudaThreadSynchronize();
+	//cudaDeviceSynchronize();
 	//compare_result_gpu(0,Mh_ chi,h_3D_SIZE[0]);
 	//compare_result_gpu(1,Mh_ chi,h_3D_SIZE[0]);
 	//compare_result_gpu(2,Mh_ fyz,h_3D_SIZE[0]);
 	
 	compute_rhs_ss_part5<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 	sub_fdderivs_shc(sst,Mh_ Lap,Mh_ fh,Mh_ fxx,Mh_ fxy,Mh_ fxz,Mh_ fyy,Mh_ fyz,Mh_ fzz,sss);
 	
 	compute_rhs_ss_part6<<<GRID_DIM,BLOCK_DIM>>>();
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	
 #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
 	sub_fderivs_shc(sst,Mh_ chi,Mh_ fh, Mh_ dtSfx_rhs, Mh_ dtSfy_rhs, Mh_ dtSfz_rhs,sss);
@@ -2423,7 +2423,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	}
 	if(co == 0){
 		compute_rhs_ss_part7<<<GRID_DIM,BLOCK_DIM>>>();
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 
 		sub_fderivs_shc(sst,Mh_ Axx,Mh_ fh,Mh_ gxxx,Mh_ gxxy,Mh_ gxxz,sss);
 		sub_fderivs_shc(sst,Mh_ Axy,Mh_ fh,Mh_ gxyx,Mh_ gxyy,Mh_ gxyz,aas);
@@ -2432,7 +2432,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 		sub_fderivs_shc(sst,Mh_ Ayz,Mh_ fh,Mh_ gyzx,Mh_ gyzy,Mh_ gyzz,saa);
 		sub_fderivs_shc(sst,Mh_ Azz,Mh_ fh,Mh_ gzzx,Mh_ gzzy,Mh_ gzzz,sss);
 		compute_rhs_ss_part8<<<GRID_DIM,BLOCK_DIM>>>();
-		cudaThreadSynchronize();
+		cudaDeviceSynchronize();
 	}
 	
 #if (ABV == 1)
@@ -2512,7 +2512,7 @@ int gpu_rhs_ss(RHS_SS_PARA)
 	//test kodis
 	//sub_kodis_sh(sst,Msh_ drhodx,Mh_ fh2,Msh_ drhody,sss);
 #ifdef TIMING
-	cudaThreadSynchronize();
+	cudaDeviceSynchronize();
 	gettimeofday(&tv2, NULL);
    	cout<<"MPI rank is: "<<mpi_rank<<" GPU TIME is"<<TimeBetween(tv1, tv2)<<" (s)."<<endl; 
 #endif

AMSS_NCKU_source/bssn_step_gpu.C

@@ -1676,8 +1676,11 @@ void bssn_class::Step_GPU(int lev, int YN)
 #endif // PSTR == ?
 
 //--------------------------With Shell--------------------------
+// Note: SHStep() implementation is in bssn_gpu_class.C
 
 #ifdef WithShell
+#if 0
+// This SHStep() implementation has been moved to bssn_gpu_class.C to avoid duplicate definition
 void bssn_class::SHStep()
 {
   int lev = 0;
@@ -1938,5 +1941,5 @@ void bssn_class::SHStep()
     sPp = sPp->next;
   }
 }
-d
+#endif // #if 0
 #endif // withshell

AMSS_NCKU_source/enforce_algebra.f90

@@ -17,62 +17,50 @@
   real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Axx,Axy,Axz
   real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz
 
-!~~~~~~~> Local variable:
+!~~~~~~~> Local variable:
-
+  
-  integer :: i,j,k
+  real*8, dimension(ex(1),ex(2),ex(3)) :: trA,detg
-  real*8 :: lgxx,lgyy,lgzz,ldetg
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
-  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
-  real*8 :: ltrA,lscale
+  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
-  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
+
-
+!~~~~~~>
-!~~~~~~>
+
-
+  gxx = dxx + ONE
-  do k=1,ex(3)
+  gyy = dyy + ONE
-  do j=1,ex(2)
+  gzz = dzz + ONE
-  do i=1,ex(1)
+
-
+  detg =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-    lgxx = dxx(i,j,k) + ONE
+          gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
-    lgyy = dyy(i,j,k) + ONE
+  gupxx =   ( gyy * gzz - gyz * gyz ) / detg
-    lgzz = dzz(i,j,k) + ONE
+  gupxy = - ( gxy * gzz - gyz * gxz ) / detg
-
+  gupxz =   ( gxy * gyz - gyy * gxz ) / detg
-    ldetg =  lgxx * lgyy * lgzz &
+  gupyy =   ( gxx * gzz - gxz * gxz ) / detg
-           + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) &
+  gupyz = - ( gxx * gyz - gxy * gxz ) / detg
-           + gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) &
+  gupzz =   ( gxx * gyy - gxy * gxy ) / detg
-           - gxz(i,j,k) * lgyy * gxz(i,j,k) &
+
-           - gxy(i,j,k) * gxy(i,j,k) * lgzz &
+  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
-           - lgxx * gyz(i,j,k) * gyz(i,j,k)
+       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
-
+
-    lgupxx =   ( lgyy * lgzz - gyz(i,j,k) * gyz(i,j,k) ) / ldetg
+  Axx = Axx - F1o3 * gxx * trA
-    lgupxy = - ( gxy(i,j,k) * lgzz - gyz(i,j,k) * gxz(i,j,k) ) / ldetg
+  Axy = Axy - F1o3 * gxy * trA
-    lgupxz =   ( gxy(i,j,k) * gyz(i,j,k) - lgyy * gxz(i,j,k) ) / ldetg
+  Axz = Axz - F1o3 * gxz * trA
-    lgupyy =   ( lgxx * lgzz - gxz(i,j,k) * gxz(i,j,k) ) / ldetg
+  Ayy = Ayy - F1o3 * gyy * trA
-    lgupyz = - ( lgxx * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / ldetg
+  Ayz = Ayz - F1o3 * gyz * trA
-    lgupzz =   ( lgxx * lgyy - gxy(i,j,k) * gxy(i,j,k) ) / ldetg
+  Azz = Azz - F1o3 * gzz * trA
-
+
-    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
+  detg = ONE / ( detg ** F1o3 ) 
-                 + lgupzz * Azz(i,j,k) &
+  
-         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
+  gxx = gxx * detg
-                 + lgupyz * Ayz(i,j,k))
+  gxy = gxy * detg
-
+  gxz = gxz * detg
-    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
+  gyy = gyy * detg
-    Axy(i,j,k) = Axy(i,j,k) - F1o3 * gxy(i,j,k) * ltrA
+  gyz = gyz * detg
-    Axz(i,j,k) = Axz(i,j,k) - F1o3 * gxz(i,j,k) * ltrA
+  gzz = gzz * detg
-    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
+
-    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * gyz(i,j,k) * ltrA
+  dxx = gxx - ONE
-    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
+  dyy = gyy - ONE
-
+  dzz = gzz - ONE
-    lscale = ONE / ( ldetg ** F1o3 )
-
-    dxx(i,j,k) = lgxx * lscale - ONE
-    gxy(i,j,k) = gxy(i,j,k) * lscale
-    gxz(i,j,k) = gxz(i,j,k) * lscale
-    dyy(i,j,k) = lgyy * lscale - ONE
-    gyz(i,j,k) = gyz(i,j,k) * lscale
-    dzz(i,j,k) = lgzz * lscale - ONE
-
-  enddo
-  enddo
-  enddo
 
   return
 
@@ -93,72 +81,52 @@
   real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Axx,Axy,Axz
   real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: Ayy,Ayz,Azz
 
-!~~~~~~~> Local variable:
+!~~~~~~~> Local variable:
-
+  
-  integer :: i,j,k
+  real*8, dimension(ex(1),ex(2),ex(3)) :: trA
-  real*8 :: lgxx,lgyy,lgzz,lscale
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
-  real*8 :: lgxy,lgxz,lgyz
+  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
-  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
+  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
-  real*8 :: ltrA
+
-  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0
+!~~~~~~>
-
+
-!~~~~~~>
+  gxx = dxx + ONE
-
+  gyy = dyy + ONE
-  do k=1,ex(3)
+  gzz = dzz + ONE
-  do j=1,ex(2)
+! for g
-  do i=1,ex(1)
+  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-
+           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
-! for g: normalize determinant first
+
-    lgxx = dxx(i,j,k) + ONE
+  gupzz = ONE / ( gupzz ** F1o3 ) 
-    lgyy = dyy(i,j,k) + ONE
+  
-    lgzz = dzz(i,j,k) + ONE
+  gxx = gxx * gupzz
-    lgxy = gxy(i,j,k)
+  gxy = gxy * gupzz
-    lgxz = gxz(i,j,k)
+  gxz = gxz * gupzz
-    lgyz = gyz(i,j,k)
+  gyy = gyy * gupzz
-
+  gyz = gyz * gupzz
-    lscale =  lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
+  gzz = gzz * gupzz
-            + lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
+
-            - lgxy * lgxy * lgzz - lgxx * lgyz * lgyz
+  dxx = gxx - ONE
-
+  dyy = gyy - ONE
-    lscale = ONE / ( lscale ** F1o3 )
+  dzz = gzz - ONE
-
+! for A  
-    lgxx = lgxx * lscale
+
-    lgxy = lgxy * lscale
+  gupxx =   ( gyy * gzz - gyz * gyz )
-    lgxz = lgxz * lscale
+  gupxy = - ( gxy * gzz - gyz * gxz )
-    lgyy = lgyy * lscale
+  gupxz =   ( gxy * gyz - gyy * gxz )
-    lgyz = lgyz * lscale
+  gupyy =   ( gxx * gzz - gxz * gxz )
-    lgzz = lgzz * lscale
+  gupyz = - ( gxx * gyz - gxy * gxz )
-
+  gupzz =   ( gxx * gyy - gxy * gxy )
-    dxx(i,j,k) = lgxx - ONE
+
-    gxy(i,j,k) = lgxy
+  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
-    gxz(i,j,k) = lgxz
+       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
-    dyy(i,j,k) = lgyy - ONE
+
-    gyz(i,j,k) = lgyz
+  Axx = Axx - F1o3 * gxx * trA
-    dzz(i,j,k) = lgzz - ONE
+  Axy = Axy - F1o3 * gxy * trA
-
+  Axz = Axz - F1o3 * gxz * trA
-! for A: trace-free using normalized metric (det=1, no division needed)
+  Ayy = Ayy - F1o3 * gyy * trA
-    lgupxx =   ( lgyy * lgzz - lgyz * lgyz )
+  Ayz = Ayz - F1o3 * gyz * trA
-    lgupxy = - ( lgxy * lgzz - lgyz * lgxz )
+  Azz = Azz - F1o3 * gzz * trA
-    lgupxz =   ( lgxy * lgyz - lgyy * lgxz )
-    lgupyy =   ( lgxx * lgzz - lgxz * lgxz )
-    lgupyz = - ( lgxx * lgyz - lgxy * lgxz )
-    lgupzz =   ( lgxx * lgyy - lgxy * lgxy )
-
-    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
-                 + lgupzz * Azz(i,j,k) &
-         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
-                 + lgupyz * Ayz(i,j,k))
-
-    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
-    Axy(i,j,k) = Axy(i,j,k) - F1o3 * lgxy * ltrA
-    Axz(i,j,k) = Axz(i,j,k) - F1o3 * lgxz * ltrA
-    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
-    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * lgyz * ltrA
-    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
-
-  enddo
-  enddo
-  enddo
 
   return
 

AMSS_NCKU_source/fmisc.f90

@@ -324,7 +324,8 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
    enddo
@@ -349,7 +350,8 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
       funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
@@ -377,7 +379,8 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
       funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
@@ -883,20 +886,17 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+  funcc = 0.d0
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+   do i=0,ord-1
-   do i=0,ord-1
+      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
-      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
+   enddo
-   enddo
+   do i=0,ord-1
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
-   do i=0,ord-1
+   enddo
-      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
+   do i=0,ord-1
-   enddo
+      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
-!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
+   enddo
-   do i=0,ord-1
-      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
-   enddo
 
 end subroutine symmetry_bd
 
@@ -912,7 +912,8 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
       funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -940,7 +941,8 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+  funcc = 0.d0
+  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
    do i=0,ord-1
       funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
       funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -1111,353 +1113,151 @@ end subroutine d2dump
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! common code for cell and vertex
+! common code for cell and vertex
-!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
-! Lagrangian polynomial interpolation
+! Lagrangian polynomial interpolation
-!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
-#ifndef POLINT6_USE_BARYCENTRIC
+
-#define POLINT6_USE_BARYCENTRIC 1
+  subroutine polint(xa,ya,x,y,dy,ordn)
-#endif
+
-
+  implicit none
-!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
+
-  subroutine polint6_neville(xa, ya, x, y, dy)
+!~~~~~~> Input Parameter:
-  implicit none
+  integer,intent(in) :: ordn
-
+  real*8, dimension(ordn), intent(in) :: xa,ya
-  real*8, dimension(6), intent(in) :: xa, ya
+  real*8, intent(in) :: x
-  real*8, intent(in) :: x
+  real*8, intent(out) :: y,dy
-  real*8, intent(out) :: y, dy
+
-
+!~~~~~~> Other parameter:
-  integer :: i, m, ns, n_m
+
-  real*8, dimension(6) :: c, d, ho
+  integer :: m,n,ns
-  real*8 :: dif, dift, hp, h, den_val
+  real*8, dimension(ordn) :: c,d,den,ho
-
+  real*8 :: dif,dift
-  c = ya
+
-  d = ya
+!~~~~~~>
-  ho = xa - x
+
-
+  n=ordn
-  ns = 1
+  m=ordn
-  dif = abs(x - xa(1))
+
-
+  c=ya
-  do i = 2, 6
+  d=ya
-    dift = abs(x - xa(i))
+  ho=xa-x
-    if (dift < dif) then
+
-      ns = i
+  ns=1
-      dif = dift
+  dif=abs(x-xa(1))
-    end if
+  do m=1,n
-  end do
+   dift=abs(x-xa(m))
-
+   if(dift < dif) then
-  y = ya(ns)
+    ns=m
-  ns = ns - 1
+    dif=dift
-
+   end if
-  do m = 1, 5
+  end do
-    n_m = 6 - m
+
-    do i = 1, n_m
+  y=ya(ns)
-      hp = ho(i)
+  ns=ns-1
-      h  = ho(i+m)
+  do m=1,n-1
-      den_val = hp - h
+    den(1:n-m)=ho(1:n-m)-ho(1+m:n)
-
+    if (any(den(1:n-m) == 0.0))then
-      if (den_val == 0.0d0) then
+      write(*,*) 'failure in polint for point',x
-        write(*,*) 'failure in polint for point',x
+      write(*,*) 'with input points: ',xa
-        write(*,*) 'with input points: ',xa
+      stop
-        stop
+    endif
-      end if
+    den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
-
+    d(1:n-m)=ho(1+m:n)*den(1:n-m)
-      den_val = (c(i+1) - d(i)) / den_val
+    c(1:n-m)=ho(1:n-m)*den(1:n-m)
-
+    if (2*ns < n-m) then
-      d(i) = h * den_val
+      dy=c(ns+1)
-      c(i) = hp * den_val
+    else
-    end do
+      dy=d(ns)
-
+      ns=ns-1
-    if (2 * ns < n_m) then
+    end if
-      dy = c(ns + 1)
+    y=y+dy
-    else
+  end do
-      dy = d(ns)
+
-      ns = ns - 1
+  return
-    end if
+
-    y = y + dy
+  end subroutine polint
-  end do
+!------------------------------------------------------------------------------
-
+!
-  return
+! interpolation in 2 dimensions, follow yx order
-  end subroutine polint6_neville
+!
-
+!------------------------------------------------------------------------------
-!DIR$ ATTRIBUTES FORCEINLINE :: polint6_barycentric
+  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
-  subroutine polint6_barycentric(xa, ya, x, y, dy)
+
-  implicit none
+  implicit none
-
+
-  real*8, dimension(6), intent(in) :: xa, ya
+!~~~~~~> Input parameters:
-  real*8, intent(in) :: x
+  integer,intent(in) :: ordn
-  real*8, intent(out) :: y, dy
+  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
-
+  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
-  integer :: i, j
+  real*8, intent(in) :: x1,x2
-  logical :: is_uniform
+  real*8, intent(out) :: y,dy
-  real*8, dimension(6) :: lambda
+
-  real*8 :: dx, den_i, term, num, den, step, tol
+!~~~~~~> Other parameters:
-  real*8, parameter :: c_uniform(6) = (/ -1.d0, 5.d0, -10.d0, 10.d0, -5.d0, 1.d0 /)
+
-
+  integer  :: i,m
-  do i = 1, 6
+  real*8, dimension(ordn) :: ymtmp
-    if (x == xa(i)) then
+  real*8, dimension(ordn) :: yntmp
-      y = ya(i)
+
-      dy = 0.d0
+  m=size(x1a)
-      return
+  
-    end if
+  do i=1,m
-  end do
+
-
+    yntmp=ya(i,:)
-  step = xa(2) - xa(1)
+    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  is_uniform = (step /= 0.d0)
+
-  if (is_uniform) then
+  end do
-    tol = 64.d0 * epsilon(1.d0) * max(1.d0, abs(step))
+
-    do i = 3, 6
+  call polint(x1a,ymtmp,x1,y,dy,ordn)
-      if (abs((xa(i) - xa(i-1)) - step) > tol) then
+
-        is_uniform = .false.
+  return
-        exit
+
-      end if
+  end subroutine polin2
-    end do
+!------------------------------------------------------------------------------
-  end if
+!
-
+! interpolation in 3 dimensions, follow zyx order
-  if (is_uniform) then
+!
-    num = 0.d0
+!------------------------------------------------------------------------------
-    den = 0.d0
+  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
-    do i = 1, 6
+
-      term = c_uniform(i) / (x - xa(i))
+  implicit none
-      num = num + term * ya(i)
+
-      den = den + term
+!~~~~~~> Input parameters:
-    end do
+  integer,intent(in) :: ordn
-    y = num / den
+  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
-    dy = 0.d0
+  real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
-    return
+  real*8, intent(in) :: x1,x2,x3
-  end if
+  real*8, intent(out) :: y,dy
-
+
-  do i = 1, 6
+!~~~~~~> Other parameters:
-    den_i = 1.d0
+
-    do j = 1, 6
+  integer  :: i,j,m,n
-      if (j /= i) then
+  real*8, dimension(ordn,ordn) :: yatmp
-        dx = xa(i) - xa(j)
+  real*8, dimension(ordn) :: ymtmp
-        if (dx == 0.0d0) then
+  real*8, dimension(ordn) :: yntmp
-          write(*,*) 'failure in polint for point',x
+  real*8, dimension(ordn) :: yqtmp
-          write(*,*) 'with input points: ',xa
+
-          stop
+  m=size(x1a)
-        end if
+  n=size(x2a)
-        den_i = den_i * dx
+  
-      end if
+  do i=1,m
-    end do
+   do j=1,n
-    lambda(i) = 1.d0 / den_i
+
-  end do
+    yqtmp=ya(i,j,:)
-
+    call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
-  num = 0.d0
+
-  den = 0.d0
+   end do
-  do i = 1, 6
+
-    term = lambda(i) / (x - xa(i))
+    yntmp=yatmp(i,:)
-    num = num + term * ya(i)
+    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-    den = den + term
+
-  end do
+  end do
-
+
-  y = num / den
+  call polint(x1a,ymtmp,x1,y,dy,ordn)
-  dy = 0.d0
+
-
+  return
-  return
+
-  end subroutine polint6_barycentric
+  end subroutine polin3
-
-!DIR$ ATTRIBUTES FORCEINLINE :: polint
-  subroutine polint(xa, ya, x, y, dy, ordn)
-  implicit none
-
-  integer, intent(in) :: ordn
-  real*8, dimension(ordn), intent(in) :: xa, ya
-  real*8, intent(in) :: x
-  real*8, intent(out) :: y, dy
-
-  integer :: i, m, ns, n_m
-  real*8, dimension(ordn) :: c, d, ho
-  real*8 :: dif, dift, hp, h, den_val
-
-  if (ordn == 6) then
-#if POLINT6_USE_BARYCENTRIC
-    call polint6_barycentric(xa, ya, x, y, dy)
-#else
-    call polint6_neville(xa, ya, x, y, dy)
-#endif
-    return
-  end if
-
-  c = ya
-  d = ya
-  ho = xa - x
-
-  ns = 1
-  dif = abs(x - xa(1))
-
-  do i = 2, ordn
-    dift = abs(x - xa(i))
-    if (dift < dif) then
-      ns = i
-      dif = dift
-    end if
-  end do
-
-  y = ya(ns)
-  ns = ns - 1
-
-  do m = 1, ordn - 1
-    n_m = ordn - m
-    do i = 1, n_m
-      hp = ho(i)
-      h  = ho(i+m)
-      den_val = hp - h
-
-      if (den_val == 0.0d0) then
-        write(*,*) 'failure in polint for point',x
-        write(*,*) 'with input points: ',xa
-        stop
-      end if
-
-      den_val = (c(i+1) - d(i)) / den_val
-
-      d(i) = h * den_val
-      c(i) = hp * den_val
-    end do
-
-    if (2 * ns < n_m) then
-      dy = c(ns + 1)
-    else
-      dy = d(ns)
-      ns = ns - 1
-    end if
-    y = y + dy
-  end do
-
-  return
-  end subroutine polint
-!------------------------------------------------------------------------------
-! Compute Lagrange interpolation basis weights for one target point.
-!------------------------------------------------------------------------------
-!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
-  subroutine polint_lagrange_weights(xa, x, w, ordn)
-  implicit none
-
-  integer, intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: xa
-  real*8, intent(in) :: x
-  real*8, dimension(1:ordn), intent(out) :: w
-
-  integer :: i, j
-  real*8 :: num, den, dx
-
-  do i = 1, ordn
-    num = 1.d0
-    den = 1.d0
-    do j = 1, ordn
-      if (j /= i) then
-        dx = xa(i) - xa(j)
-        if (dx == 0.0d0) then
-          write(*,*) 'failure in polint for point',x
-          write(*,*) 'with input points: ',xa
-          stop
-        end if
-        num = num * (x - xa(j))
-        den = den * dx
-      end if
-    end do
-    w(i) = num / den
-  end do
-
-  return
-  end subroutine polint_lagrange_weights
-!------------------------------------------------------------------------------
-!
-! interpolation in 2 dimensions, follow yx order
-!
-!------------------------------------------------------------------------------
-  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
-  implicit none
-
-  integer,intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
-  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
-  real*8, intent(in) :: x1,x2
-  real*8, intent(out) :: y,dy
-
-#ifdef POLINT_LEGACY_ORDER
-  integer  :: i,m
-  real*8, dimension(ordn) :: ymtmp
-  real*8, dimension(ordn) :: yntmp
-
-  m=size(x1a)
-  do i=1,m
-    yntmp=ya(i,:)
-    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  end do
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
-#else
-  integer  :: j
-  real*8, dimension(ordn) :: ymtmp
-  real*8 :: dy_temp
-
-  do j=1,ordn
-    call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
-  end do
-  call polint(x2a, ymtmp, x2, y, dy, ordn)
-#endif
-
-  return
-  end subroutine polin2
-!------------------------------------------------------------------------------
-!
-! interpolation in 3 dimensions, follow zyx order
-!
-!------------------------------------------------------------------------------
-  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
-  implicit none
-
-  integer,intent(in) :: ordn
-  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
-  real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
-  real*8, intent(in) :: x1,x2,x3
-  real*8, intent(out) :: y,dy
-
-#ifdef POLINT_LEGACY_ORDER
-  integer  :: i,j,m,n
-  real*8, dimension(ordn,ordn) :: yatmp
-  real*8, dimension(ordn) :: ymtmp
-  real*8, dimension(ordn) :: yntmp
-  real*8, dimension(ordn) :: yqtmp
-
-  m=size(x1a)
-  n=size(x2a)
-  do i=1,m
-   do j=1,n
-    yqtmp=ya(i,j,:)
-    call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
-   end do
-    yntmp=yatmp(i,:)
-    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-  end do
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
-#else
-  integer  :: i, j, k
-  real*8, dimension(ordn) :: w1, w2
-  real*8, dimension(ordn) :: ymtmp
-  real*8 :: yx_sum, x_sum
-
-  call polint_lagrange_weights(x1a, x1, w1, ordn)
-  call polint_lagrange_weights(x2a, x2, w2, ordn)
-
-  do k = 1, ordn
-    yx_sum = 0.d0
-    do j = 1, ordn
-      x_sum = 0.d0
-      do i = 1, ordn
-        x_sum = x_sum + w1(i) * ya(i,j,k)
-      end do
-      yx_sum = yx_sum + w2(j) * x_sum
-    end do
-    ymtmp(k) = yx_sum
-  end do
-
-  call polint(x3a, ymtmp, x3, y, dy, ordn)
-#endif
-
-  return
-  end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm  
   subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
@@ -1476,9 +1276,7 @@ end subroutine d2dump
   real*8            :: dX, dY, dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
+  integer::i,j,k
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
 
   dX = X(2) - X(1)
   dY = Y(2) - Y(1)
@@ -1502,91 +1300,15 @@ if(dabs(X(1)-xmin) < dX) imin = 1
 if(dabs(Y(1)-ymin) < dY) jmin = 1
 if(dabs(Z(1)-zmin) < dZ) kmin = 1
 
-  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
-  allocate(f_flat(n_elements))
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
-  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
-  deallocate(f_flat)
 
 f_out = f_out*dX*dY*dZ
 
   return
 
   end subroutine l2normhelper
-!--------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------
-  subroutine l2normhelper7(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
+! calculate L2norm especially for shell Blocks
-                           f1,f2,f3,f4,f5,f6,f7,f_out,gw)
-
-  implicit none
-!~~~~~~> Input parameters:
-  integer,intent(in ):: ex(1:3)
-  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
-  integer,intent(in)::gw
-  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f1,f2,f3,f4,f5,f6,f7
-  real*8, intent(out) :: f_out(7)
-!~~~~~~> Other variables:
-
-  real*8            :: dX, dY, dZ
-  integer::imin,jmin,kmin
-  integer::imax,jmax,kmax
-  integer::i,j,k
-  real*8 :: s1,s2,s3,s4,s5,s6,s7
-
-  dX = X(2) - X(1)
-  dY = Y(2) - Y(1)
-  dZ = Z(2) - Z(1)
-
-   imin = gw+1
-   jmin = gw+1
-   kmin = gw+1
-
-   imax = ex(1) - gw
-   jmax = ex(2) - gw
-   kmax = ex(3) - gw
-
-if(dabs(X(ex(1))-xmax) < dX) imax = ex(1)
-if(dabs(Y(ex(2))-ymax) < dY) jmax = ex(2)
-if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)
-if(dabs(X(1)-xmin) < dX) imin = 1
-if(dabs(Y(1)-ymin) < dY) jmin = 1
-if(dabs(Z(1)-zmin) < dZ) kmin = 1
-
-  s1 = 0.d0
-  s2 = 0.d0
-  s3 = 0.d0
-  s4 = 0.d0
-  s5 = 0.d0
-  s6 = 0.d0
-  s7 = 0.d0
-
-  do k=kmin,kmax
-    do j=jmin,jmax
-!DIR$ SIMD REDUCTION(+:s1,s2,s3,s4,s5,s6,s7)
-      do i=imin,imax
-        s1 = s1 + f1(i,j,k)*f1(i,j,k)
-        s2 = s2 + f2(i,j,k)*f2(i,j,k)
-        s3 = s3 + f3(i,j,k)*f3(i,j,k)
-        s4 = s4 + f4(i,j,k)*f4(i,j,k)
-        s5 = s5 + f5(i,j,k)*f5(i,j,k)
-        s6 = s6 + f6(i,j,k)*f6(i,j,k)
-        s7 = s7 + f7(i,j,k)*f7(i,j,k)
-      enddo
-    enddo
-  enddo
-
-  f_out(1) = s1*dX*dY*dZ
-  f_out(2) = s2*dX*dY*dZ
-  f_out(3) = s3*dX*dY*dZ
-  f_out(4) = s4*dX*dY*dZ
-  f_out(5) = s5*dX*dY*dZ
-  f_out(6) = s6*dX*dY*dZ
-  f_out(7) = s7*dX*dY*dZ
-
-  return
-
-  end subroutine l2normhelper7
-!--------------------------------------------------------------------------------------
-! calculate L2norm especially for shell Blocks
   subroutine l2normhelper_sh(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                           f,f_out,gw,ogw,Symmetry)
 
@@ -1603,9 +1325,7 @@ if(dabs(Z(1)-zmin) < dZ) kmin = 1
   real*8            :: dX, dY, dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
+  integer::i,j,k
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
 
   real*8 :: PIo4
 
@@ -1668,11 +1388,7 @@ if(Symmetry==2)then
   if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif
 
-  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
-  allocate(f_flat(n_elements))
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
-  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
-  deallocate(f_flat)
 
 f_out = f_out*dX*dY*dZ
 
@@ -1699,9 +1415,7 @@ f_out = f_out*dX*dY*dZ
   real*8            :: dX, dY, dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
-  integer::i,j,k
+  integer::i,j,k
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
 
   real*8 :: PIo4
 
@@ -1764,11 +1478,11 @@ if(Symmetry==2)then
   if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif
 
-Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
-  allocate(f_flat(Nout))
+
-  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [Nout])
+f_out = f_out
-  f_out = DDOT(Nout, f_flat, 1, f_flat, 1)
+
-  deallocate(f_flat)
+Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 
   return
 
@@ -1869,12 +1583,9 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 !       ^
 ! f=3/8*f_1 + 3/4*f_2 - 1/8*f_3
 
-  real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
+  real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
-  integer :: i,j,k
+
-
+  fout = C1*f1+C2*f2+C3*f3
-  do concurrent (k=1:ext(3), j=1:ext(2), i=1:ext(1))
-    fout(i,j,k) = C1*f1(i,j,k)+C2*f2(i,j,k)+C3*f3(i,j,k)
-  end do
 
   return
 
@@ -1968,8 +1679,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   real*8, dimension(ORDN,ORDN,ORDN) :: ya
   real*8, dimension(ORDN,ORDN) :: tmp2
   real*8, dimension(ORDN) :: tmp1
-  real*8, dimension(3) :: SoAh
+  real*8, dimension(3) :: SoAh
-  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   cxB = inds+1  
@@ -2015,7 +1725,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
     tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
   enddo
 
-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo
 
   return
 
@@ -2044,8 +1757,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   integer,dimension(2) :: cxB,cxT
   real*8, dimension(ORDN,ORDN) :: ya
   real*8, dimension(ORDN) :: tmp1
-  real*8, dimension(2) :: SoAh
+  real*8, dimension(2) :: SoAh
-  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   cxB = inds(1:2)+1  
@@ -2080,7 +1792,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
     tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
   enddo
 
-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo
 
   return
 
@@ -2106,12 +1821,11 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 !~~~~~~> Other parameters:
 
   real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,ex(3)) :: fh
-  integer :: m
+  integer :: m
-  integer :: cxB,cxT
+  integer :: cxB,cxT
-  real*8, dimension(ORDN) :: ya
+  real*8, dimension(ORDN) :: ya
-  real*8 :: SoAh
+  real*8 :: SoAh
-  integer,dimension(3) :: inds
+  integer,dimension(3) :: inds
-  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   inds = indsi + 1
@@ -2172,7 +1886,10 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
           write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
   endif
 
-  f_int = DDOT(ORDN, coef, 1, ya, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*ya(m)
+  enddo
 
   return
 
@@ -2408,28 +2125,20 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   implicit none
   integer,intent(in) :: N
 
-  real*8 :: gont
+  real*8 :: gont
-
+
-  integer :: i
+  integer :: i
-  real*8, parameter, dimension(0:20) :: fact_table = [ &
-    1.d0, 1.d0, 2.d0, 6.d0, 24.d0, 120.d0, 720.d0, 5040.d0, 40320.d0, &
-    362880.d0, 3628800.d0, 39916800.d0, 479001600.d0, 6227020800.d0, &
-    87178291200.d0, 1307674368000.d0, 20922789888000.d0, &
-    355687428096000.d0, 6402373705728000.d0, 121645100408832000.d0, &
-    2432902008176640000.d0 ]
 
 ! sanity check
-  if(N < 0)then
+  if(N < 0)then
-     write(*,*) "ffact: error input for factorial"
+     write(*,*) "ffact: error input for factorial"
-     gont = 1.d0
+     return
-     return
+  endif
-  endif
+
-
+  gont = 1.d0
-  if(N <= 20)then
+  do i=1,N
-     gont = fact_table(N)
+     gont = gont*i
-  else
+  enddo
-     gont = exp(log_gamma(dble(N+1)))
-  endif
 
   return
 

AMSS_NCKU_source/fmisc.h

@@ -12,10 +12,9 @@
 #define f_global_interpind global_interpind
 #define f_global_interpind2d global_interpind2d
 #define f_global_interpind1d global_interpind1d
-#define f_l2normhelper l2normhelper
+#define f_l2normhelper l2normhelper
-#define f_l2normhelper7 l2normhelper7
+#define f_l2normhelper_sh l2normhelper_sh
-#define f_l2normhelper_sh l2normhelper_sh
+#define f_l2normhelper_sh_rms l2normhelper_sh_rms
-#define f_l2normhelper_sh_rms l2normhelper_sh_rms
 #define f_average average
 #define f_average3 average3
 #define f_average2 average2
@@ -42,10 +41,9 @@
 #define f_global_interpind GLOBAL_INTERPIND
 #define f_global_interpind2d GLOBAL_INTERPIND2D
 #define f_global_interpind1d GLOBAL_INTERPIND1D
-#define f_l2normhelper L2NORMHELPER
+#define f_l2normhelper L2NORMHELPER
-#define f_l2normhelper7 L2NORMHELPER7
+#define f_l2normhelper_sh L2NORMHELPER_SH
-#define f_l2normhelper_sh L2NORMHELPER_SH
+#define f_l2normhelper_sh_rms L2NORMHELPER_SH_RMS
-#define f_l2normhelper_sh_rms L2NORMHELPER_SH_RMS
 #define f_average AVERAGE
 #define f_average3 AVERAGE3
 #define f_average2 AVERAGE2
@@ -72,10 +70,9 @@
 #define f_global_interpind global_interpind_
 #define f_global_interpind2d global_interpind2d_
 #define f_global_interpind1d global_interpind1d_
-#define f_l2normhelper l2normhelper_
+#define f_l2normhelper l2normhelper_
-#define f_l2normhelper7 l2normhelper7_
+#define f_l2normhelper_sh l2normhelper_sh_
-#define f_l2normhelper_sh l2normhelper_sh_
+#define f_l2normhelper_sh_rms l2normhelper_sh_rms_
-#define f_l2normhelper_sh_rms l2normhelper_sh_rms_
 #define f_average average_
 #define f_average3 average3_
 #define f_average2 average2_
@@ -159,30 +156,21 @@ extern "C"
 							  int *, double *, int &, int &);
 }
 
-extern "C"
+extern "C"
-{
+{
-	void f_l2normhelper(int *, double *, double *, double *,
+	void f_l2normhelper(int *, double *, double *, double *,
-						double &, double &, double &,
+						double &, double &, double &,
-						double &, double &, double &,
+						double &, double &, double &,
-						double *, double &, int &);
+						double *, double &, int &);
-}
+}
-
+
-extern "C"
+extern "C"
-{
+{
-	void f_l2normhelper7(int *, double *, double *, double *,
+	void f_l2normhelper_sh(int *, double *, double *, double *,
-						 double &, double &, double &,
+						   double &, double &, double &,
-						 double &, double &, double &,
+						   double &, double &, double &,
-						 double *, double *, double *, double *,
+						   double *, double &, int &, int &, int &);
-						 double *, double *, double *, double *, int &);
+}
-}
-
-extern "C"
-{
-	void f_l2normhelper_sh(int *, double *, double *, double *,
-						   double &, double &, double &,
-						   double &, double &, double &,
-						   double *, double &, int &, int &, int &);
-}
 
 extern "C"
 {

AMSS_NCKU_source/macrodef.h

@@ -2,7 +2,7 @@
 #ifndef MICRODEF_H
 #define MICRODEF_H
 
-#include "macrodef.fh"
+#include "microdef.fh"
 
 // application parameters
 

AMSS_NCKU_source/makefile

@@ -1,25 +1,11 @@
-
+
-
+
-include makefile.inc
+include makefile.inc
-
+
-## polint(ordn=6) kernel selector:
+.SUFFIXES: .o .f90 .C .for .cu
-##   1 (default): barycentric fast path
+
-##   0          : fallback to Neville path
+.f90.o:
-POLINT6_USE_BARY ?= 1
+	$(f90) $(f90appflags) -c $< -o $@
-POLINT6_FLAG = -DPOLINT6_USE_BARYCENTRIC=$(POLINT6_USE_BARY)
-
-ARCH_OPT = -march=x86-64-v4
-CXXAPPFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
-TP_OPTFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
-              -Dfortran3 -Dnewc -I${MKLROOT}/include
-
-.SUFFIXES: .o .f90 .C .for .cu
-
-.f90.o:
-	$(f90) $(f90appflags) -c $< -o $@
 
 .C.o:
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
@@ -27,14 +13,8 @@ TP_OPTFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
 .for.o:
 	$(f77) -c $< -o $@
 
-.cu.o:
+.cu.o:
-	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
+	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
-
-TwoPunctures.o: TwoPunctures.C
-	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
-
-TwoPunctureABE.o: TwoPunctureABE.C
-	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
 
 # Input files
 C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
@@ -115,8 +95,8 @@ ABE: $(C++FILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
 ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)
 
-TwoPunctureABE: $(TwoPunctureFILES)
+TwoPunctureABE: $(TwoPunctureFILES)
-	$(CLINKER) $(TP_OPTFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(TwoPunctureFILES) $(LDLIBS)
 
 clean:
 	rm *.o ABE ABEGPU TwoPunctureABE make.log -f

AMSS_NCKU_source/makefile.inc

@@ -1,34 +1,22 @@
 
-## GCC version (commented out)
+filein  = -I/usr/include -I/usr/include/openmpi-x86_64 -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
-## filein  = -I/usr/include -I/usr/lib/x86_64-linux-gnu/mpich/include -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
-## filein  = -I/usr/include/ -I/usr/include/openmpi-x86_64/ -I/usr/lib/x86_64-linux-gnu/openmpi/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
-## LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran
 
-## Intel oneAPI version with oneMKL
+##filein  = -I/usr/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/ -I/usr/lib/cuda/include
-filein  = -I/usr/include/ -I${MKLROOT}/include
 
-## Use sequential oneMKL to avoid introducing extra OpenMP behavior into ABE.
+LDLIBS  = -L/usr/lib64/openmpi/lib -Wl,-rpath,/usr/lib64/openmpi/lib -lmpi -lgfortran -L/usr/local/cuda-13.1/lib64 -Wl,-rpath,/usr/local/cuda-13.1/lib64 -lcudart -lcuda 
-LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl -liomp5
+##LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -L/usr/lib/cuda/lib64 -lcudart -lmpi -lgfortran
 
-## Optional Intel oneTBB allocator, kept aligned with main's build environment.
+CXXAPPFLAGS  = -O3 -Wno-deprecated -Dfortran3 -Dnewc
-USE_TBBMALLOC ?= 1
+#f90appflags = -O3 -fpp
-TBBMALLOC_SO ?= /home/intel/oneapi/2025.3/lib/libtbbmalloc.so
+f90appflags  = -O3 -x f95-cpp-input
-ifneq ($(wildcard $(TBBMALLOC_SO)),)
+f90          = gfortran 
-TBBMALLOC_LIBS = -Wl,--no-as-needed $(TBBMALLOC_SO) -Wl,--as-needed
+f77          = gfortran 
-else
+CXX          = g++
-TBBMALLOC_LIBS = -Wl,--no-as-needed -ltbbmalloc -Wl,--as-needed
+CC           = gcc
-endif
+CLINKER      = mpic++ 
-ifeq ($(USE_TBBMALLOC),1)
-LDLIBS := $(TBBMALLOC_LIBS) $(LDLIBS)
-endif
 
-f90          = ifx
+Cu = /usr/local/cuda-13.1/bin/nvcc
-f77          = ifx
+CUDA_LIB_PATH = -L/usr/local/cuda-13.1/lib64 -I/usr/include -I/usr/local/cuda-13.1/include
-CXX          = icpx
-CC           = icx
-CLINKER      = mpiicpx
-
-Cu = nvcc
-CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include
 #CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -arch compute_13 -code compute_13,sm_13 -Dfortran3 -Dnewc
-CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -Dfortran3 -Dnewc
+# RTX 4050 uses Ada Lovelace architecture (compute capability 8.9)
+CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -arch=sm_89 -Dfortran3 -Dnewc

makefile_and_run.py

@@ -28,7 +28,7 @@ def makefile_ABE():
 
     ## Build command
     if (input_data.GPU_Calculation == "no"):
-        makefile_command  = "make -j96" + " ABE"
+        makefile_command  = "make -j4" + " ABE"
     elif (input_data.GPU_Calculation == "yes"):
         makefile_command  = "make -j4" + " ABEGPU"
     else:

parallel_plot_helper.py

@@ -1,12 +0,0 @@
-import multiprocessing
-
-
-def run_plot_task(task):
-    func, args = task
-    return func(*args)
-
-
-def run_plot_tasks_parallel(plot_tasks):
-    ctx = multiprocessing.get_context('fork')
-    with ctx.Pool() as pool:
-        pool.map(run_plot_task, plot_tasks)

plot_GW_strain_amplitude_xiaoqu.py

@@ -8,13 +8,11 @@
 ##
 #################################################
 
-import numpy                               ## numpy for array operations
+import numpy                               ## numpy for array operations
-import scipy                               ## scipy for interpolation and signal processing
+import scipy                               ## scipy for interpolation and signal processing
-import math
+import math
-import matplotlib
+import matplotlib.pyplot    as     plt     ## matplotlib for plotting
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
+import os                                  ## os for system/file operations
-import matplotlib.pyplot    as     plt     ## matplotlib for plotting
-import os                                  ## os for system/file operations
 
 import AMSS_NCKU_Input as input_data
 

plot_binary_data.py

@@ -6,22 +6,17 @@
 ## Author: Xiaoqu
 ## Dates: 2024/10/01 --- 2025/09/14
 ##
-#################################################
+#################################################
-
+
-## Restrict OpenMP to one thread per process so that parallel
+import numpy
-## subprocess plotting does not multiply BLAS thread counts.
+import scipy
-import os
+import matplotlib.pyplot    as     plt
-os.environ.setdefault("OMP_NUM_THREADS", "1")
+from   matplotlib.colors    import LogNorm
-
+from   mpl_toolkits.mplot3d import Axes3D
-import numpy
+## import torch
-import scipy
+import AMSS_NCKU_Input      as input_data
-import matplotlib
+
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
+import os
-import matplotlib.pyplot    as     plt
-from   matplotlib.colors    import LogNorm
-from   mpl_toolkits.mplot3d import Axes3D
-## import torch
-import AMSS_NCKU_Input      as input_data
 
 
 #########################################################################################
@@ -97,9 +92,9 @@ def plot_binary_data( filename, binary_outdir, figure_outdir ):
 
 
 
-####################################################################################
+####################################################################################
-
+
-# Plot a single binary dataset (2D slices and 3D surface)
+# Plot a single binary dataset (2D slices and 3D surface)
 
 def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):
 
@@ -193,15 +188,7 @@ def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):
     plt.savefig( os.path.join(figure_surfaceplot_outdir, figure_title + " time = " + str(time) + " surface_plot.pdf") )   # save figure
     plt.close()
 
-    return
+    return
-
+
-####################################################################################
+####################################################################################
-
-## Allow standalone subprocess execution for parallel binary-data plotting.
-if __name__ == '__main__':
-    import sys
-    if len(sys.argv) != 4:
-        print(f"Usage: {sys.argv[0]} <filename> <binary_outdir> <figure_outdir>")
-        sys.exit(1)
-    plot_binary_data(sys.argv[1], sys.argv[2], sys.argv[3])
 

plot_xiaoqu.py

@@ -6,20 +6,15 @@
 ## 2024/10/01 --- 2025/09/14
 ##
 #################################################
-
+
-import numpy                               ## numpy for array operations
+import numpy                               ## numpy for array operations
-import matplotlib
+import matplotlib.pyplot    as     plt     ## matplotlib for plotting
-matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
+from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
-import matplotlib.pyplot    as     plt     ## matplotlib for plotting
+import glob
-from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
+import os                                  ## operating system utilities
-import glob
+
-import os                                  ## operating system utilities
+import plot_binary_data
-
+import AMSS_NCKU_Input as input_data
-import plot_binary_data
-import AMSS_NCKU_Input as input_data
-import subprocess
-import sys
-import multiprocessing
 
 # plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots
 
@@ -55,37 +50,13 @@ def generate_binary_data_plot( binary_outdir, figure_outdir ):
         file_list.append(x)
         print(x)
 
-    ## Plot each file in parallel using subprocesses.
+    ## Plot each file in the list
-    ## Each subprocess starts with BLAS thread limits in plot_binary_data.py.
+    for filename in file_list:
-    script = os.path.join( os.path.dirname(__file__), "plot_binary_data.py" )
+        print(filename)
-    max_workers = min( multiprocessing.cpu_count(), len(file_list) ) if file_list else 0
+        plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)
-
+
-    running = []
+    print(                        )
-    failed  = []
+    print( " Binary Data Plot Has been Finished " )
-    for filename in file_list:
-        print(filename)
-        proc = subprocess.Popen(
-            [sys.executable, script, filename, binary_outdir, figure_outdir],
-        )
-        running.append( (proc, filename) )
-        if len(running) >= max_workers:
-            p, fn = running.pop(0)
-            p.wait()
-            if p.returncode != 0:
-                failed.append(fn)
-
-    for p, fn in running:
-        p.wait()
-        if p.returncode != 0:
-            failed.append(fn)
-
-    if failed:
-        print( " WARNING: the following binary data plots failed:" )
-        for fn in failed:
-            print( "   ", fn )
-
-    print(                        )
-    print( " Binary Data Plot Has been Finished " )
     print(                                        )
 
     return