202504 update

AMSS_NCKU_Program.py

-
+ 
 ##################################################################
 ##
 ## AMSS-NCKU 数值相对论启动程序
@@ -145,10 +145,11 @@ print(" AMSS-NCKU 程序的宏文件 macrodef.fh 已生成 ")
 
 ##  根据用户的要求编译 AMSS-NCKU 程序
 
-print(                                                         )
-print( " 准备根据要求编译并运行 AMSS-NCKU 程序，按回车继续！！！ " )
-inputvalue = input()           
-print()
+print(                                         )
+print( " 准备根据要求编译并运行 AMSS-NCKU 程序 " )
+print(                                         )
+#inputvalue = input()           
+#print()
 
 AMSS_NCKU_source_path = "AMSS_NCKU_source"
 AMSS_NCKU_source_copy = os.path.join(File_directionary, "AMSS_NCKU_source_copy")
@@ -365,6 +366,15 @@ shutil.copy( AMSS_NCKU_ADM_data,        os.path.join(output_directionary, "bssn_
 shutil.copy( AMSS_NCKU_psi4_data,       os.path.join(output_directionary, "bssn_psi4.dat"      ) )
 shutil.copy( AMSS_NCKU_constraint_data, os.path.join(output_directionary, "bssn_constraint.dat") )
 
+## 程序其它输出
+if (input_data.Equation_Class == "BSSN-EM"):
+    AMSS_NCKU_phi1_data = os.path.join(binary_results_directionary, "bssn_phi1.dat" )
+    AMSS_NCKU_phi2_data = os.path.join(binary_results_directionary, "bssn_phi2.dat" )
+    shutil.copy( AMSS_NCKU_phi1_data, os.path.join(output_directionary, "bssn_phi1.dat" ) )
+    shutil.copy( AMSS_NCKU_phi2_data, os.path.join(output_directionary, "bssn_phi2.dat" ) )
+elif (input_data.Equation_Class == "BSSN-EScalar"):
+    AMSS_NCKU_maxs_data = os.path.join(binary_results_directionary, "bssn_maxs.dat" )
+    shutil.copy( AMSS_NCKU_maxs_data, os.path.join(output_directionary, "bssn_maxs.dat" ) )
 
 ##################################################################
 
@@ -381,6 +391,9 @@ import plot_xiaoqu
 plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directionary, figure_directionary )
 plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directionary, figure_directionary )
 
+## 画出黑洞间距随时间的变化图
+plot_xiaoqu.generate_puncture_distence_plot( binary_results_directionary, figure_directionary )
+
 ## 画出引力波波形图
 for i in range(input_data.Detector_Number):
     plot_xiaoqu.generate_gravitational_waveform_plot( binary_results_directionary, figure_directionary, i )

AMSS_NCKU_source/Ansorg.h

-//$Id: Ansorg.h,v 1.2 2012/04/03 10:49:40 zjcao Exp $
-#ifndef Ansorg_H
-#define Ansorg_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>
-
-#define PI M_PI
-
-class Ansorg
-{
-protected:
-  int n1, n2, n3, ntotal;
-  int order;
-  double *coordA, *coordB, *coordphi;
-  int ps_rxx, ps_rxy, ps_ryx, ps_ryy;
-  double ps_b, ps_dx;
-  double PIh;
-  double *pu_ps;
-  int myrank;
-
-public:
-  Ansorg(char *filename, int orderi);
-  ~Ansorg();
-  double ps_u_at_xyz(double x, double y, double z);
-  void set_ABp();
-  void xyz_to_ABp(double x, double y, double z,
-                  double *A, double *B, double *phi);
-  double interpolate_tri_bar(double x, double y, double z,
-                             int n1, int n2, int n3,
-                             double *x1, double *x2, double *x3, double *yp);
-  int find_point_bisection(double x, int n, double *xp, int o);
-  void barycentric_omega(int n, int s, double *x, double *omega);
-  double barycentric(double x0, int n, int s, double *x, double *y,
-                     double *omega);
-};
-#endif /* Ansorg_H */
+ 
+#ifndef Ansorg_H
+#define Ansorg_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>
+
+#define PI M_PI
+
+class Ansorg
+{
+protected:
+  int n1, n2, n3, ntotal;
+  int order;
+  double *coordA, *coordB, *coordphi;
+  int ps_rxx, ps_rxy, ps_ryx, ps_ryy;
+  double ps_b, ps_dx;
+  double PIh;
+  double *pu_ps;
+  int myrank;
+
+public:
+  Ansorg(char *filename, int orderi);
+  ~Ansorg();
+  double ps_u_at_xyz(double x, double y, double z);
+  void set_ABp();
+  void xyz_to_ABp(double x, double y, double z,
+                  double *A, double *B, double *phi);
+  double interpolate_tri_bar(double x, double y, double z,
+                             int n1, int n2, int n3,
+                             double *x1, double *x2, double *x3, double *yp);
+  int find_point_bisection(double x, int n, double *xp, int o);
+  void barycentric_omega(int n, int s, double *x, double *omega);
+  double barycentric(double x0, int n, int s, double *x, double *y,
+                     double *omega);
+};
+#endif /* Ansorg_H */

