Initial commit

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

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

+//$Id: Block.C,v 1.4 2013/05/18 09:17:31 zjcao Exp $
+#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

+//$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 */

AMSS_NCKU_source/DataCT.C

+//$Id: DataCT.C,v 1.2 2012/04/03 10:49:42 zjcao Exp $
+//-----------------------------------------------------------------------
+// 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

+!$Id: FFT.f90,v 1.1 2012/07/23 05:49:10 zjcao Exp $
+
+#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

+//$Id: MPatch.h,v 1.3 2013/03/30 02:31:17 zjcao Exp $
+#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

+//$Id: MyList.h,v 1.3 2012/06/10 03:09:39 zjcao Exp $
+#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

+//$Id: NullEvol.h,v 1.5 2013/08/20 11:49:04 zjcao Exp $
+#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

+//$Id: NullNews.h,v 1.4 2013/08/20 11:49:04 zjcao Exp $
+#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/NullShellPatch.h

+//$Id: NullShellPatch.h,v 1.5 2013/01/26 11:23:09 zjcao Exp $
+#ifndef NULLSHELLPATCH_H
+#define NULLSHELLPATCH_H
+
+#ifdef newc
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cstdlib>
+#include <cstdio>
+#include <string>
+#include <cmath>
+#include <complex>
+using namespace std;
+#else
+#include <iostream.h>
+#include <iomanip.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+#endif
+
+#include <mpi.h>
+#include "MyList.h"
+#include "Block.h"
+#include "Parallel.h"
+#include "ShellPatch.h"
+#include "var.h"
+#include "macrodef.h" //need dim here; Vertex or Cell; ghost_width
+
+#if (dim != 3)
+#error NullShellPatch only supports 3 dimensional stuff yet
+#endif
+
+class xp_npatch : public ss_patch
+{
+public:
+   xp_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 2; };
+};
+
+class xm_npatch : public ss_patch
+{
+public:
+   xm_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 3; };
+};
+class yp_npatch : public ss_patch
+{
+public:
+   yp_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 4; };
+};
+
+class ym_npatch : public ss_patch
+{
+public:
+   ym_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 5; };
+};
+class zp_npatch : public ss_patch
+{
+public:
+   zp_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 0; };
+};
+
+class zm_npatch : public ss_patch
+{
+public:
+   zm_npatch(int ingfsi, int fngfsi, int *shapei, double *bboxi, int myranki) : ss_patch(ingfsi, fngfsi, shapei, bboxi, myranki) { sst = 1; };
+};
+
+class NullShellPatch
+{
+
+public:
+   struct pointstru
+   {
+      double gpox[dim]; // global cordinate
+      double lpox[dim]; // local cordinate
+      Block *Bg;
+      int ssst; //-1: cardisian, others as sst of ss_patch source sst
+      int tsst; //-1: cardisian, others as sst of ss_patch target sst
+      double *coef;
+      int *sind;
+      int dumyd;            // the dimension which has common lines, only useful in interdata_packer
+      complex<double> swtf; // exp(i gamma) of Eq.(26) of CQG 24 S327
+   };
+
+   var *FXZEO;
+   var *gx, *gy, *gz;
+   // we always assume the number of VarList = 2* the number of Varwt
+   // so VarList must apear with pairs, either components of complex number or a fake pair
+   var *beta, *W;
+   var *Rnu, *Inu, *Rk, *Ik, *RB, *IB;
+   var *RQ, *IQ, *RU, *IU, *RTheta, *ITheta;
+   var *KK, *HKK, *KKx, *HKKx;
+   var *RJo, *IJo, *omegao;
+   var *RJ0, *IJ0, *omega0;
+   var *RJ, *IJ, *omega;
+   var *RJ1, *IJ1, *omega1;
+   var *RJ_rhs, *IJ_rhs, *omega_rhs;
+
+   var *quR1, *quR2, *quI1, *quI2;
+   var *qlR1, *qlR2, *qlI1, *qlI2;
+   var *gR, *gI;
+   var *dquR1, *dquR2, *dquI1, *dquI2;
+   var *bdquR1, *bdquR2, *bdquI1, *bdquI2;
+   var *dgR, *dgI;
+   var *bdgR, *bdgI;
+
+   var *RNews, *INews;
+
+   MyList<var> *StateList, *SynchList_pre, *SynchList_cor, *RHSList;
+   MyList<var> *OldStateList, *DumpList, *CheckList;
+
+   MyList<var> *betaList, *QUList, *WTheList, *TheList, *JrhsList, *J1List;
+   int betawt[1], QUwt[2], WThewt[2];
+
+   int myrank;
+   int shape[dim]; // for (rho, sigma, X), for rho and sigma means number of points for every pi/2
+   double Rmin, xmin, xmax;
+   int Symmetry;
+   int ingfs, fngfs;
+
+   MyList<ss_patch> *PatL;
+
+   MyList<pointstru> **ss_src, **ss_dst;
+   MyList<pointstru> **cs_src, **cs_dst;
+
+public:
+   NullShellPatch(int *shapei, double Rmini, double xmini, double xmaxi, int Symmetry, int myranki);
+
+   ~NullShellPatch();
+
+   void destroypsuList(MyList<pointstru> *ct);
+   void fill_symmetric_boundarybuffer(MyList<var> *VarList, int *Varwt);
+   MyList<Block> *compose_sh(int cpusize);
+   int getdumydimension(int acsst, int posst);
+   void Setup_dyad();
+   void Setup_Initial_Data(bool checkrun, double PhysTime);
+   void eth_derivs(var *Rv, var *Iv, var *ethRv, var *ethIv, int s, int e);
+   void eth_dderivs(var *Rv, var *Iv, var *ethRv, var *ethIv, int s, int e1, int e2);
+   void getlocalpox_ss(int isst, double ix, double iy, double iz, int &sst, double &lx, double &ly, double &lz);
+   void getlocalpox_fake(double x, double y, double z, int &sst, double &lx, double &ly, double &lz);
+   void getlocalpox(double x, double y, double z, int &sst, double &lx, double &ly, double &lz);
+   void getlocalpoxsst_ss(int isst, double ix, double iy, double iz, int lsst, double &lx, double &ly, double &lz);
+   void getlocalpoxsst(double x, double y, double z, int sst, double &lx, double &ly, double &lz);
+   void getglobalpox(double &x, double &y, double &z, int sst, double lx, double ly, double lz);
+   complex<double> get_swtf(double *pox, int tsst, int ssst);
+   void prolongpointstru(MyList<pointstru> *&psul, MyList<ss_patch> *sPpi, double DH[dim],
+                         MyList<Patch> *Ppi, double CDH[dim], MyList<pointstru> *pss);
+   bool prolongpointstru(MyList<pointstru> *&psul, bool ssyn, int tsst, MyList<ss_patch> *sPp, double DH[dim],
+                         MyList<Patch> *Pp, double CDH[dim], double x, double y, double z, int Symmetry, int rank_in);
+   bool prolongpointstru_ss(MyList<pointstru> *&psul, int tsst, MyList<ss_patch> *sPp, double DH[dim],
+                            MyList<Patch> *Pp, double CDH[dim], double x, double y, double z, int Symmetry, int rank_in);
+   void setupintintstuff(int cpusize, MyList<Patch> *CPatL, int Symmetry);
+   void checkPatch();
+   void checkBlock(int sst);
+   double getdX(int dir);
+   void shellname(char *sn, int i);
+   void Dump_xyz(char *tag, double time, double dT);
+   void Dump_Data(MyList<var> *DumpListi, char *tag, double time, double dT);
+   void intertransfer(MyList<pointstru> **src, MyList<pointstru> **dst,
+                      MyList<var> *VarList1 /* source */, MyList<var> *VarList2 /*target */,
+                      int Symmetry, int *Varwt);
+   int interdata_packer(double *data, MyList<pointstru> *src, MyList<pointstru> *dst, int rank_in, int dir,
+                        MyList<var> *VarLists /* source */, MyList<var> *VarListd /* target */, int Symmetry, int *Varwt);
+   void Synch(MyList<var> *VarList, int Symmetry, int *Varwt);
+   void CS_Inter(MyList<var> *VarList, int Symmetry, int *Varwt);
+   void check_pointstrul(MyList<pointstru> *pp, bool first_only);
+   void check_pointstrul2(MyList<pointstru> *pp, int first_last_only);
+   void matchcheck(MyList<Patch> *CPatL);
+   void Interp_Points(MyList<var> *VarList,
+                      int NN, double **XX, /*input global Cartesian coordinate*/
+                      double *Shellf, int Symmetry);
+   void Interp_Points_2D(MyList<var> *VarList,
+                         int NN, double **XX, /*input global Cartesian coordinate*/
+                         double *Shellf, int Symmetry);
+   void Step(double dT, double PhysTime, monitor *ErrorMonitor);
+   void Null_Boundary(double PhysTime);
+   void HyperSlice(double dT, double PhysTime, monitor *ErrorMonitor, int RK_count);
+   double News_Error_Check(double PhysTime, double dT, bool dp);
+   double Error_Check(double PhysTime, double dT, bool dp);
+   double EqTheta_Check(double PhysTime, double dT, bool dp);
+   void Compute_News(double PhysTime, double dT, bool dp);
+   void Check_News(double PhysTime, double dT, bool dp);
+};
+
+#endif /* NULLSHELLPATCH_H */