From 953a8286ebbdc40cd3afb1481b1612fecde12e23 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 26 Feb 2020 18:54:40 -0500 Subject: [PATCH] Multiple changed mostly focusing on active particle search * Added optional active particle search via GRAPite. * Fixed ETICS_DTSCF in the Makefile. * Disabled ETICS_CEP by default. * Renamed and improved the initialization Python script --- Makefile | 4 +- gen-init.py | 74 ------------------------------------ init.py | 105 ++++++++++++++++++++++++++++++++++++++++++++++++++++ phi-GRAPE.c | 31 +++++++++++++++- 4 files changed, 136 insertions(+), 78 deletions(-) delete mode 100644 gen-init.py create mode 100644 init.py diff --git a/Makefile b/Makefile index 8921623..f8fe666 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CUDAHOME ?= /usr/local/cuda CPPFLAGS += -DYEBISU -DETICS OPTIMIZATION ?= 3 -DTSCF ?= 0.015625 +ETICS_DTSCF ?= 0.015625 CUDAINC = -I$(CUDAHOME)/include -I$(CUDAHOME)/samples/common/inc/ CUDALIB = -L$(CUDAHOME)/lib64 -lcudart -lcudadevrt @@ -21,7 +21,7 @@ MPICC ?= mpicc EXECUTABLE ?= phi-GRAPE.exe default: - $(MPICC) $(CPPFLAGS) $(CFLAGS) $(INC) -DDTSCF=$(DTSCF) phi-GRAPE.c -o $(EXECUTABLE) $(LIB) + $(MPICC) $(CPPFLAGS) $(CFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) phi-GRAPE.c -o $(EXECUTABLE) $(LIB) yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS)) yebisu: default diff --git a/gen-init.py b/gen-init.py deleted file mode 100644 index 899df4e..0000000 --- a/gen-init.py +++ /dev/null @@ -1,74 +0,0 @@ -import numpy as np -import argparse - -def gen_plum(N, seed=None, RMAX=10): - if not seed is None: np.random.seed(seed) - particle_list = np.zeros((N,6)) - i = 0 - while (i < N): - X1, X2, X3 = np.random.random(3) - R = 1/np.sqrt(X1**(-2/3) - 1) - if (R > RMAX): continue - Z = (1 - 2*X2)*R - X = np.sqrt(R**2 - Z**2) * np.cos(2*np.pi*X3) - Y = np.sqrt(R**2 - Z**2) * np.sin(2*np.pi*X3) - - Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25); - - X4, X5 = 0, 0 - while 0.1*X5 >= X4**2*(1-X4**2)**3.5: - X4, X5 = np.random.random(2) - - V = Ve*X4 - X6, X7 = np.random.random(2) - - Vz = (1 - 2*X6)*V; - Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7); - Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7); - - X, Y, Z = np.array([X, Y, Z])*3*np.pi/16 - Vx, Vy, Vz = np.array([Vx, Vy, Vz])/np.sqrt(3*np.pi/16) - - particle_list[i,:] = [X, Y, Z, Vx, Vy, Vz] - i += 1 - return particle_list - -parser = argparse.ArgumentParser(description='Process some integers.') -parser.add_argument('N', type=str, help='The number of particles (follow by k to multiply by 1024)') -parser.add_argument('--seed', type=int, default=42, help='Random seed') -parser.add_argument('--eps', type=np.double, default=1E-4, help='Plummer softening parameter (can be even 0)') -parser.add_argument('--t_end', type=np.double, default=4, help='end time of calculation') -parser.add_argument('--dt_disk', type=np.double, default=1, help='interval of snapshot files output (0xxx.dat)') -parser.add_argument('--dt_contr', type=np.double, default=.125, help='interval for the energy control output (contr.dat)') -parser.add_argument('--dt_bh', type=np.double, default=.125, help='interval for BH output (bh.dat & bh_nb.dat)') -parser.add_argument('--eta', type=np.double, default=.01, help='parameter for timestep determination (0.02 or 0.01)') - - -args = parser.parse_args() -try: - N = int(args.N) -except ValueError: - if args.N[-1]=='k': N = int(args.N[:-1])*1024 - else: raise ValueError('Unable to parse N.') - -f = open('phi-GRAPE.cfg', 'w') -f.write('%.4E %.4E %.4E %.4E %.4E %.4E data.con\n' % (args.eps, args.t_end, args.dt_disk, args.dt_contr, args.dt_bh, args.eta)) -f.close() - - -particle_list = gen_plum(N) -m = np.empty(N) -m[:] = 1/N - -f = open('data.con', 'w') -f.write('000000\n') -f.write('%d\n' % N) -f.write('0.0000000000E+00\n') -for i in range(N): - f.write('%06d%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e\n' % (i, m[i], particle_list[i,0], particle_list[i,1], particle_list[i,2], particle_list[i,3], particle_list[i,4], particle_list[i,5])) -f.close() - -f = open('grapite.mask', 'w') -for i in range(N): - f.write('%06d %d\n' % (i, 0)) -f.close() diff --git a/init.py b/init.py new file mode 100644 index 0000000..e5aa76c --- /dev/null +++ b/init.py @@ -0,0 +1,105 @@ +import numpy as np +import argparse + +def gen_plum(N, seed=None, RMAX=10): + if not seed is None: np.random.seed(seed) + particle_list = np.zeros((N,6)) + i = 0 + while (i < N): + X1, X2, X3 = np.random.random(3) + R = 1/np.sqrt(X1**(-2/3) - 1) + if (R > RMAX): continue + Z = (1 - 2*X2)*R + X = np.sqrt(R**2 - Z**2) * np.cos(2*np.pi*X3) + Y = np.sqrt(R**2 - Z**2) * np.sin(2*np.pi*X3) + + Ve = np.sqrt(2)*(1.0 + R**2)**(-0.25); + + X4, X5 = 0, 0 + while 0.1*X5 >= X4**2*(1-X4**2)**3.5: + X4, X5 = np.random.random(2) + + V = Ve*X4 + X6, X7 = np.random.random(2) + + Vz = (1 - 2*X6)*V; + Vx = np.sqrt(V**2 - Vz**2) * np.cos(2*np.pi*X7); + Vy = np.sqrt(V**2 - Vz**2) * np.sin(2*np.pi*X7); + + X, Y, Z = np.array([X, Y, Z])*3*np.pi/16 + Vx, Vy, Vz = np.array([Vx, Vy, Vz])/np.sqrt(3*np.pi/16) + + particle_list[i,:] = [X, Y, Z, Vx, Vy, Vz] + i += 1 + return particle_list + +def write_phi_grape_config(**kargs): + if 'file_name' in kargs: file_name = kargs['file_name'] + else: file_name = 'phi-GRAPE.cfg' + f = open(file_name, 'w') + f.write('%.4E %.4E %.4E %.4E %.4E %.4E data.con\n' % (kargs['eps'], kargs['t_end'], kargs['dt_disk'], kargs['dt_contr'], kargs['dt_bh'], kargs['eta'])) + f.close() + +def write_particles(particle_list, m=None, file_name='data.con'): + N = particle_list.shape[0] + if m is None: + m = np.empty(N) + m[:] = 1/N + f = open(file_name, 'w') + f.write('000000\n') + f.write('%d\n' % N) + f.write('0.0000000000E+00\n') + for i in range(N): + f.write('%06d%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e\n' % (i, m[i], particle_list[i,0], particle_list[i,1], particle_list[i,2], particle_list[i,3], particle_list[i,4], particle_list[i,5])) + f.close() + +def gen_mask(particle_list, frac): + N = particle_list.shape[0] + if frac==0: + mask = np.ones(N, dtype=int) + elif frac==1: + mask = np.zeros(N, dtype=int) + else: + X = particle_list[:,:3] + V = particle_list[:,3:] + L = np.linalg.norm(np.cross(X, V), axis=1) + L_sorted = L.copy() + L_sorted.sort() + L_cutoff = L_sorted[int(N*frac)] + mask = np.array(L > L_cutoff, dtype=int) + return mask + +def write_mask(mask, file_name='grapite.mask'): + f = open(file_name, 'w') + for i in range(len(mask)): + f.write('%06d %d\n' % (i, mask[i])) + f.close() + +if __name__=='__main__': + parser = argparse.ArgumentParser(description='Process some integers.') + parser.add_argument('N', type=str, help='number of particles (follow by k to multiply by 1024)') + parser.add_argument('--seed', type=int, default=42, help='random seed') + parser.add_argument('--eps', type=np.double, default=1E-4, help='Plummer softening parameter (can be even 0)') + parser.add_argument('--t_end', type=np.double, default=4, help='end time of calculation') + parser.add_argument('--dt_disk', type=np.double, default=1, help='interval of snapshot files output (0xxx.dat)') + parser.add_argument('--dt_contr', type=np.double, default=.125, help='interval for the energy control output (contr.dat)') + parser.add_argument('--dt_bh', type=np.double, default=.125, help='interval for BH output (bh.dat & bh_nb.dat)') + parser.add_argument('--eta', type=np.double, default=.01, help='parameter for timestep determination (0.02 or 0.01)') + parser.add_argument('--frac', type=np.double, default=0, help='fraction of collisional particles (by angular momentum)') + args = parser.parse_args() + + try: + N = int(args.N) + except ValueError: + if args.N[-1]=='k': N = int(args.N[:-1])*1024 + else: raise ValueError('Unable to parse N.') + + write_phi_grape_config(**vars(args)) + + particle_list = gen_plum(N, seed=args.seed) + + write_particles(particle_list) + + mask = gen_mask(particle_list, args.frac) + + write_mask(mask) diff --git a/phi-GRAPE.c b/phi-GRAPE.c index f4e6981..df15d03 100644 --- a/phi-GRAPE.c +++ b/phi-GRAPE.c @@ -134,8 +134,10 @@ #ifdef ETICS #include "grapite.h" // why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations. -#define ETICS_CEP -#define ETICS_DTSCF 0.125 +//#define ETICS_CEP +#ifndef ETICS_DTSCF +#error "ETICS_DTSCF must be defined" +#endif #endif #define TIMING @@ -163,6 +165,10 @@ // #define ACT_DEF_LL +#if defined(ACT_DEF_LL) && defined(ACT_DEF_GRAPITE) +#error "Contradicting preprocessor flags!" +#endif + /****************************************************************************/ #include #include @@ -5630,6 +5636,9 @@ for(i=0; i 0) + for(int i=0; i