AMSS_NCKU_source/BH_diagnostics.h

-#ifndef BH_DIAGNOSTICS_H
-#define BH_DIAGNOSTICS_H
-namespace AHFinderDirect
-{
-
-	struct BH_diagnostics
-	{
-	public:
-		// mean x,y,z
-		fp centroid_x, centroid_y, centroid_z;
-
-		// these are quadrupole moments about the centroid, i.e.
-		// mean(xi*xj) - centroid_i*centroid_j
-		fp quadrupole_xx, quadrupole_xy, quadrupole_xz,
-			quadrupole_yy, quadrupole_yz,
-			quadrupole_zz;
-
-		// min,max,mean surface radius about local coordinate origin
-		fp min_radius, max_radius, mean_radius;
-
-		// xyz bounding box
-		fp min_x, max_x,
-			min_y, max_y,
-			min_z, max_z;
-
-		// proper circumference
-		// (computed using induced metric along these local-coordinate planes)
-		fp circumference_xy,
-			circumference_xz,
-			circumference_yz;
-
-		// surface area (computed using induced metric)
-		// and quantities derived from it
-		fp area, irreducible_mass, areal_radius;
-
-		double Px, Py, Pz, Sx, Sy, Sz;
-
-	public:
-		// position of diagnostics in buffer and number of diagnostics
-		enum
-		{
-			posn__centroid_x = 0,
-			posn__centroid_y,
-			posn__centroid_z,
-			posn__quadrupole_xx,
-			posn__quadrupole_xy,
-			posn__quadrupole_xz,
-			posn__quadrupole_yy,
-			posn__quadrupole_yz,
-			posn__quadrupole_zz,
-			posn__min_radius,
-			posn__max_radius,
-			posn__mean_radius,
-
-			posn__min_x,
-			posn__max_x,
-			posn__min_y,
-			posn__max_y,
-			posn__min_z,
-			posn__max_z,
-
-			posn__circumference_xy,
-			posn__circumference_xz,
-			posn__circumference_yz,
-
-			posn__area,
-			posn__irreducible_mass,
-			posn__areal_radius,
-
-			N_buffer // no comma	// size of buffer
-		};
-
-		// copy diagnostics to/from buffer
-		void copy_to_buffer(double buffer[N_buffer]) const;
-		void copy_from_buffer(const double buffer[N_buffer]);
-
-	public:
-		void compute(patch_system &ps);
-
-		void compute_signature(patch_system &ps, const double dT);
-
-		FILE *setup_output_file(int N_horizons, int hn)
-			const;
-
-		void output(FILE *fileptr, double time)
-			const;
-
-		BH_diagnostics();
-
-	private:
-		static double surface_integral(const patch_system &ps,
-									   int src_gfn, bool src_gfn_is_even_across_xy_plane,
-									   bool src_gfn_is_even_across_xz_plane,
-									   bool src_gfn_is_even_across_yz_plane,
-									   enum patch::integration_method method);
-	};
-
-	//******************************************************************************
-
-} // namespace AHFinderDirect
-#endif /* BH_DIAGNOSTICS_H */
+
+#ifndef BH_DIAGNOSTICS_H
+#define BH_DIAGNOSTICS_H
+namespace AHFinderDirect
+{
+
+	struct BH_diagnostics
+	{
+	public:
+		// mean x,y,z
+		fp centroid_x, centroid_y, centroid_z;
+
+		// these are quadrupole moments about the centroid, i.e.
+		// mean(xi*xj) - centroid_i*centroid_j
+		fp quadrupole_xx, quadrupole_xy, quadrupole_xz,
+			quadrupole_yy, quadrupole_yz,
+			quadrupole_zz;
+
+		// min,max,mean surface radius about local coordinate origin
+		fp min_radius, max_radius, mean_radius;
+
+		// xyz bounding box
+		fp min_x, max_x,
+			min_y, max_y,
+			min_z, max_z;
+
+		// proper circumference
+		// (computed using induced metric along these local-coordinate planes)
+		fp circumference_xy,
+			circumference_xz,
+			circumference_yz;
+
+		// surface area (computed using induced metric)
+		// and quantities derived from it
+		fp area, irreducible_mass, areal_radius;
+
+		double Px, Py, Pz, Sx, Sy, Sz;
+
+	public:
+		// position of diagnostics in buffer and number of diagnostics
+		enum
+		{
+			posn__centroid_x = 0,
+			posn__centroid_y,
+			posn__centroid_z,
+			posn__quadrupole_xx,
+			posn__quadrupole_xy,
+			posn__quadrupole_xz,
+			posn__quadrupole_yy,
+			posn__quadrupole_yz,
+			posn__quadrupole_zz,
+			posn__min_radius,
+			posn__max_radius,
+			posn__mean_radius,
+
+			posn__min_x,
+			posn__max_x,
+			posn__min_y,
+			posn__max_y,
+			posn__min_z,
+			posn__max_z,
+
+			posn__circumference_xy,
+			posn__circumference_xz,
+			posn__circumference_yz,
+
+			posn__area,
+			posn__irreducible_mass,
+			posn__areal_radius,
+
+			N_buffer // no comma	// size of buffer
+		};
+
+		// copy diagnostics to/from buffer
+		void copy_to_buffer(double buffer[N_buffer]) const;
+		void copy_from_buffer(const double buffer[N_buffer]);
+
+	public:
+		void compute(patch_system &ps);
+
+		void compute_signature(patch_system &ps, const double dT);
+
+		FILE *setup_output_file(int N_horizons, int hn)
+			const;
+
+		void output(FILE *fileptr, double time)
+			const;
+
+		BH_diagnostics();
+
+	private:
+		static double surface_integral(const patch_system &ps,
+									   int src_gfn, bool src_gfn_is_even_across_xy_plane,
+									   bool src_gfn_is_even_across_xz_plane,
+									   bool src_gfn_is_even_across_yz_plane,
+									   enum patch::integration_method method);
+	};
+
+	//******************************************************************************
+
+} // namespace AHFinderDirect
+#endif /* BH_DIAGNOSTICS_H */

