Initial commit

.gitignore

+*.c
+*.so
+*.o

README.md

+# pyloc3d
+
+Python bindings for Sarah's `AHloc3d`.
+
+# Instructions
+
+- Clone this repo and the `AHloc3d` repo.
+- Compile `AHloc3d` into a shared library (I think only dependency is `MPI` which is not used
+  by `pyloc3d` but still needs to be installed)
+```sh
+cd ahloc3d
+make lib
+```
+- Generate the glue code for python. This might require you to update the arguments
+of the function call `ffibuilder.set_source` inside `wrap_ahloc3d.py` to conform with
+your installation. In particular for `mpi` you will have to specify the right
+`include_dirs` which you should be able to query from a command line using
+```sh
+pkg-config --cflags mpi
+```
+Then run to following to wrap it:
+```sh
+cd pyloc3d
+./wrap_ahloc3d
+```
+- You need to set `LD_LIBRARY_PATH` so that python can find `libAHloc3d`
+```
+export LD_LIBRARY_PATH="<absolute-path-to-ahloc3d-git-folder>"
+```
+- Test it
+```sh
+./pyloc3d
+```
+
+
+# Hack notes
+
+- If `AHloc3d` changes one needs to rerun the instructions above.

pyloc3d.py

+#!/usr/bin/env python3
+
+from _ahloc_3d import ffi, lib
+from ctypes import create_string_buffer
+
+parfile = "../ahloc3d/test.par"
+cstr_parfile = parfile.encode('utf8')
+lib.setup_output(cstr_parfile)
+lib.ReadParfile(cstr_parfile)
+wus = lib.read_work()

test.par

+###############################################################################
+# AHloc parameter file                                                        #
+###############################################################################
+
+###############################################################################
+# Symmetry mode                                                               #
+###############################################################################
+
+sym.mode = AXISYM
+
+###############################################################################
+# input/output                                                                #
+###############################################################################
+
+test.Kerr = 0
+test.Kerr.a = 0.0
+
+AHloc.datadir = /home/sarah/data/new_ah_test/old_data/output_ah
+AHloc.outdir = old_data_new_ahloc
+
+data.legacy = 1
+
+output.flowstatus = 1
+output.flowstatus.every = 100
+output.flowsurf = 0
+output.flowsurf.every = 100
+
+###############################################################################
+# horizon parameters                                                          #
+###############################################################################
+
+horizon.res.theta = 50
+
+horizon.max.exp = 0.0000000001
+
+###############################################################################
+# Algorithm parameters                                                        #
+###############################################################################
+
+flow.initial = 25.0
+flow.max_exp = 10000
+flow.dtfactor = 0.5
+flow.handoff = 0.1
+newton.eps = 0.0000000001
+
+###############################################################################
+# Search along z-axis                                                         #
+###############################################################################
+
+offset.z.min = 0.0;
+offset.z.max = 3.0;
+offset.z.n = 1;

wrap_ahloc3d.py

