更新

AMSS_NCKU_Input.py

+
+#################################################
+##
+## 这个文件包含了数值相对论所需要的输入
+## 小曲
+## 2024/03/19 --- 2025/02/27
+##
+#################################################
+
+import numpy   ## 导入 numpy 包 
+
+#################################################
+
+## 设置程序运行目录和计算资源
+
+File_directionary   = "xiaoqu_Results_testBBH"      ## 程序运行目录
+Output_directionary = "output_file"                 ## 存放二进制数据的子目录
+MPI_processes       = 64                            ## 想要调用的进程数目
+
+GPU_Calculation     = "no"                          ## 是否开启 GPU 计算，可选 yes 或 no
+CPU_Part            = 0.5
+GPU_Part            = 0.5
+
+#################################################
+
+
+#################################################
+
+## 设置程序计算方法
+
+Symmetry                 = "equatorial-symmetry"   ## 系统对称性，可选 equatorial-symmetry、no-symmetry
+Equation_Class           = "BSSN"                  ## 设置方程形式，可选 BSSN、Z4C、BSSN-EScalar、BSSN-EM
+                                                   ##       BSSN 和 Z4C   适合于 GR 旋转黑洞的真空计算
+                                                   ##       BSSN-EM      涉及 GR 带电黑洞的真空计算
+                                                   ##       BSSN-EScalar 涉及到标量张量-F(R) 理论的计算，需要在后面设定额外参数
+                                                   ## 注意：GPU 计算仅支持 BSSN
+                                                   ## 这里选择 BSSN-EScalar，需要在后面设定 F(R) 理论的参数
+Initial_Data_Method      = "Ansorg-TwoPuncture"    ## 设置求解数值相对论初值的方法
+                                                   ## 可选 Ansorg-TwoPuncture、Lousto-Analytical、KerrSchild-Analytical
+                                                   ## 注意：当前 BSSN-EM 的计算不支持用解析公式 Lousto-Analytical
+Time_Evolution_Method    = "runge-kutta-45"        ## 时间演化方法，可选 runge-kutta-45
+Finite_Diffenence_Method = "6th-order"             ## 有限差分方法，可选 2nd-order、4th-order、6th-order、8th-order
+
+#################################################
+
+
+#################################################
+
+## 设置时间演化信息
+
+Start_Evolution_Time     = 0.0                    ## 起始演化时间
+Final_Evolution_Time     = 1500.0                  ## 最终演化时间
+Check_Time               = 1000.0
+Dump_Time                = 50.0                   ## 每隔一定时间间隔储存数据
+D2_Dump_Time             = 300.0
+Analysis_Time            = 0.1
+Evolution_Step_Number    = 10000000               ## 时间迭代次数
+Courant_Factor           = 0.5                    ## Courant 因子（决定每一步时间演化的时间间隔）
+Dissipation              = 0.2                    ## 耗散因子
+
+#################################################
+
+
+#################################################
+
+## 设置多层格点信息
+
+basic_grid_set    = "Patch"                          ## 设定网格类型，可选 Patch 和 Shell-Patch
+grid_center_set   = "Cell"                           ## 网格中心设置，可选 Cell 和 Vertex
+
+grid_level        = 11                               ## 设置格点的总层数
+static_grid_level = 6                                ## 设置静态格点的层数
+moving_grid_level = grid_level - static_grid_level   ## 可移动格点的层数
+
+analysis_level    = 0
+refinement_level  = 4                                ## 从该层开始进行时间细化
+
+grid_maxvalue     =  500.0                           ## 设置最外层格点的坐标最大值
+grid_minvalue     = - grid_maxvalue                  ## 设置最外层格点的坐标最小值ֵ
+
+static_grid_number = 96                              ## 设置固定格点每一层每一维数的格点数目
+moving_grid_number = 48                              ## 设置可移动格点每一层每一维数的格点数目
+shell_grid_number  = [32, 32, 100]                    ## 设置最外层球状网格（shell patch）的格点数目
+                                                     ## 以 phi、theta、r 的顺序给定
+devide_factor      = 2.0                             ## 设置相邻两层网格分辨率的比例（不要轻易改变）
+static_grid_type   = 'Linear'                        ## 设置固定格点的类型，可选 'Linear'
+moving_grid_type   = 'Linear'                        ## 设置固定格点的类型，可选 'Linear'
+
+quarter_sphere_number = 64                           ## 1/4 球面积分的格点数目
+
+#################################################
+
+
+#################################################
+
+## 设置黑洞 puncture （穿刺法）的信息
+
+puncture_number       = 2                                     ## 设置 puncture 的数目
+
+position_BH           = numpy.zeros( (puncture_number, 3) )   ## 初始化每个黑洞的初始位置
+parameter_BH          = numpy.zeros( (puncture_number, 3) )   ## 初始化每个黑洞的参数
+dimensionless_spin_BH = numpy.zeros( (puncture_number, 3) )   ## 初始化每个黑洞的无量纲自旋
+momentum_BH           = numpy.zeros( (puncture_number, 3) )   ## 初始化每个黑洞的动量
+
+puncture_data_set     = "TwoPuncture"                            ## 设置双星轨道坐标的方式，可选 Manually 和 TwoPuncture
+
+#---------------------------------------------
+
+## 如果设置双星初始轨道坐标的方式选为 TwoPuncture，只需要给定黑洞参数，偏心率，距离即可
+
+## 这一步要与初值求解中的 Ansorg-TwoPuncture 配合使用，否则可能会报错
+## （也就是说，用 Ansorg-TwoPuncture 求解初值，轨道坐标设置可以设置 Manually 和 TwoPuncture，
+##  但双星轨道坐标如果设置为 Manually 而不是 TwoPuncture，最好使用其它方法求解数值初值，否则可能会使 TwoPuncture 程序无法正确读入输入而报错）
+
+Distance = 11.0
+e0       = 0.0
+
+## 设置每个黑洞的参数 (M Q* a*)  
+## 质量  无量纲电荷  无量纲自旋
+parameter_BH[0] = [ 0.5,  0.0,  0.0 ]   
+parameter_BH[1] = [ 0.5,  0.0, -0.0 ]  
+## 注意，如果求解数值相对论初值的方法选为 Ansorg-TwoPuncture ，第一个黑洞必须为质量较大的那个，且则黑洞总质量会自动 rescale 为 M=1 （其它情况下必须手动 rescale）
+
+## 设置每个黑洞的无量纲自旋
+## 无对称性时 ，需要手动给 3 个方向的自旋角动量
+dimensionless_spin_BH[0] = [ 0.0,  0.0,  0.0 ]   
+dimensionless_spin_BH[1] = [ 0.0,  0.0, -0.0 ]  
+
+## 注意，如果设置双星初始轨道坐标的方式选为 TwoPuncture，则程序自动调整将较大质量黑洞放在 y 轴正向，将较小质量黑洞放在 y 轴负向
+##       如果设置双星初始轨道坐标的方式选为 Manually，且则需要手动调整到 y 轴方向 
+## use Brugmann's convention
+##  -----0-----> y
+##   -      +     
+
+#---------------------------------------------
+
+## 如果设置 puncture data 的方式选为 Manually，需要手动给定所有黑洞参数
+
+## 设置每个黑洞的初始位置
+position_BH[0]  = [  0.0,  +5.5,  0.0 ]  
+position_BH[1]  = [  0.0,  -5.5,  0.0 ]  
+
+## 设置每个黑洞的动量信息  
+momentum_BH[0]  = [ -0.09011, -0.000704, 0.0 ]
+momentum_BH[1]  = [ +0.09011, +0.000704, 0.0 ] 
+
+
+#################################################
+
+
+#################################################
+
+## 设置引力波和探测器的相关信息
+
+GW_L_max        = 4                      ## 引力波最大的 L
+GW_M_max        = 4                      ## 引力波最大的 M
+Detector_Number = 11                     ## 探测器的数目
+Detector_Rmin   = 50.0                   ## 最近探测器的距离
+Detector_Rmax   = 150.0                  ## 最远探测器的距离
+
+#################################################
+
+
+#################################################
+
+## 设置表观视界的参数
+
+AHF_Find       = "no"                    ## 是否开启表观视界计算，可选 yes 或 no
+
+AHF_Find_Every = 24
+AHF_Dump_Time  = 20.0
+
+#################################################
+                                   