AMSS_NCKU_source/Block.h

-//$Id: Block.h,v 1.5 2013/05/18 09:17:31 zjcao Exp $
-#ifndef BLOCK_H
-#define BLOCK_H
-
-#include <mpi.h>
-#include "macrodef.h" //need dim here; Vertex or Cell
-#include "var.h"
-#include "MyList.h"
-class Block
-{
-
-public:
-   int shape[dim];
-   double bbox[2 * dim];
-   double *X[dim];
-   int rank; // where the real data locate in
-   int lev, cgpu;
-   int ingfs, fngfs;
-   int *(*igfs);
-   double *(*fgfs);
-
-public:
-   Block() {};
-   Block(int DIM, int *shapei, double *bboxi, int ranki, int ingfsi, int fngfs, int levi, const int cgpui = 0);
-
-   ~Block();
-
-   void checkBlock();
-
-   double getdX(int dir);
-   void swapList(MyList<var> *VarList1, MyList<var> *VarList2, int myrank);
-};
-
-#endif /* BLOCK_H */
+ 
+#ifndef BLOCK_H
+#define BLOCK_H
+
+#include <mpi.h>
+#include "macrodef.h" //need dim here; Vertex or Cell
+#include "var.h"
+#include "MyList.h"
+class Block
+{
+
+public:
+   int shape[dim];
+   double bbox[2 * dim];
+   double *X[dim];
+   int rank; // where the real data locate in
+   int lev, cgpu;
+   int ingfs, fngfs;
+   int *(*igfs);
+   double *(*fgfs);
+
+public:
+   Block() {};
+   Block(int DIM, int *shapei, double *bboxi, int ranki, int ingfsi, int fngfs, int levi, const int cgpui = 0);
+
+   ~Block();
+
+   void checkBlock();
+
+   double getdX(int dir);
+   void swapList(MyList<var> *VarList1, MyList<var> *VarList2, int myrank);
+};
+
+#endif /* BLOCK_H */

AMSS_NCKU_source/IntPnts0.C

-
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdarg.h>
-#include <string.h>
-
-#include <mpi.h>
-
-#include "myglobal.h"
-
-int CCTK_VInfo(const char *thorn, const char *format, ...)
-{
-   int myrank;
-   MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
-   if (myrank !=0) return 0;
-   
-   va_list ap;
-   va_start (ap, format);
-   fprintf (stdout, "INFO (%s): ", thorn);
-   vfprintf (stdout, format, ap);
-   fprintf (stdout, "\n");
-   va_end (ap);
-   return 0;
-}
-int CCTK_VWarn (int level,
-                int line,
-                const char *file,
-                const char *thorn,
-                const char *format,
-                ...) 
-{  
-   int myrank;
-   MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
-   if (myrank !=0) return 0;
-   
-   va_list ap;
-   va_start (ap, format);
-   fprintf (stdout, "WARN (%s): ", thorn);
-   vfprintf (stdout, format, ap);
-   fprintf (stdout, "\n");
-   va_end (ap);
-   return 0;
-}
+ 
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+
+#include <mpi.h>
+
+#include "myglobal.h"
+
+int CCTK_VInfo(const char *thorn, const char *format, ...)
+{
+   int myrank;
+   MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+   if (myrank !=0) return 0;
+   
+   va_list ap;
+   va_start (ap, format);
+   fprintf (stdout, "INFO (%s): ", thorn);
+   vfprintf (stdout, format, ap);
+   fprintf (stdout, "\n");
+   va_end (ap);
+   return 0;
+}
+int CCTK_VWarn (int level,
+                int line,
+                const char *file,
+                const char *thorn,
+                const char *format,
+                ...) 
+{  
+   int myrank;
+   MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+   if (myrank !=0) return 0;
+   
+   va_list ap;
+   va_start (ap, format);
+   fprintf (stdout, "WARN (%s): ", thorn);
+   vfprintf (stdout, format, ap);
+   fprintf (stdout, "\n");
+   va_end (ap);
+   return 0;
+}