+#!/usr/bin/env python3
+
+
+from cffi import FFI
+ffibuilder = FFI()
+
+
+ahloc3d_h = """
+// Symmetry mode
+typedef enum tSym_mode{ FULL3D = 0, AXISYM = 1 } tSym_mode;
+
+// Flag for using data from <=4.0
+typedef enum tGridConvention { CURRENT = 0, LEGACY = 1 } tGridConvention;
+
+// Domain type
+typedef enum tDomainType{ CUBEDSPHERE = 0, CUBEDBALL = 1} tDomainType;
+
+// Grid type
+typedef enum tGridType{
+  CC=0,     ///< cube to cube (deformed cube)
+  CS=1,     ///< cube to sphere
+  SS=2,     ///< sphere to sphere
+  SC=3,     ///< sphere to cube
+  CU=4,     ///< cube
+  DF=5      ///< dual foliation
+} tGridType;
+
+// Grid orientation
+typedef enum tGridOrientation{ XP=0, YP=1, ZP=2,
+                               XM=3, YM=4, ZM=5, CENTER=6 } tGridOrientation;
+
+typedef struct tTimeStepInfo {
+  tDomainType domain;
+  double cu_max;
+  double cs_max;
+  double ss_max;
+  int reflect_x;
+  int reflect_y;
+  int reflect_z;
+  int n_grids;
+  int evolve_step;
+  double time;
+} tTimeStepInfo;
+
+typedef struct tGridInfo {
+  int n_grid;
+  tGridType type;
+  tGridOrientation orientation;
+  int nu;
+  int nv;
+  int nw;
+  double u_min;
+  double u_max;
+  double v_min;
+  double v_max;
+  double w_min;
+  double w_max;
+} tGridInfo;
+
+// Names of Cartesian tensor components
+typedef enum tFields{
+  GXX = 0, GXY = 1, GXZ = 2, GYY = 3, GYZ = 4,  GZZ = 5,
+  KXX = 6, KXY = 7, KXZ = 8, KYY = 9, KYZ = 10, KZZ = 11
+} tFields;
+
+// Names of needed derivatives of metric components
+typedef enum tDerivatives{
+  DGXX_DX = 0,  DGXY_DX = 1,  DGXZ_DX = 2,  DGYY_DX = 3,  DGYZ_DX = 4,  DGZZ_DX = 5,
+  DGXX_DY = 6,  DGXY_DY = 7,  DGXZ_DY = 8,  DGYY_DY = 9,  DGYZ_DY = 10, DGZZ_DY = 11,
+  DGXX_DZ = 12, DGXY_DZ = 13, DGXZ_DZ = 14, DGYY_DZ = 15, DGYZ_DZ = 16, DGZZ_DZ = 17
+} tDerivatives;
+
+// Names of Cartesian tensor components for axisymmetry
+typedef enum tFields2d{
+  GXX2D = 0, GXZ2D = 1, GYY2D = 2,  GZZ2D = 3,
+  KXX2D = 4, KXZ2D = 5, KYY2D = 6,  KZZ2D = 7,
+} tFields2d;
+
+// Names of needed derivatives of metric components for axisymmetry
+typedef enum tDerivatives2d{
+  DGXX_DX2D = 0,  DGXZ_DX2D = 1,
+  DGYY_DX2D = 2,  DGZZ_DX2D = 3,
+  DGXX_DZ2D = 4,  DGXZ_DZ2D = 5,
+  DGYY_DZ2D = 6,  DGZZ_DZ2D = 7,
+} tDerivatives2d;
+
+// Distinguish between different test cases
+typedef enum tTestCase{
+  NONE = 0, KERR = 1, BL = 2
+} tTestCase;
+
+typedef enum tAxis {
+  X_AXIS = 0,
+  Z_AXIS = 1
+} tAxis;
+
+/******************************************************************************
+ * Structures
+ *****************************************************************************/
+
+// Grid structure: Defines a patch in global grid
+
+typedef struct {
+
+  // Grid number
+  int ng;
+
+  // Grid type
+  tGridType type;
+
+  // Grid orientation
+  tGridOrientation orientation;
+
+  //Pointer to data for cartesian variables components
+  double *fields_cart;
+
+  //Pointer to data for derivatives of variables components
+  double *derivatives;
+
+  //Number of grid points in x and y direction and total number of points
+  int nu, nv, nw, nn;
+
+  //Range of u and w coordinates
+  double umax, umin, vmax, vmin, wmax, wmin;
+
+  // x-Boundaries of higher-level coordinate domain the grid belongs to.
+  // Needed for coordinate transformation in CS type grids
+  // double patch_xmin, patch_xmax;
+
+  // Pointer to 3d array for cartesian coordinates
+  double *xcart, *ycart, *zcart;
+
+  // Pointer to 3d array for local u,v,w coordinates
+  double *ucoord, *vcoord, *wcoord;
+
+  // Pointer to 1d array of u,v,w coordinates
+  double *ucoord_1d, *vcoord_1d, *wcoord_1d;
+
+  // Pointer to 1d array of spectral weights in u,v,w directions
+  double *uweigh_1d, *vweigh_1d, *wweigh_1d;
+
+  // Jacobians D(u,v,w)/D(x,y,z)
+  double *J_du_dx, *J_du_dy, *J_du_dz;
+  double *J_dv_dx, *J_dv_dy, *J_dv_dz;
+  double *J_dw_dx, *J_dw_dy, *J_dw_dz;
+
+  // Derivatives matrices in u,v,w coordinates
+  double *DM_u, *DM_v, *DM_w;
+
+} tGrid;
+
+
+// Horizon structure: Holds the function h(theta,phi) that defines the horizon,
+// as well as possible offsets to its position.
+
+typedef struct {
+
+  double *h;
+  double *tvalues;
+  double *pvalues;
+  int nt;
+  int np;
+  double offset_x;
+  double offset_y;
+  double offset_z;
+
+} tSurface;
+
+
+// Work unit structure: A work unit represents a time slice and holds the
+// grid information and the data which is needed to search for a horizon within
+// this time slice
+
+typedef struct {
+  // Unique id of work unit
+  int id;
+
+  // Flag to denote whether associated data is loaded or not
+  int data_loaded;
+
+  // type of Domain
+  tDomainType domain;
+
+  // Array of pointers to the grids
+  tGrid **grids;
+
+  // Number of grids
+  int ngrids;
+
+  // Directory of grid data (for legacy data)
+  char gridDir[256];
+
+  // Filename of associated data file
+  char data_filename[2048];
+
+  // Boundaries of cube, cube-to-sphere and sphere-to-sphere grids.
+  double cu_max, cu_min;
+  double cs_max, cs_min;
+  double ss_max, ss_min;
+
+  // Time and number of bamps evolution_step of the slice
+  double time;
+  int evolve_step;
+
+  tSurface **horizons;
+  int n_horizons;
+
+} tWorkUnit;
+
+
+/******************************************************************************
+ * Global variables. Defined and set in variables.c
+ * Note that these are set by the code and not in the par file
+ *****************************************************************************/
+
+// MPI rank
+extern int mpi_rank;
+
+// Number of Cartesian variables
+extern int n_fields_cart;
+
+// Number of variables which should be derivated. So the total number
+// Derivatives is 3 times this number because of x, y and z direction,
+// or 2 times if we're in axisymmetry
+extern int n_derivatives;
+
+// Array of pointers to all work units
+extern tWorkUnit **work;
+
+// Number of work units
+extern int n_work;
+
+// Symmetry flag needed for octant
+extern int odd_fields_x[], odd_fields_y[], odd_fields_z[];
+
+extern int odd_derivs_x[], odd_derivs_y[], odd_derivs_z[];
+
+extern int odd_fields_x_axisym[];
+extern int odd_fields_z_axisym[];
+
+extern int odd_derivs_x_axisym[];
+extern int odd_derivs_z_axisym[];
+
+// Symmetry flags for reflection symmetry. All == 1 is equivalent to octant symmetry
+extern int reflect_x;
+extern int reflect_y;
+extern int reflect_z;
+
+// Name of axes for status printing (adressable by tAxis)
+extern const char* axis_name[];
+
+// Whether the surface only spans the octant
+//extern int oct_surf; //obsolete?
+
+/******************************************************************************
+ * Global parameters. These can be set and changed by par file
+ *****************************************************************************/
+// Symmetry mode
+extern tSym_mode sym_mode;
+extern int search_x_axis;
+
+// Options for backwards compatibility
+extern tGridConvention gridConvention;
+extern int legacy_data;
+
+
+// Location of data directory
+extern char datadir[1024];
+
+// Output directory and output options
+extern char outdir[1024];
+extern int out_flow_status;
+extern int out_flow_status_int;
+extern int out_flow_log;
+extern int out_flow_log_int;
+extern int out_flow_surf;
+extern int out_flow_surf_int;
+
+// Horizon finder resolution
+extern int hor_res_theta;
+extern int hor_res_phi;
+extern int flow_res_theta;
+extern int flow_res_phi;
+
+// Search area
+extern double offset_min;
+extern double offset_max;
+extern int offset_n;
+
+// Horizon finding parameters
+extern double flow_dtfactor;
+extern double newton_eps;
+extern double final_max_exp;
+extern double flow_handoff;
+extern double flow_initial;
+extern double recenter_threshold;
+extern int check_bifurcation_int;
+extern double minimum_radius;
+extern int check_minimum_radius_int;
+extern double flow_max_exp;
+extern int recenter;
+
+// Testing parameters
+// extern int testcase;
+// extern int useSS;
+// extern int useKerr;
+// extern double Kerr_a;
+// extern int useBL;
+
+/******************************************************************************
+ * Function headers (Explaining comments to function above function header in
+ * source files)
+ *****************************************************************************/
+
+// Ahloc.c ////////////////////////////////////////////////////////////////////
+
+// Nothing to see here, move along.
+
+// data.c /////////////////////////////////////////////////////////////////////
+
+int deserialize_read_int(unsigned char *buffer, int *x);
+
+int deserialize_read_double(unsigned char *buffer, double *x);
+
+int deserialize_read_timestepinfo(unsigned char *buffer, tTimeStepInfo *info);
+
+int deserialize_read_gridinfo(unsigned char *buffer, tGridInfo *info);
+
+void load_data(tWorkUnit *wu);
+
+// derivative.c ///////////////////////////////////////////////////////////////
+
+double *ChebyDerivMatrix(int n);
+
+void SetDerivativeMatrices(tWorkUnit *wu);
+
+void PrintDeriativeMatrices(int nw);
+
+void ComputeDerivatives(tGrid *g);
+
+void ComputeDerivatives_axisym(tGrid *g);
+
+void mm2d(double *A, double *M, double* B, int n1, int n2, int dim);
+
+void mm3d(double *A, double *M, double *B, int n1, int n2, int n3, int dim);
+
+void matrix_sym(double *Mnn, double *Mp, double *Mm, int nn);
+
+// fields.c ///////////////////////////////////////////////////////////////////
+
+void load_data_and_derivatives(tWorkUnit *wu);
+
+double *getField(int var, double *f, int nn);
+
+void ReadFieldsCart_legacy(tWorkUnit *wu, tGrid *g);
+
+void FreeFieldsAndDerivatives(tWorkUnit *wu);
+
+// grids.c ////////////////////////////////////////////////////////////////////
+
+int SetupGrids_legacy(tWorkUnit *wu);
+
+tGrid **ReadAllocInitGrids(tWorkUnit *wu);
+
+void set_UVW_coords(tWorkUnit *wu, tGrid *g);
+
+void set_UW_coords(tWorkUnit *wu, tGrid *g);
+
+void set_cart_jac(tWorkUnit *wu, tGrid *g);
+
+void set_cart_jac_axisym(tWorkUnit *wu, tGrid *g);
+
+void PrintCartJac(int nw);
+
+void PrintGrids(int nw);
+
+void PrintUWCoords(int nw);
+
+void FreeGrids(tWorkUnit *wu);
+
+// horizon.c //////////////////////////////////////////////////////////////////
+
+int EvaluateExpansion(double *Exp, double *DFNorm, tWorkUnit *wu, tSurface *surf);
+
+int FindHorizon(tSurface ***horizons, tWorkUnit *wu, tSurface *prev, double offx, double offy, double offz, double r);
+
+int HorizonFlow(tWorkUnit *wu, tSurface *surf, double eps);
+
+int HorizonNewtonSolve(tWorkUnit *wu, tSurface *surf, double eps);
+
+tSurface* RecenterMidpoint(tWorkUnit *wu, tSurface *prev, tSurface *hor, double r);
+
+int timestepICN(tWorkUnit *wu, tSurface *surf, double dt);
+
+void SplitHorizon(tWorkUnit *wu, tSurface *old, double *hMin, double *tMin, double *pMin);
+
+double getAHMass(tWorkUnit *wu, tSurface *hor);
+
+// horizon_axisym.c ///////////////////////////////////////////////////////////
+
+int FindHorizon_axisym(tSurface ***horizons, tWorkUnit *wu, tSurface *prev, double x, double z, double r);
+
+int SweepInterval_axisym(tSurface ***horizons, tWorkUnit *wu, double min, double max, int n, tAxis axis);
+
+tSurface* RecenterMidpoint_axisym(tWorkUnit *wu, tSurface *prev, tSurface *hor, double r);
+
+double HorizonMin_x(tSurface *hor);
+
+double HorizonMax_x(tSurface *hor);
+
+double HorizonMin_z(tSurface *hor);
+
+double HorizonMax_z(tSurface *hor);
+
+int EvaluateExpansion_axisym(double *Exp, double *DFNorm, tWorkUnit *wu, tSurface *surf);
+
+int timestepICN_axisym(tWorkUnit *wu, tSurface *surf, double dt);
+
+int HorizonFlow_axisym(tWorkUnit *wu, tSurface *surf, double eps);
+
+int HorizonNewtonSolve_axisym(tWorkUnit *wu, tSurface *surf, double eps);
+
+double getAHMass_axisym(tWorkUnit *wu, tSurface *hor);
+
+double findSplit_axisym(tWorkUnit *wu, tSurface *surf);
+
+// interpolate.c //////////////////////////////////////////////////////////////
+
+typedef tGrid* (*tFindPointInGrids_axisym)(tWorkUnit *wu, double r_b, double theta_b,
+                          double *uc, double *wc);
+
+extern tFindPointInGrids_axisym FindPointInGrids_axisym;
+
+tGrid *FindPointInGrids_axisym_cubedball(tWorkUnit *wu, double r_b, double theta_b,
+                          double *uc, double *wc);
+
+tGrid *FindPointInGrids_axisym_cubedsphere(tWorkUnit *wu, double r_b, double theta_b,
+                          double *uc, double *wc);
+
+tGrid *FindPointInGrids(tWorkUnit *wu, double x, double y, double z, double *u, double *v, double *w);
+
+double InterpolateUVW(tGrid *g, int var, double *data, double u, double v, double w);
+
+void Interpolate1d(double *uin, double *uout, double x, int n1, int n2,
+                   int n3, double *xb, double *wb);
+
+double InterpolateUW(tGrid *g, int var, double *data, double u, double w);
+
+void Interpolate1d_axisym(double *uin, double *uout, double x, int n1, int n2,
+                   double *xb, double *wb);
+
+
+void Interpolate1d_axisym_sym(double *uin, double *uout, double x, int n1, int n1_sym,
+                       int n2, double *xb, double *wb, int odd);
+
+void grid_coordtrans_xbarofx(const int sphere0, const int sphere1,
+                             double u0, double u1,
+                             double x, double y, double z,
+                             double *pu, double *pv, double *pw);
+
+// output.c ///////////////////////////////////////////////////////////////////
+
+void OutputDivider();
+
+void print2d(double *f, int n1, int n2);
+
+void print2d_nonflat(double **f, int n1, int n2);
+
+void PrintWorkGnuplot2D(tWorkUnit *wu, int res, char* fname);
+
+void OutputGrids(tWorkUnit *wu, char* dname);
+
+void OutputGrids_axisym(tWorkUnit *wu, char* dname);
+
+void OutputHorizon(tSurface *hor, char *fname);
+
+void OutputHorizon_axisym(tSurface *hor, char *fname);
+
+void OutputArray(double *a, int n, char *fname);
+
+void OutputArray3(double *a, double *b, double *c, int n, char *fname);
+
+void OutputMasses(tWorkUnit *wu, tSurface **hors, int n, char *fname);
+
+void OutputMasses_axisym(tWorkUnit *wu, tSurface **hors, int n, char *fname);
+
+// parameters.c ///////////////////////////////////////////////////////////////
+
+void ReadParfile(char *fname);
+
+// testing.c //////////////////////////////////////////////////////////////////
+
+void Kerr_set_all(double g[3][3], double dg[3][3][3], double Chr[3][3][3], double K[3][3], double x, double y, double z, double a);
+
+void Kerr_set_metric(double g[3][3], double x, double y, double z, double a);
+
+void Kerr_set_Hl(double Hl[3], double x, double y, double z, double a);
+
+void BL_set_all(double g[3][3], double dg[3][3][3], double Chr[3][3][3], double K[3][3], double x, double y, double z);
+
+void BL_set_metric(double g[3][3], double x, double y, double z);
+
+void InterpolationTest(tWorkUnit *wu, double xMax, double yMax, double zMax, int nx, int ny, int nz, int var, char *fname);
+
+void timesteppingTestFun(double *u, double *dtu);
+
+void testTimestepper();
+
+void get3KS(tWorkUnit *wu,double ***KS, double *xs, double *ys, double *zs, int nx, int ny, int nz);
+
+// util.c /////////////////////////////////////////////////////////////////////
+
+void AddHorizonToWU(tWorkUnit *wu, tSurface *new);
+
+void AddHorizon(tSurface ***horizons, int n, tSurface *new);
+
+void CopySurface(tSurface *src, tSurface **tgt);
+
+void FreeHorizon(tSurface *hor);
+
+void HorizonMinMax(tSurface *hor, double *xMin, double *yMin, double *zMin, double *xMax, double *yMax, double *zMax);
+
+void invertLU3(double a[3][3], double b[3][3], int preserve);
+
+void lubacksub(double** a, int n, int* indices, double *b);
+
+void lubacksub3(double a[3][3], int* indices, double b[]);
+
+void ludecomp(double** a, int n, int* indices, int* d);
+
+void ludecomp3(double a[3][3], int* indices, int* d);
+
+double MaxInAbs(double* a, int n);
+
+tSurface *NewSurface(int nt, int np, double x, double y, double z, double r);
+
+tSurface *loadRawSurface(char *fname);
+
+void saveRawSurface(tSurface *surf, char *fname);
+
+int isConcave(tSurface *surf, int i, int j);
+
+int isConcave_axisym(tSurface *surf, int i, int j);
+
+void OutputConcave(tSurface *surf, char *fname);
+
+void OutputConcave_axisym(tSurface *surf, char *fname);
+
+void MergeResults(tSurface **src1, tSurface **src2, tSurface ***tgt, int n1, int n2);
+
+double maxOf(double x, double y);
+
+double minOf(double x, double y);
+
+double minRad(tSurface *surf);
+
+void interpolateSurface(tSurface *src, tSurface *tgt);
+
+tSurface *newSurfaceInterpolated(tSurface *src, int nt, int np);
+
+int isInfoFile(struct dirent *entry);
+
+tGridType parseType(char *str);
+
+tGridOrientation parseOrient(char *str);
+
+void getTypeString(char* str, tGridType type, int bufSize);
+
+void getOrientString(char *str, tGridOrientation orient, int bufSize);
+
+int count_newlines(char *filename);
+
+void setup_output(char *filename);
+
+// work.c /////////////////////////////////////////////////////////////////////
+
+void ReadWork_legacy(void);
+
+void read_work();
+
+void PrintWork(void);
+
+void FreeWork(void);
+
+"""
+
+ffibuilder.cdef(ahloc3d_h)
+ffibuilder.set_source("_ahloc_3d", """#include "./../ahloc3d/src/AHloc.h" """,
+                      libraries=['AHloc','mpi'],
+                      library_dirs=['../ahloc3d/', '../ahloc3d/src', '../ahloc3d/objs'],
+                      include_dirs=['/usr/lib/x86_64-linux-gnu/openmpi/include',
+                                    '/usr/lib/x86_64-linux-gnu/openmpi/include/openmpi'])
+
+if __name__ == "__main__":
+        ffibuilder.compile(verbose=True)