AMSS_NCKU_source/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 */

AMSS_NCKU_source/Block.C

+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cstdlib>
+#include <cstdio>
+#include <string>
+#include <cmath>
+#include <new>
+using namespace std;
+
+#include "Block.h"
+#include "misc.h"
+
+Block::Block(int DIM, int *shapei, double *bboxi, int ranki, int ingfsi, int fngfsi, int levi, const int cgpui) : rank(ranki), ingfs(ingfsi), fngfs(fngfsi), lev(levi), cgpu(cgpui)
+{
+  for (int i = 0; i < dim; i++)
+    X[i] = 0;
+
+  if (DIM != dim)
+  {
+    cout << "dimension is not consistent in Block construction" << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+  bool flag = false;
+  for (int i = 0; i < dim; i++)
+  {
+    shape[i] = shapei[i];
+    if (shape[i] <= 0)
+      flag = true;
+    bbox[i] = bboxi[i];
+    bbox[dim + i] = bboxi[dim + i];
+  }
+
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  if (flag)
+  {
+    cout << "myrank: " << myrank << ", on rank: " << rank << endl;
+    cout << "error shape in Block construction: (" << shape[0] << "," << shape[1] << "," << shape[2] << ")" << endl;
+    cout << "box boundary: (" << bbox[0] << ":" << bbox[3] << "," << bbox[1] << ":" << bbox[4] << "," << bbox[2] << ":" << bbox[5] << ")" << endl;
+    cout << "belong to level " << lev << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+#ifndef FAKECHECK
+  if (myrank == rank)
+  {
+    for (int i = 0; i < dim; i++)
+    {
+      X[i] = new double[shape[i]];
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+      double h = (bbox[dim + i] - bbox[i]) / (shape[i] - 1);
+      for (int j = 0; j < shape[i]; j++)
+        X[i][j] = bbox[i] + j * h;
+#else
+#ifdef Cell
+      double h = (bbox[dim + i] - bbox[i]) / shape[i];
+      for (int j = 0; j < shape[i]; j++)
+        X[i][j] = bbox[i] + (j + 0.5) * h;
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+    }
+
+    int nn = shape[0] * shape[1] * shape[2];
+    fgfs = new double *[fngfs];
+    for (int i = 0; i < fngfs; i++)
+    {
+      fgfs[i] = (double *)malloc(sizeof(double) * nn);
+      if (!(fgfs[i]))
+      {
+        cout << "on node#" << rank << ", out of memory when constructing Block." << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+      memset(fgfs[i], 0, sizeof(double) * nn);
+    }
+
+    igfs = new int *[ingfs];
+    for (int i = 0; i < ingfs; i++)
+    {
+      igfs[i] = (int *)malloc(sizeof(int) * nn);
+      if (!(igfs[i]))
+      {
+        cout << "on node#" << rank << ", out of memory when constructing Block." << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
+      }
+      memset(igfs[i], 0, sizeof(int) * nn);
+    }
+  }
+#endif
+}
+Block::~Block()
+{
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  if (myrank == rank)
+  {
+    for (int i = 0; i < dim; i++)
+      delete[] X[i];
+    for (int i = 0; i < ingfs; i++)
+      free(igfs[i]);
+    delete[] igfs;
+    for (int i = 0; i < fngfs; i++)
+      free(fgfs[i]);
+    delete[] fgfs;
+    X[0] = X[1] = X[2] = 0;
+    igfs = 0;
+    fgfs = 0;
+  }
+}
+void Block::checkBlock()
+{
+  int myrank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+  if (myrank == 0)
+  {
+    cout << "belong to level " << lev << endl;
+    cout << "shape: [";
+    for (int i = 0; i < dim; i++)
+    {
+      cout << shape[i];
+      if (i < dim - 1)
+        cout << ",";
+      else
+        cout << "]";
+    }
+    cout << " resolution: [";
+    for (int i = 0; i < dim; i++)
+    {
+      cout << getdX(i);
+      if (i < dim - 1)
+        cout << ",";
+      else
+        cout << "]" << endl;
+    }
+    cout << "locate on node " << rank << ", at (includes ghost zone):" << endl;
+    cout << "(";
+    for (int i = 0; i < dim; i++)
+    {
+      cout << bbox[i] << ":" << bbox[dim + i];
+      if (i < dim - 1)
+        cout << ",";
+      else
+        cout << ")" << endl;
+    }
+    cout << "has " << ingfs << " int type grids functions," << fngfs << " double type grids functions" << endl;
+  }
+}
+double Block::getdX(int dir)
+{
+  if (dir < 0 || dir >= dim)
+  {
+    cout << "Block::getdX: error input dir = " << dir << ", this Block has direction (0," << dim - 1 << ")" << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  double h;
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+  if (shape[dir] == 1)
+  {
+    cout << "Block::getdX: for direction " << dir << ", this Block has only one point. Can not determine dX for vertex center grid." << endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+  h = (bbox[dim + dir] - bbox[dir]) / (shape[dir] - 1);
+#else
+#ifdef Cell
+  h = (bbox[dim + dir] - bbox[dir]) / shape[dir];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+  return h;
+}
+void Block::swapList(MyList<var> *VarList1, MyList<var> *VarList2, int myrank)
+{
+  if (rank == myrank)
+  {
+    MyList<var> *varl1 = VarList1, *varl2 = VarList2;
+    while (varl1 && varl2)
+    {
+      misc::swap<double *>(fgfs[varl1->data->sgfn], fgfs[varl2->data->sgfn]);
+      varl1 = varl1->next;
+      varl2 = varl2->next;
+    }
+    if (varl1 || varl2)
+    {
+      cout << "error in Block::swaplist, var lists does not match." << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+  }
+}

AMSS_NCKU_source/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/DataCT.C

+
+//-----------------------------------------------------------------------
+// Read binary files and do fancy things with them...
+//-----------------------------------------------------------------------
+#ifdef newc
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+#include <cstdlib>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+using namespace std;
+#else
+#include <math.h>
+#include <iostream.h>
+#include <iomanip.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <fstream.h>
+#endif
+
+#include "microdef.fh"
+
+int main(int argc, char *argv[])
+{
+  //
+  // USE:  DataCT flag file1 [ file2 ]
+  //
+  // where: - flag can be XY,XZ,YZ
+  //
+  void set_fname(char *fname);
+
+  if (argc < 3)
+  {
+    cout << "\aUsage: DataCT flag binaryfile1 [ binaryfile2 ] \n "
+         << "   where: - flag can be XY,XZ,YZ"
+         << endl;
+    exit(1);
+  }
+  ifstream infile1;
+  infile1.open(argv[2]);
+  if (!infile1)
+  {
+    cerr << "\a Can't open " << argv[2] << " for input." << endl;
+    exit(1);
+  }
+
+  /* read properties of the binary file */
+  double time;
+  int nx, ny, nz;
+  double xmin, xmax, ymin, ymax, zmin, zmax;
+  infile1.seekg(0, ios::beg);
+  infile1.read((char *)&time, sizeof(double));
+  infile1.read((char *)&nx, sizeof(int));
+  infile1.read((char *)&ny, sizeof(int));
+  infile1.read((char *)&nz, sizeof(int));
+  infile1.read((char *)&xmin, sizeof(double));
+  infile1.read((char *)&xmax, sizeof(double));
+  infile1.read((char *)&ymin, sizeof(double));
+  infile1.read((char *)&ymax, sizeof(double));
+  infile1.read((char *)&zmin, sizeof(double));
+  infile1.read((char *)&zmax, sizeof(double));
+
+  /* get rid of any 4 character suffix */
+  set_fname(argv[2]);
+
+  /* sanity check */
+  if (nx != ny || nx != nz)
+  {
+    cout << "\n"
+         << endl;
+    cout << " nx, ny and nz do not agree! Using a symmetry?... ";
+    cout << "\n"
+         << endl;
+  }
+
+  cout << "\n Reading file : " << argv[2] << endl;
+  cout << "\n  Time        : " << time << endl;
+  cout << "  Dimensions  : " << setw(16) << nx << setw(16) << ny << setw(16) << nz << endl;
+  cout << "   xmin, xmax : " << setw(16) << xmin << setw(16) << xmax << endl;
+  cout << "   ymin, ymax : " << setw(16) << ymin << setw(16) << ymax << endl;
+  cout << "   zmin, zmax : " << setw(16) << zmin << setw(16) << zmax << endl;
+  cout << "\n";
+
+  double *data;
+  data = new double[nx * ny * nz];
+  int i = 0, j = 0, k = 0;
+  infile1.read((char *)data, nx * ny * nz * sizeof(double));
+  infile1.close();
+  //
+  //
+  // if second file given, open second file and subtract from first one!
+  //
+  //
+  if (argc == 4)
+  {
+    infile1.open(argv[3]);
+    if (!infile1)
+    {
+      cerr << "\a Can't open " << argv[3] << " for input." << endl;
+      exit(1);
+    }
+    double *indata;
+    indata = new double[nx * ny * nz];
+    // read in header
+    infile1.seekg(0, ios::beg);
+    int nxin, nyin, nzin;
+    infile1.read((char *)&time, sizeof(double));
+    infile1.read((char *)&nxin, sizeof(int));
+    infile1.read((char *)&nyin, sizeof(int));
+    infile1.read((char *)&nzin, sizeof(int));
+    infile1.read((char *)&xmin, sizeof(double));
+    infile1.read((char *)&xmax, sizeof(double));
+    infile1.read((char *)&ymin, sizeof(double));
+    infile1.read((char *)&ymax, sizeof(double));
+    infile1.read((char *)&zmin, sizeof(double));
+    infile1.read((char *)&zmax, sizeof(double));
+    if (nxin != nx || nyin != ny || nzin != nz)
+    {
+      cerr << "\a Number of indices do not agree! " << endl;
+      exit(1);
+    }
+    cout << " Comparing with data at time " << time << "\n"
+         << endl;
+    infile1.read((char *)indata, nx * ny * nz * sizeof(double));
+    infile1.close();
+    for (i = 0; i < nx * ny * nz; i++)
+      data[i] -= indata[i];
+  }
+
+  double *X, *Y, *Z;
+  X = new double[nx];
+  Y = new double[ny];
+  Z = new double[nz];
+  double dd;
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+  dd = (xmax - xmin) / (nx - 1);
+  for (i = 0; i < nx; i++)
+    X[i] = xmin + i * dd;
+  dd = (ymax - ymin) / (ny - 1);
+  for (j = 0; j < ny; j++)
+    Y[j] = ymin + j * dd;
+  dd = (zmax - zmin) / (nz - 1);
+  for (k = 0; k < nz; k++)
+    Z[k] = zmin + k * dd;
+#else
+#ifdef Cell
+  dd = (xmax - xmin) / nx;
+  for (i = 0; i < nx; i++)
+    X[i] = xmin + (i + 0.5) * dd;
+  dd = (ymax - ymin) / ny;
+  for (j = 0; j < ny; j++)
+    Y[j] = ymin + (j + 0.5) * dd;
+  dd = (zmax - zmin) / nz;
+  for (k = 0; k < nz; k++)
+    Z[k] = zmin + (k + 0.5) * dd;
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+
+  int ext[3];
+  ext[0] = nx;
+  ext[1] = ny;
+  ext[2] = nz;
+  void writefile(int *ext, double *XX, double *YY, double *ZZ, double *datain,
+                 char *filename, const char *flag);
+  writefile(ext, X, Y, Z, data, argv[2], argv[1]);
+
+  delete[] data;
+  delete[] X;
+  delete[] Y;
+  delete[] Z;
+}
+
+/*-----------------------------------*/
+/* get rid of any 4 character suffix */
+/*-----------------------------------*/
+void set_fname(char *fname)
+{
+  int len = strlen(fname) - 4;
+  char *n_fname;
+  n_fname = new char[len];
+
+  for (int i = 0; i < len; ++i)
+  {
+    n_fname[i] = fname[i];
+    //     cout << n_fname[i] << " " << i << endl;
+  }
+  n_fname[len] = '\0';
+
+  //      cout << "n_fname: " << n_fname << " fname: " << fname << ", "
+  // 	  << len << endl;
+
+  strcpy(fname, n_fname); /* Send back the old pointer */
+  delete n_fname;
+}
+//|----------------------------------------------------------------------------
+//  writefile
+//|----------------------------------------------------------------------------
+void writefile(int *ext, double *XX, double *YY, double *ZZ, double *datain,
+               char *filename, const char *flag)
+{
+  int nx = ext[0], ny = ext[1], nz = ext[2];
+  int i, j, k;
+  char filename_h[50];
+  //|--->open out put file
+  ofstream outfile;
+
+  if (!strcmp(flag, "YZ"))
+  {
+    for (i = 0; i < nx; i++)
+    {
+      sprintf(filename_h, "%s_%d.dat", filename, i);
+      outfile.open(filename_h);
+      outfile << "# CT along X at " << i << endl;
+      for (k = 0; k < nz; k++)
+      {
+        for (j = 0; j < ny; j++)
+        {
+          outfile << setw(10) << setprecision(10) << YY[j] << " "
+                  << setw(10) << setprecision(10) << ZZ[k] << " "
+                  << datain[i + j * nx + k * nx * ny] << " "
+                  << endl;
+        }
+        outfile << "\n"; /* blanck line for gnuplot */
+      }
+      outfile.close();
+    }
+  }
+  else if (!strcmp(flag, "XZ"))
+  {
+    for (j = 0; j < ny; j++)
+    {
+      sprintf(filename_h, "%s_%d.dat", filename, j);
+      outfile.open(filename_h);
+      outfile << "# CT along Y at " << j << endl;
+      for (k = 0; k < nz; k++)
+      {
+        for (i = 0; i < nx; i++)
+        {
+          outfile << setw(10) << setprecision(10) << XX[i] << " "
+                  << setw(10) << setprecision(10) << ZZ[k] << " "
+                  << datain[i + j * nx + k * nx * ny] << " "
+                  << endl;
+        }
+        outfile << "\n"; /* blanck line for gnuplot */
+      }
+      outfile.close();
+    }
+  }
+  else if (!strcmp(flag, "XY"))
+  {
+    for (k = 0; k < nz; k++)
+    {
+      sprintf(filename_h, "%s_%d.dat", filename, k);
+      outfile.open(filename_h);
+      outfile << "# CT along Z at " << k << endl;
+      for (j = 0; j < ny; j++)
+      {
+        for (i = 0; i < nx; i++)
+        {
+          outfile << setw(10) << setprecision(10) << XX[i] << " "
+                  << setw(10) << setprecision(10) << YY[j] << " "
+                  << datain[i + j * nx + k * nx * ny] << " "
+                  << endl;
+        }
+        outfile << "\n"; /* blanck line for gnuplot */
+      }
+      outfile.close();
+    }
+  }
+  else
+  {
+    cout << "In output_data: not recognized flag-->" << flag << endl;
+    exit(0);
+  }
+}

AMSS_NCKU_source/FFT.f90

+
+
+#if 0
+program checkFFT
+use dfport
+implicit none
+double precision::x
+integer,parameter::N=256
+double precision,dimension(N*2)::p
+double precision,dimension(N/2)::s
+integer::ncount,j,idum
+character(len=8)::tt
+tt=clock()
+idum=iachar(tt(8:8))-48
+p=0.0
+open(77,file='prime.dat',status='unknown')
+loop1:do ncount=1,N
+   x=ran(idum)
+   p(2*ncount-1)=x
+   write(77,'(f15.3)')x
+enddo loop1
+close(77)
+call four1(p,N,1)
+do j=1,N/2
+  s(j)=p(2*j)*p(2*j)+p(2*j-1)*p(2*j-1)
+enddo
+x=0.0
+do j=1,N/2
+  x=x+s(j)
+enddo
+s=s/x
+open(77,file='power.dat',status='unknown')
+do j=1,N/2
+  write(77,'(2(1x,f15.3))')dble(j-1)/dble(N),s(j)
+enddo
+close(77)
+end program checkFFT
+#endif
+
+!-------------
+SUBROUTINE four1(dataa,nn,isign)
+implicit none
+INTEGER::isign,nn
+double precision,dimension(2*nn)::dataa
+INTEGER::i,istep,j,m,mmax,n
+double precision::tempi,tempr
+DOUBLE PRECISION::theta,wi,wpi,wpr,wr,wtemp
+n=2*nn
+j=1
+do i=1,n,2
+  if(j.gt.i)then
+     tempr=dataa(j)
+     tempi=dataa(j+1)
+     dataa(j)=dataa(i)
+     dataa(j+1)=dataa(i+1)
+     dataa(i)=tempr
+     dataa(i+1)=tempi
+  endif
+  m=nn
+1 if ((m.ge.2).and.(j.gt.m)) then
+  j=j-m
+  m=m/2
+goto 1
+  endif
+j=j+m
+enddo
+mmax=2
+2  if (n.gt.mmax) then
+     istep=2*mmax
+     theta=6.28318530717959d0/(isign*mmax)
+     wpr=-2.d0*sin(0.5d0*theta)**2
+     wpi=sin(theta)
+     wr=1.d0
+     wi=0.d0
+     do m=1,mmax,2
+       do i=m,n,istep
+         j=i+mmax
+         tempr=sngl(wr)*dataa(j)-sngl(wi)*dataa(j+1)
+         tempi=sngl(wr)*dataa(j+1)+sngl(wi)*dataa(j)
+         dataa(j)=dataa(i)-tempr
+         dataa(j+1)=dataa(i+1)-tempi
+         dataa(i)=dataa(i)+tempr
+         dataa(i+1)=dataa(i+1)+tempi
+       enddo
+          wtemp=wr
+          wr=wr*wpr-wi*wpi+wr
+          wi=wi*wpr+wtemp*wpi+wi
+     enddo
+mmax=istep
+goto 2
+endif
+return
+END SUBROUTINE four1

AMSS_NCKU_source/IntPnts.C

+//$Id: IntPnts.C,v 1.1 2012/04/03 10:49:42 zjcao Exp $
+
+#include "macrodef.h"
+#ifdef With_AHF
+
+#include <math.h>
+#include <stdio.h>
+
+#include <iostream>
+using namespace std;
+
+#include "myglobal.h"
+
+namespace AHFinderDirect
+{
+  extern struct state state;
+  int globalInterpGFL(double *X, double *Y, double *Z, int Ns,
+                      double *Data)
+  {
+    if (Ns == 0)
+      return 0;
+    int n;
+    double *pox[3];
+    for (int i = 0; i < 3; i++)
+      pox[i] = new double[Ns];
+    for (n = 0; n < Ns; n++)
+    {
+      pox[0][n] = X[n];
+      pox[1][n] = Y[n];
+      pox[2][n] = Z[n];
+    }
+
+    const int InList = 35;
+
+    double *datap;
+    datap = new double[Ns * InList];
+    if (!(state.ADM->AH_Interp_Points(state.AHList, Ns, pox, datap, state.Symmetry)))
+      return 0;
+    // reform data
+    for (int pnt = 0; pnt < Ns; pnt++)
+      for (int ii = 0; ii < InList; ii++)
+      {
+        if (ii == 0 || ii == 12 || ii == 20)
+          Data[pnt + ii * Ns] = datap[ii + pnt * InList] + 1;
+        else if (ii == 24) // from chi-1 to psi
+          Data[pnt + ii * Ns] = pow(datap[ii + pnt * InList] + 1, -0.25);
+        else if (ii == 25 || ii == 26 || ii == 27) // from chi,i to psi,i
+          Data[pnt + ii * Ns] = -pow(datap[24 + pnt * InList] + 1, -1.25) / 4 * datap[ii + pnt * InList];
+        else
+          Data[pnt + ii * Ns] = datap[ii + pnt * InList];
+      }
+    delete[] datap;
+
+    delete[] pox[0];
+    delete[] pox[1];
+    delete[] pox[2];
+
+    return 1;
+  }
+  // inerpolate lapse and shift
+  int globalInterpGFLlash(double *X, double *Y, double *Z, int Ns,
+                          double *Data)
+  {
+    if (Ns == 0)
+      return 0;
+    int n;
+    double *pox[3];
+    for (int i = 0; i < 3; i++)
+      pox[i] = new double[Ns];
+    for (n = 0; n < Ns; n++)
+    {
+      pox[0][n] = X[n];
+      pox[1][n] = Y[n];
+      pox[2][n] = Z[n];
+    }
+
+    double SYM = 1.0, ANT = -1.0;
+    const int InList = 4;
+
+    double *datap;
+    datap = new double[Ns * InList];
+    state.ADM->AH_Interp_Points(state.GaugeList, Ns, pox, datap, state.Symmetry);
+    // reform data
+    for (int pnt = 0; pnt < Ns; pnt++)
+      for (int ii = 0; ii < InList; ii++)
+        Data[pnt + ii * Ns] = datap[ii + pnt * InList];
+
+    delete[] datap;
+    delete[] pox[0];
+    delete[] pox[1];
+    delete[] pox[2];
+
+    return 1;
+  }
+
+} // namespace AHFinderDirect
+#endif

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;
+}

AMSS_NCKU_source/Jacobian.C

+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <string.h>
+
+#include "util_Table.h"
+#include "cctk.h"
+
+#include "config.h"
+#include "stdc.h"
+
+#include "util.h"
+#include "array.h"
+#include "cpm_map.h"
+#include "linear_map.h"
+
+#include "coords.h"
+#include "tgrid.h"
+#include "fd_grid.h"
+#include "patch.h"
+#include "patch_edge.h"
+#include "patch_interp.h"
+#include "ghost_zone.h"
+#include "patch_system.h"
+
+#include "Jacobian.h"
+#include "ilucg.h"
+// all the code in this file is inside this namespace
+namespace AHFinderDirect
+{
+	// this represents a single element stored in the matrix for
+	// sort_row_into_column_order()  and  sort_row_into_column_order__cmp()
+	struct matrix_element
+	{
+		int JA;
+		fp A;
+	};
+
+	Jacobian::Jacobian(patch_system &ps)
+		: ps_(ps),
+		  N_rows_(ps.N_grid_points()),
+		  N_nonzeros_(0), current_N_rows_(0), N_nonzeros_allocated_(0),
+		  IA_(new integer[N_rows_ + 1]), JA_(NULL), A_(NULL),
+		  itemp_(NULL), rtemp_(NULL)
+	{
+		IO_ = 1;
+		zero_matrix();
+	}
+
+	Jacobian::~Jacobian()
+	{
+		if (A_)
+			delete[] A_;
+		if (JA_)
+			delete[] JA_;
+		if (IA_)
+			delete[] IA_;
+		if (rtemp_)
+			delete[] rtemp_;
+		if (itemp_)
+			delete[] itemp_;
+	}
+
+	double Jacobian::element(int II, int JJ)
+		const
+	{
+		const int posn = find_element(II, JJ);
+		return (posn >= 0) ? A_[posn] : 0.0;
+	}
+
+	void Jacobian::zero_matrix()
+	{
+
+		N_nonzeros_ = 0;
+		current_N_rows_ = 0;
+		IA_[0] = IO_;
+	}
+
+	void Jacobian::set_element(int II, int JJ, fp value)
+	{
+		const int posn = find_element(II, JJ);
+		if (posn >= 0)
+			then A_[posn] = value;
+		else
+			insert_element(II, JJ, value);
+	}
+
+	void Jacobian::sum_into_element(int II, int JJ, fp value)
+	{
+		const int posn = find_element(II, JJ);
+		if (posn >= 0)
+			then A_[posn] += value;
+		else
+			insert_element(II, JJ, value);
+	}
+
+	int Jacobian::find_element(int II, int JJ)
+		const
+	{
+		if (II >= current_N_rows_)
+			then return -1; // this row not defined yet
+
+		const int start = IA_[II] - IO_;
+		const int stop = IA_[II + 1] - IO_;
+		for (int posn = start; posn < stop; ++posn)
+		{
+			if (JA_[posn] - IO_ == JJ)
+				then return posn; // found
+		}
+
+		return -1; // not found
+	}
+
+	int Jacobian::insert_element(int II, int JJ, double value)
+	{
+		if (!((II == current_N_rows_ - 1) || (II == current_N_rows_)))
+		{
+			printf(
+				"***** row_sparse_Jacobian::insert_element(II=%d, JJ=%d, value=%g):\n"
+				"        attempt to insert element elsewhere than {last row, last row+1}!\n"
+				"        N_rows_=%d   current_N_rows_=%d   IO_=%d\n"
+				"        N_nonzeros_=%d   N_nonzeros_allocated_=%d\n",
+				II, JJ, double(value),
+				N_rows_, current_N_rows_, IO_,
+				N_nonzeros_, N_nonzeros_allocated_);
+			abort();
+		}
+
+		// start a new row if necessary
+		if (II == current_N_rows_)
+			then
+			{
+				assert(current_N_rows_ < N_rows_);
+				IA_[current_N_rows_ + 1] = IA_[current_N_rows_];
+				++current_N_rows_;
+			}
+
+		// insert into current row
+		assert(II == current_N_rows_ - 1);
+		if (IA_[II + 1] - IO_ >= N_nonzeros_allocated_)
+			then grow_arrays();
+		const int posn = IA_[II + 1] - IO_;
+		assert(posn < N_nonzeros_allocated_);
+		JA_[posn] = JJ + IO_;
+		A_[posn] = value;
+		++IA_[II + 1];
+		++N_nonzeros_;
+
+		return posn;
+	}
+
+	void Jacobian::grow_arrays()
+	{
+		N_nonzeros_allocated_ += base_growth_amount + (N_nonzeros_allocated_ >> 1);
+
+		int *const new_JA = new int[N_nonzeros_allocated_];
+		double *const new_A = new double[N_nonzeros_allocated_];
+		for (int posn = 0; posn < N_nonzeros_; ++posn)
+		{
+			new_JA[posn] = JA_[posn];
+			new_A[posn] = A_[posn];
+		}
+		delete[] A_;
+		delete[] JA_;
+		JA_ = new_JA;
+		A_ = new_A;
+	}
+
+	int compare_matrix_elements(const void *x, const void *y)
+	{
+		const struct matrix_element *const px = static_cast<const struct matrix_element *>(x);
+		const struct matrix_element *const py = static_cast<const struct matrix_element *>(y);
+
+		return px->JA - py->JA;
+	}
+
+	void Jacobian::sort_each_row_into_column_order()
+	{
+		// buffer must be big enough to hold the largest row
+		int max_N_in_row = 0;
+		{
+			for (int II = 0; II < N_rows_; ++II)
+			{
+				max_N_in_row = max(max_N_in_row, IA_[II + 1] - IA_[II]);
+			}
+		}
+
+		// contiguous buffer for sorting
+		struct matrix_element *const buffer = new struct matrix_element[max_N_in_row];
+
+		{
+			for (int II = 0; II < N_rows_; ++II)
+			{
+				const int N_in_row = IA_[II + 1] - IA_[II];
+
+				// copy this row's JA_[] and A_[] values to the buffer
+				const int start = IA_[II] - IO_;
+				for (int p = 0; p < N_in_row; ++p)
+				{
+					const int posn = start + p;
+					buffer[p].JA = JA_[posn];
+					buffer[p].A = A_[posn];
+				}
+
+				// sort the buffer
+				qsort(static_cast<void *>(buffer), N_in_row, sizeof(buffer[0]),
+					  &compare_matrix_elements);
+
+				// copy the buffer values back to this row's JA_[] and A_[]
+				for (int p = 0; p < N_in_row; ++p)
+				{
+					const int posn = start + p;
+					JA_[posn] = buffer[p].JA;
+					A_[posn] = buffer[p].A;
+				}
+			}
+		}
+
+		delete[] buffer;
+	}
+
+	double Jacobian::solve_linear_system(int rhs_gfn, int x_gfn, bool print_msg_flag)
+	{
+		assert(IO_ == Fortran_index_origin);
+		assert(current_N_rows_ == N_rows_);
+
+		if (itemp_ == NULL)
+			then
+			{
+				itemp_ = new int[3 * N_rows_ + 3 * N_nonzeros_ + 2];
+				rtemp_ = new double[4 * N_rows_ + N_nonzeros_];
+			}
+
+		// initial guess = all zeros
+		double *x = ps_.gridfn_data(x_gfn);
+		for (int II = 0; II < N_rows_; ++II)
+		{
+			x[II] = 0.0;
+		}
+
+		const int N = N_rows_;
+		const double *rhs = ps_.gridfn_data(rhs_gfn);
+		const double eps = 1e-10;
+		const int max_iterations = N_rows_;
+		int istatus;
+
+		// the actual linear solution
+		f_ilucg(N,
+				IA_, JA_, A_,
+				rhs, x,
+				itemp_, rtemp_,
+				eps, max_iterations,
+				istatus);
+
+		if (istatus < 0)
+		{
+			printf(
+				"***** row_sparse_Jacobian__ILUCG::solve_linear_system(rhs_gfn=%d, x_gfn=%d):\n"
+				"        error return from [sd]ilucg() routine!\n"
+				"        istatus=%d < 0 ==> bad matrix structure, eg. zero diagonal element!\n",
+				rhs_gfn, x_gfn,
+				int(istatus));
+			abort();
+		}
+
+		return -1.0;
+	}
+
+} // namespace AHFinderDirect

AMSS_NCKU_source/Jacobian.h

+#ifndef AHFINDERDIRECT__JACOBIAN_HH
+#define AHFINDERDIRECT__JACOBIAN_HH
+
+namespace AHFinderDirect
+{
+	class Jacobian
+	{
+	public:
+		// basic meta-info
+		patch_system &my_patch_system() const { return ps_; }
+		int N_rows() const { return N_rows_; }
+
+		// convert (patch,irho,isigma) <--> row/column index
+		int II_of_patch_irho_isigma(const patch &p, int irho, int isigma)
+			const
+		{
+			return ps_.gpn_of_patch_irho_isigma(p, irho, isigma);
+		}
+		const patch &patch_irho_isigma_of_II(int II, int &irho, int &isigma)
+			const
+		{
+			return ps_.patch_irho_isigma_of_gpn(II, irho, isigma);
+		}
+
+		double element(int II, int JJ) const;
+
+		// is the matrix element (II,JJ) stored explicitly?
+		bool is_explicitly_stored(int II, int JJ) const
+		{
+			return find_element(II, JJ) > 0;
+		}
+
+		int IO() const { return IO_; }
+		enum
+		{
+			C_index_origin = 0,
+			Fortran_index_origin = 1
+		};
+
+		void zero_matrix();
+
+		void set_element(int II, int JJ, fp value);
+
+		void sum_into_element(int II, int JJ, fp value);
+
+		int find_element(int II, int JJ) const;
+
+		int insert_element(int II, int JJ, fp value);
+
+		void grow_arrays();
+
+		enum
+		{
+			base_growth_amount = 1000
+		};
+
+		void sort_each_row_into_column_order();
+
+		double solve_linear_system(int rhs_gfn, int x_gfn,
+								   bool print_msg_flag);
+
+	public:
+		Jacobian(patch_system &ps);
+		~Jacobian();
+
+	protected:
+		patch_system &ps_;
+		int N_rows_;
+
+		int IO_;
+
+		int N_nonzeros_;
+		int current_N_rows_;
+
+		int N_nonzeros_allocated_;
+
+		int *IA_;
+
+		int *JA_;
+
+		double *A_;
+
+		int *itemp_;
+		double *rtemp_;
+	};
+
+	//******************************************************************************
+
+} // namespace AHFinderDirect
+#endif /* AHFINDERDIRECT__JACOBIAN_HH */

AMSS_NCKU_source/MPatch.h

+
+#ifndef PATCH_H
+#define PATCH_H
+
+#include <mpi.h>
+#include "MyList.h"
+#include "Block.h"
+#include "var.h"
+#include "macrodef.h" //need dim here; Vertex or Cell; ghost_width
+
+class Patch
+{
+
+public:
+   int lev;
+   int shape[dim];
+   double bbox[2 * dim]; // this bbox includes buffer points
+   MyList<Block> *blb, *ble;
+   int lli[dim], uui[dim]; // denote the buffer points on each boundary
+
+public:
+   Patch() {};
+   Patch(int DIM, int *shapei, double *bboxi, int levi, bool buflog, int Symmetry);
+
+   ~Patch();
+
+   void checkPatch(bool buflog);
+   void checkPatch(bool buflog, const int out_rank);
+   void checkBlock();
+   void Interp_Points(MyList<var> *VarList,
+                      int NN, double **XX,
+                      double *Shellf, int Symmetry);
+   bool Interp_ONE_Point(MyList<var> *VarList, double *XX,
+                         double *Shellf, int Symmetry);
+   double getdX(int dir);
+
+   void Find_Maximum(MyList<var> *VarList, double *XX,
+                     double *Shellf);
+
+   bool Find_Point(double *XX);
+
+   void Interp_Points(MyList<var> *VarList,
+                      int NN, double **XX,
+                      double *Shellf, int Symmetry, MPI_Comm Comm_here);
+   bool Interp_ONE_Point(MyList<var> *VarList, double *XX,
+                         double *Shellf, int Symmetry, MPI_Comm Comm_here);
+   void Find_Maximum(MyList<var> *VarList, double *XX,
+                     double *Shellf, MPI_Comm Comm_here);
+};
+
+#endif /* PATCH_H */

AMSS_NCKU_source/MyList.h

+
+#ifndef MYLIST_H
+#define MYLIST_H
+
+// Note: There is never an implementation file (*.C) for a template class
+
+template <class T>
+class MyList
+{
+
+public:
+  MyList *next;
+  T *data;
+
+public:
+  MyList();
+  MyList(T *p);
+  ~MyList();
+  void insert(T *p);
+  void clearList();
+  void destroyList();
+  void catList(MyList<T> *p);
+  void CloneList(MyList<T> *p);
+};
+
+template <class T>
+MyList<T>::MyList()
+{
+  data = 0;
+  next = 0;
+}
+template <class T>
+MyList<T>::MyList(T *p)
+{
+  data = p;
+  next = 0;
+}
+
+template <class T>
+MyList<T>::~MyList()
+{
+}
+template <class T>
+void MyList<T>::insert(T *p)
+{
+  MyList *ct = this;
+  if (data == 0)
+  {
+    data = p;
+  }
+  else
+  {
+    while (ct->next)
+    {
+      ct = ct->next;
+    }
+    ct->next = new MyList(p);
+    ct = ct->next;
+    ct->next = 0;
+  }
+}
+template <class T>
+void MyList<T>::clearList()
+{
+  MyList *ct = this, *n;
+  while (ct)
+  {
+    n = ct->next;
+    delete ct;
+    ct = n;
+  }
+}
+template <class T>
+void MyList<T>::destroyList()
+{
+  MyList *ct = this, *n;
+  while (ct)
+  {
+    n = ct->next;
+    delete ct->data;
+    delete ct;
+    ct = n;
+  }
+}
+template <class T>
+void MyList<T>::catList(MyList<T> *p)
+{
+  MyList *ct = this;
+  while (ct->next)
+  {
+    ct = ct->next;
+  }
+  ct->next = p;
+}
+template <class T>
+void MyList<T>::CloneList(MyList<T> *p)
+{
+  MyList *ct = this;
+  p = 0;
+  while (ct)
+  {
+    if (!p)
+      p = new MyList<T>(ct->data);
+    else
+      p->insert(ct->data);
+    ct = ct->next;
+  }
+}
+#endif /* MyList_H */

AMSS_NCKU_source/NullEvol.h

+
+#ifndef NULLEVOL_H
+#define NULLEVOL_H
+
+#ifdef fortran1
+#define f_setup_dyad setup_dyad
+#define f_eth_derivs eth_derivs
+#define f_eth_dderivs eth_dderivs
+#define f_fill_symmetric_boundarybuffer fill_symmetric_boundarybuffer
+#define f_fill_symmetric_boundarybuffer2 fill_symmetric_boundarybuffer2
+#define f_calculate_K calculate_k
+#define f_NullEvol_beta nullevol_beta
+#define f_NullEvol_Q nullevol_q
+#define f_NullEvol_U nullevol_u
+#define f_NullEvol_W nullevol_w
+#define f_NullEvol_Theta nullevol_theta
+#define f_NullEvol_Theta_givenx nullevol_theta_givenx
+#define f_Eq_Theta eq_theta
+#define f_Eq_Theta_2 eq_theta_2
+#define f_NullEvol_g01 nullevol_g01
+#define f_NullEvol_pg0A nullevol_pg0a
+#define f_NullEvol_Theta2 nullevol_theta2
+#define f_NullEvol_Thetag00 nullevol_thetag00
+#endif
+#ifdef fortran2
+#define f_setup_dyad SETUP_DYAD
+#define f_eth_derivs ETH_DERIVS
+#define f_eth_dderivs ETH_DDERIVS
+#define f_fill_symmetric_boundarybuffer FILL_SYMMETRIC_BOUNDARYBUFFER
+#define f_fill_symmetric_boundarybuffer2 FILL_SYMMETRIC_BOUNDARYBUFFER2
+#define f_calculate_K CALCULATE_K
+#define f_NullEvol_beta NULLEVOL_BETA
+#define f_NullEvol_Q NULLEVOL_Q
+#define f_NullEvol_U NULLEVOL_U
+#define f_NullEvol_W NULLEVOL_W
+#define f_NullEvol_Theta NULLEVOL_THETA
+#define f_NullEvol_Theta_givenx NULLEVOL_THETA_GIVENX
+#define f_Eq_Theta EQ_THETA
+#define f_Eq_Theta_2 EQ_THETA_2
+#define f_NullEvol_g01 NULLEVOL_G01
+#define f_NullEvol_pg0A NULLEVOL_PG0A
+#define f_NullEvol_Theta2 NULLEVOL_THETA2
+#define f_NullEvol_Thetag00 NULLEVOL_THETAG00
+#endif
+#ifdef fortran3
+#define f_setup_dyad setup_dyad_
+#define f_eth_derivs eth_derivs_
+#define f_eth_dderivs eth_dderivs_
+#define f_fill_symmetric_boundarybuffer fill_symmetric_boundarybuffer_
+#define f_fill_symmetric_boundarybuffer2 fill_symmetric_boundarybuffer2_
+#define f_calculate_K calculate_k_
+#define f_NullEvol_beta nullevol_beta_
+#define f_NullEvol_Q nullevol_q_
+#define f_NullEvol_U nullevol_u_
+#define f_NullEvol_W nullevol_w_
+#define f_NullEvol_Theta nullevol_theta_
+#define f_NullEvol_Theta_givenx nullevol_theta_givenx_
+#define f_Eq_Theta eq_theta_
+#define f_Eq_Theta_2 eq_theta_2_
+#define f_NullEvol_g01 nullevol_g01_
+#define f_NullEvol_pg0A nullevol_pg0a_
+#define f_NullEvol_Theta2 nullevol_theta2_
+#define f_NullEvol_Thetag00 nullevol_thetag00_
+#endif
+
+extern "C"
+{
+	void f_setup_dyad(int *, double *, double *, double *,
+					  double *, double *, double *, double *,
+					  double *, double *, double *, double *,
+					  double *, double *,
+					  double *, double *, double *, double *,
+					  double *, double *, double *, double *,
+					  double *, double *, double *, double *,
+					  double *, double *, double *,
+					  int &, double &);
+}
+
+extern "C"
+{
+	void f_eth_derivs(int *, double *, double *,
+					  double *, double *,
+					  double *, double *,
+					  int &, int &,
+					  double *, double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	void f_eth_dderivs(int *, double *, double *,
+					   double *, double *,
+					   double *, double *,
+					   int &, int &, int &,
+					   double *, double *, double *, double *, double *, double *,
+					   double *, double *, double *, double *,
+					   double *, double *, double *, double *,
+					   double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	void f_fill_symmetric_boundarybuffer(int *, double *, double *, double *,
+										 double &, double &,
+										 double *, double *, double *, double *, double *, double *, double *, double *,
+										 double *, double *, int &, int &, int &);
+}
+
+extern "C"
+{
+	void f_fill_symmetric_boundarybuffer2(int *, double *, double *, double *,
+										  double &, double &,
+										  double *, int &, int &, double *);
+}
+
+extern "C"
+{
+	void f_calculate_K(int *, double *, double *, double *,
+					   double *, double *,
+					   double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_NullEvol_beta(int *, double *, double *, double *,
+						double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_NullEvol_Q(int *, double *, double *, double *,
+					 double *, double *, double *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_NullEvol_U(int *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *, double *,
+					 double *, double *, double &);
+}
+
+extern "C"
+{
+	int f_NullEvol_W(int *, double *, double *, double *,
+					 double *, double *, double *, double *, double *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *, double *, double *, double &,
+					 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_NullEvol_Theta(int *, double *, double *, double *,
+						 double *, double *, double *, double *, double *, double *, double *,
+						 double *, double *, double *, double *, double *, double *, double *,
+						 double *, double *, double *, double *,
+						 double &,
+						 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_NullEvol_Theta_givenx(int *, double *, double *, double *,
+								double *, double *, double *, double *, double *, double *, double *,
+								double *, double *, double *, double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double &,
+								double *, double *, double *, double *, double *, double *, double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double &, int &);
+}
+
+extern "C"
+{
+	int f_Eq_Theta(int *, double *, double *, double *,
+				   double *, double *, double *, double *, double *, double *, double *,
+				   double *, double *, double *, double *, double *, double *, double *, double &,
+				   double *, double *, double *, double *, double *, double *, double *, double *, double *, double *);
+}
+
+extern "C"
+{
+	int f_Eq_Theta_2(int *, double *, double *, double *,
+					 double *, double *, double *, double *, double *, double *, double *,
+					 double *, double *, double *, double *, double *, double *, double *, double &,
+					 double *, double *, double *, double *, double *, double *, double *, double *, double *, double *,
+					 double &, int &);
+}
+
+extern "C"
+{
+	int f_NullEvol_g01(int *, double *, double *, double *,
+					   double *, double *, double *, double *,
+					   double &);
+}
+
+extern "C"
+{
+	int f_NullEvol_pg0A(int *, double *, double *, double *,
+						double *, double *, double *, double *,
+						double *, double *, double *, double *,
+						double &);
+}
+
+extern "C"
+{
+	int f_NullEvol_Theta2(int *, double *, double *, double *,
+						  double *, double *, double *, double *, double *, double *, double *, double *, double *,
+						  double *, double *, double *,
+						  double &);
+}
+
+extern "C"
+{
+	int f_NullEvol_Thetag00(int *, double *, double *, double *,
+							double *, double *, double *, double *, double *, double *, double *, double *, double *,
+							double *, double *, double *,
+							double &);
+}
+#endif /* NULLEVOL_H */

AMSS_NCKU_source/NullNews.h

+
+#ifndef NULLNEWS_H
+#define NULLNEWS_H
+
+#ifdef fortran1
+#define f_drive_null_news drive_null_news
+#define f_get_null_news2 get_null_news2
+#define f_drive_null_news_diff drive_null_news_diff
+#define f_omega_rhs omega_rhs
+#define f_get_exact_omega get_exact_omega
+#define f_get_omega_and_dtomega_pre get_omega_and_dtomega_pre
+#define f_get_omega_and_dtomega_LN get_omega_and_dtomega_ln
+#define f_get_dtomega get_dtomega
+#endif
+#ifdef fortran2
+#define f_drive_null_news DRIVE_NULL_NEWS
+#define f_get_null_news2 GET_NULL_NEWS2
+#define f_drive_null_news_diff DRIVE_NULL_NEWS_DIFF
+#define f_omega_rhs OMEGA_RHS
+#define f_get_exact_omega GET_EXACT_OMEGA
+#define f_get_omega_and_dtomega_pre GET_OMEGA_AND_DTOMEGA_PRE
+#define f_get_omega_and_dtomega_LN GET_OMEGA_AND_DTOMEGA_LN
+#define f_get_dtomega GET_DTOMEGA
+#endif
+#ifdef fortran3
+#define f_drive_null_news drive_null_news_
+#define f_get_null_news2 get_null_news2_
+#define f_drive_null_news_diff drive_null_news_diff_
+#define f_omega_rhs omega_rhs_
+#define f_get_exact_omega get_exact_omega_
+#define f_get_omega_and_dtomega_pre get_omega_and_dtomega_pre_
+#define f_get_omega_and_dtomega_LN get_omega_and_dtomega_ln_
+#define f_get_dtomega get_dtomega_
+#endif
+
+extern "C"
+{
+	void f_drive_null_news(int *, double *, double *, double *,
+						   double *, double *, double *, double *, double *, double *, double *, double *,
+						   double *, double *, double *, double *,
+						   double *, double *, double *, double *,
+						   double *, double *,
+						   double *, double *, double *, double *,
+						   double *, double *, double *, double *,
+						   double *, double *, double *, double *,
+						   double *, double *, double &, int &);
+}
+
+extern "C"
+{
+	void f_drive_null_news_diff(int *, double *, double *, double *,
+								double *, double *, double *, double *, double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double *, double *,
+								double *, double *, double &, int &, double &);
+}
+
+extern "C"
+{
+	void f_omega_rhs(int *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *, double *, double *,
+					 double *, double *);
+}
+
+extern "C"
+{
+	void f_get_exact_omega(int *, double *, double *, double *,
+						   double *,
+						   int &, double &, double &);
+}
+
+extern "C"
+{
+	void f_get_null_news2(int *, double *, double *, double *,
+						  double *, double *,
+						  double *, double *, double *, double *,
+						  double *, double *, double *,
+						  double *, double *, double *,
+						  double *, double *, double &, int &);
+}
+
+extern "C"
+{
+	void f_get_omega_and_dtomega_pre(int *, double *, double *, double *,
+									 double *, double *, double *,
+									 double *, double *, double &);
+}
+
+extern "C"
+{
+	void f_get_dtomega(int *, double *, double *, double *,
+					   double *, double *, double *,
+					   double *, double *, double &);
+}
+
+extern "C"
+{
+	void f_get_omega_and_dtomega_LN(double &, int *, double *, double *, double *,
+									double *, double *, double &, int &);
+}
+#endif /* NULLNEWS_H */

AMSS_NCKU_source/Parallel_bam.h

+
+#ifndef PARALLEL_BAM_H
+#define PARALLEL_BAM_H
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cstdlib>
+#include <cstdio>
+#include <string>
+#include <cmath>
+#include <new>
+using namespace std;
+
+#include "var.h"
+#include "MPatch.h"
+#include "Block.h"
+#include "MyList.h"
+#include "macrodef.h"
+namespace Parallel
+{
+  struct pointstru_bam
+  {
+    double pox[dim]; // cordinate
+    Block *Bgs;      // interplate from
+    Block *Bgd;      // interplate for
+    double *coef;    // interpolation coefficients
+    int sind[dim];   // interpolation starting array index
+  };
+  void destroypsuList_bam(MyList<pointstru_bam> *ct);
+  void OutBdLow2Hi_bam(MyList<Patch> *PLc, MyList<Patch> *PLf,
+                       MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
+                       int Symmetry);
+  void OutBdLow2Hi_bam(MyList<Patch> *PLc, MyList<Patch> *PLf,
+                       MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
+                       MyList<Parallel::pointstru_bam> *bdsul, int Symmetry);
+  void Constr_pointstr_OutBdLow2Hi(MyList<Patch> *PLf, MyList<Patch> *PLc,
+                                   MyList<Parallel::pointstru_bam> *&bdsul);
+  void Restrict_bam(MyList<Patch> *PLc, MyList<Patch> *PLf,
+                    MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
+                    int Symmetry);
+  void Restrict_bam(MyList<Patch> *PLc, MyList<Patch> *PLf,
+                    MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /* target */,
+                    MyList<Parallel::pointstru_bam> *rsul, int Symmetry);
+  void Constr_pointstr_Restrict(MyList<Patch> *PLf, MyList<Patch> *PLc,
+                                MyList<Parallel::pointstru_bam> *&rsul);
+  void intertransfer(MyList<Parallel::pointstru_bam> *&sul,
+                     MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
+                     int Symmetry);
+  int interdata_packer(double *data, MyList<Parallel::pointstru_bam> *sul, int myrank, int node, int dir,
+                       MyList<var> *VarLists /* source */, MyList<var> *VarListd /* target */, int Symmetry);
+}
+#endif /*PARALLEL_BAM_H */

AMSS_NCKU_source/Z4c_class.h

+
+#ifndef Z4c_CLASS_H
+#define Z4c_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 "cgh.h"
+#include "ShellPatch.h"
+#include "misc.h"
+#include "var.h"
+#include "MyList.h"
+#include "monitor.h"
+#include "surface_integral.h"
+
+#include "macrodef.h"
+
+#ifdef USE_GPU
+#include "bssn_gpu_class.h"
+#else
+#include "bssn_class.h"
+#endif
+
+class Z4c_class : public bssn_class
+{
+public:
+     Z4c_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);
+     ~Z4c_class();
+
+     void Initialize();
+     void Check_extrop();
+     // Since we have set zero to variables at very begining
+     // we can neglect TZ for initial data setting
+     void Step(int lev, int YN);
+     void Interp_Constraint();
+     void Constraint_Out();
+     void Compute_Constraint();
+
+protected:
+     var *TZo;
+     var *TZ0;
+     var *TZ;
+     var *TZ1;
+     var *TZ_rhs;
+};
+#endif /* Z4c_CLASS_H */

AMSS_NCKU_source/bssn2adm.f90

+
+!-------------------------------------------------------------------------------!
+! 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