AMSS_NCKU_source/bssn2adm.f90

-!$Id: bssn2adm.f90,v 1.1 2012/12/01 04:03:15 zjcao Exp $
-!-------------------------------------------------------------------------------!
-! convert bssn variables to ADM variables                                       !
-!-------------------------------------------------------------------------------!
-  subroutine bssn2adm(ex,chi,trK, &
-               gxx,gxy,gxz,gyy,gyz,gzz, &
-               Axx,Axy,Axz,Ayy,Ayz,Azz, &
-               adm_gxx,adm_gxy,adm_gxz,adm_gyy,adm_gyz,adm_gzz, &
-               Kxx,Kxy,Kxz,Kyy,Kyz,Kzz)
-
-  implicit none
-!~~~~~~> Input parameters:
-
-  integer,intent(in ):: ex(1:3)
-  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::chi,trK
-  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::gxx,gxy,gxz,gyy,gyz,gzz
-  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::Axx,Axy,Axz,Ayy,Ayz,Azz
-
-  real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: adm_gxx,adm_gxy,adm_gxz,adm_gyy,adm_gyz,adm_gzz
-  real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Kxx,Kxy,Kxz,Kyy,Kyz,Kzz
-
-  real*8, parameter :: F1o3=1.d0/3.d0
-
-  adm_gxx = gxx/chi
-  adm_gxy = gxy/chi
-  adm_gxz = gxz/chi
-  adm_gyy = gyy/chi
-  adm_gyz = gyz/chi
-  adm_gzz = gzz/chi
-
-  Kxx = Axx/chi+F1o3*trK*adm_gxx
-  Kxy = Axy/chi+F1o3*trK*adm_gxy
-  Kxz = Axz/chi+F1o3*trK*adm_gxz
-  Kyy = Ayy/chi+F1o3*trK*adm_gyy
-  Kyz = Ayz/chi+F1o3*trK*adm_gyz
-  Kzz = Azz/chi+F1o3*trK*adm_gzz
-
-  return
-
-  end subroutine bssn2adm
+ 
+!-------------------------------------------------------------------------------!
+! convert bssn variables to ADM variables                                       !
+!-------------------------------------------------------------------------------!
+  subroutine bssn2adm(ex,chi,trK, &
+               gxx,gxy,gxz,gyy,gyz,gzz, &
+               Axx,Axy,Axz,Ayy,Ayz,Azz, &
+               adm_gxx,adm_gxy,adm_gxz,adm_gyy,adm_gyz,adm_gzz, &
+               Kxx,Kxy,Kxz,Kyy,Kyz,Kzz)
+
+  implicit none
+!~~~~~~> Input parameters:
+
+  integer,intent(in ):: ex(1:3)
+  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::chi,trK
+  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::gxx,gxy,gxz,gyy,gyz,gzz
+  double precision,intent(in),dimension(ex(1),ex(2),ex(3))::Axx,Axy,Axz,Ayy,Ayz,Azz
+
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: adm_gxx,adm_gxy,adm_gxz,adm_gyy,adm_gyz,adm_gzz
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Kxx,Kxy,Kxz,Kyy,Kyz,Kzz
+
+  real*8, parameter :: F1o3=1.d0/3.d0
+
+  adm_gxx = gxx/chi
+  adm_gxy = gxy/chi
+  adm_gxz = gxz/chi
+  adm_gyy = gyy/chi
+  adm_gyz = gyz/chi
+  adm_gzz = gzz/chi
+
+  Kxx = Axx/chi+F1o3*trK*adm_gxx
+  Kxy = Axy/chi+F1o3*trK*adm_gxy
+  Kxz = Axz/chi+F1o3*trK*adm_gxz
+  Kyy = Ayy/chi+F1o3*trK*adm_gyy
+  Kyz = Ayz/chi+F1o3*trK*adm_gyz
+  Kzz = Azz/chi+F1o3*trK*adm_gzz
+
+  return
+
+  end subroutine bssn2adm