diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml new file mode 100644 index 0000000..3855b91 --- /dev/null +++ b/.github/workflows/python-package.yml @@ -0,0 +1,11 @@ +name: black-action +on: [push, pull_request] +jobs: + linter_name: + name: runner / black formatter + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: rickstaa/action-black@v1 + with: + black_args: "-l 119 ." \ No newline at end of file diff --git a/MakeCloud.py b/MakeCloud.py index 92f1045..241b67f 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -116,7 +116,11 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): if i == j: vk_new[i] += sol_weight * VK[j] vk_new[i] += ( - (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j] + (1 - 2 * sol_weight) + * freq3d[i] + * freq3d[j] + / (kSqr + 1e-300) + * VK[j] ) vk_new[:, kSqr == 0] = 0.0 VK = vk_new @@ -124,8 +128,12 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): vel = np.array( [fftpack.ifftn(vk).real for vk in VK] ) # transform back to real space - vel -= np.average(vel, axis=(1, 2, 3))[:, np.newaxis, np.newaxis, np.newaxis] - vel = vel / np.sqrt(np.sum(vel**2, axis=0).mean()) # normalize so that RMS is 1 + vel -= np.average(vel, axis=(1, 2, 3))[ + :, np.newaxis, np.newaxis, np.newaxis + ] + vel = vel / np.sqrt( + np.sum(vel**2, axis=0).mean() + ) # normalize so that RMS is 1 return np.array(vel) @@ -179,8 +187,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): print("Downloading glass file...") urllib.request.urlretrieve( - "http://www.tapir.caltech.edu/~mgrudich/glass_orig.npy", glass_path -# "https://data.obs.carnegiescience.edu/starforge/glass_orig.npy", glass_path + "http://www.tapir.caltech.edu/~mgrudich/glass_orig.npy", + glass_path, + # "https://data.obs.carnegiescience.edu/starforge/glass_orig.npy", glass_path ) if localdir: @@ -220,7 +229,8 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): ) + ("_%d" % seed) + ( - "_collision_%g_%g_%g_%s" % (impact_dist, impact_param, v_impact, impact_axis) + "_collision_%g_%g_%g_%s" + % (impact_dist, impact_param, v_impact, impact_axis) if impact_dist > 0 else "" ) @@ -233,7 +243,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): delta_m_solar = delta_m / mass_unit rho_avg = 3 * M_gas / R**3 / (4 * np.pi) if delta_m_solar < 0.1: # if we're doing something marginally IMF-resolving - softening = 3.11e-5 # ~6.5 AU, minimum sink radius is 2.8 times that (~18 AU) + softening = ( + 3.11e-5 # ~6.5 AU, minimum sink radius is 2.8 times that (~18 AU) + ) ncrit = 1e13 # ~100x the opacity limit else: # something more FIRE-like, where we rely on a sub-grid prescription turning gas into star particles softening = 0.1 @@ -256,7 +268,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): ) # ST_Energy sets the dissipation rate of SPECIFIC energy ~ v^2 / (L/v) ~ v^3/L paramsfile = str( - open(os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r").read() + open( + os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r" + ).read() ) jet_particle_mass = min(delta_m, max(1e-4, delta_m / 10.0)) @@ -295,7 +309,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): } for k, r in replacements.items(): - paramsfile = paramsfile.replace(k,(r if isinstance(r,str) else '{:.2e}'.format(r))) + paramsfile = paramsfile.replace( + k, (r if isinstance(r, str) else "{:.2e}".format(r)) + ) open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write(paramsfile) if makebox: @@ -306,12 +322,15 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): replacements_box["TURB_MAXLAMBDA"] = int(100 * L * 2) / 100 paramsfile = str( open( - os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r" + os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), + "r", ).read() ) for k in replacements_box.keys(): paramsfile = paramsfile.replace(k, str(replacements_box[k])) - open("params_" + filename.replace(".hdf5", "") + "_BOX.txt", "w").write(paramsfile) + open("params_" + filename.replace(".hdf5", "") + "_BOX.txt", "w").write( + paramsfile + ) if makecylinder: # Get cylinder params R_cyl = R * np.sqrt( @@ -346,12 +365,15 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): replacements_cyl["TURB_KMAX"] = int(100 * 4 * np.pi / (R_cyl) + 1) / 100.0 paramsfile = str( open( - os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r" + os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), + "r", ).read() ) for k in replacements_cyl.keys(): paramsfile = paramsfile.replace(k, str(replacements_cyl[k])) - open("params_" + filename.replace(".hdf5", "") + "_CYL.txt", "w").write(paramsfile) + open("params_" + filename.replace(".hdf5", "") + "_CYL.txt", "w").write( + paramsfile + ) if param_only: print("Parameters only run, exiting...") @@ -414,9 +436,16 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): B = np.c_[np.zeros(N_gas), np.zeros(N_gas), np.ones(N_gas)] vA_unit = ( - 3.429e8 * B_unit * (M_gas) ** -0.5 * R**1.5 * np.sqrt(4 * np.pi / 3) / v_unit + 3.429e8 + * B_unit + * (M_gas) ** -0.5 + * R**1.5 + * np.sqrt(4 * np.pi / 3) + / v_unit ) # alfven speed for unit magnetic field -uB = 0.5 * M_gas * vA_unit**2 # magnetic energy we would have for unit magnetic field +uB = ( + 0.5 * M_gas * vA_unit**2 +) # magnetic energy we would have for unit magnetic field if bfixed > 0: B = B * bfixed else: @@ -432,7 +461,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): phi += phimode * np.sin(2 * phi) / 2 x = ( r[:, np.newaxis] - * np.c_[np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)] + * np.c_[ + np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta) + ] ) if makecylinder: @@ -519,13 +550,17 @@ def ind_in_cylinder(x, L_cyl, R_cyl): x_warm = x_warm[np.max(np.abs(x_warm), axis=1) < boxsize / 2] N_warm = len(x_warm) R_warm = (x_warm * x_warm).sum(1) ** 0.5 - mgas = np.concatenate([mgas, np.clip(dm * (R_warm / (3 * R)) ** 3, dm, np.inf)]) + mgas = np.concatenate( + [mgas, np.clip(dm * (R_warm / (3 * R)) ** 3, dm, np.inf)] + ) else: x_warm = boxsize * np.random.rand(N_warm, 3) - boxsize / 2 if impact_dist == 0: x_warm = x_warm[np.sum(x_warm**2, axis=1) > R**2] N_warm = len(x_warm) - mgas = np.concatenate([mgas, np.repeat(mgas.sum() / len(mgas), N_warm)]) + mgas = np.concatenate( + [mgas, np.repeat(mgas.sum() / len(mgas), N_warm)] + ) x = np.concatenate([x, x_warm]) v = np.concatenate([v, np.zeros((N_warm, 3))]) Bmag = np.average(np.sum(B**2, axis=1)) ** 0.5 @@ -578,7 +613,14 @@ def ind_in_cylinder(x, L_cyl, R_cyl): 0, (1 if M_star > 0 else 0), ] -F["Header"].attrs["NumPart_Total"] = [len(mgas), 0, 0, 0, 0, (1 if M_star > 0 else 0)] +F["Header"].attrs["NumPart_Total"] = [ + len(mgas), + 0, + 0, + 0, + 0, + (1 if M_star > 0 else 0), +] F["Header"].attrs["BoxSize"] = boxsize F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas) @@ -591,7 +633,9 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F.create_group("PartType5") # Let's add the sink at the center F["PartType5"].create_dataset("Masses", data=np.array([M_star])) - F["PartType5"].create_dataset("Coordinates", data=[x_star]) # at the center + F["PartType5"].create_dataset( + "Coordinates", data=[x_star] + ) # at the center F["PartType5"].create_dataset("Velocities", data=[v_star]) # at rest F["PartType5"].create_dataset( "ParticleIDs", data=np.array([F["PartType0/ParticleIDs"][:].max() + 1]) @@ -633,7 +677,9 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F["PartType5"].create_dataset( "ProtoStellarRadius_inSolar", data=np.array([R_ZAMS]) ) # Sinkradius set to softening - F["PartType5"].create_dataset("StarLuminosity_Solar", data=np.array([0.0])) # dummy + F["PartType5"].create_dataset( + "StarLuminosity_Solar", data=np.array([0.0]) + ) # dummy F["PartType5"].create_dataset("Mass_D", data=np.array([0.0])) # No D left if magnetic_field > 0.0: @@ -651,7 +697,8 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas[: len(mgas)]) F["PartType0"].create_dataset( - "Coordinates", data=np.random.rand(len(mgas), 3) * F["Header"].attrs["BoxSize"] + "Coordinates", + data=np.random.rand(len(mgas), 3) * F["Header"].attrs["BoxSize"], ) F["PartType0"].create_dataset("Velocities", data=np.zeros((len(mgas), 3))) F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(len(mgas))) @@ -672,7 +719,9 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F["PartType0"].create_dataset("Masses", data=mgas) F["PartType0"].create_dataset("Coordinates", data=x_cyl) F["PartType0"].create_dataset("Velocities", data=v_cyl) - F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(N_gas + N_warm_cyl)) + F["PartType0"].create_dataset( + "ParticleIDs", data=1 + np.arange(N_gas + N_warm_cyl) + ) F["PartType0"].create_dataset("InternalEnergy", data=u) if magnetic_field > 0.0: F["PartType0"].create_dataset("MagneticField", data=B_cyl) diff --git a/MakeCloud_Bondi.py b/MakeCloud_Bondi.py index 1d97a37..aaeb10f 100644 --- a/MakeCloud_Bondi.py +++ b/MakeCloud_Bondi.py @@ -33,7 +33,8 @@ import h5py import os from docopt import docopt -#from pykdgrav import Potential + +# from pykdgrav import Potential arguments = docopt(__doc__) GMC_units = arguments["--GMC_units"] length_unit = float(arguments["--length_unit"]) @@ -43,165 +44,267 @@ length_unit = 1.0 mass_unit = 1.0 v_unit = 1.0 - + R_over_Rsonic = float(arguments["--R_over_Rsonic"]) Rs_over_Rsonic = float(arguments["--Rs_over_Rsonic"]) -rho_gas = float(arguments["--rho_gas"])/(mass_unit/(length_unit**3.0)) #msolar/pc^3 to code units -N_gas = int(float(arguments["--N"])+0.5) -M_BH = float(arguments["--MBH"])/mass_unit #mstar to code units -R_cut_out = float(arguments["--R_cut_out"])/length_unit +rho_gas = float(arguments["--rho_gas"]) / ( + mass_unit / (length_unit**3.0) +) # msolar/pc^3 to code units +N_gas = int(float(arguments["--N"]) + 0.5) +M_BH = float(arguments["--MBH"]) / mass_unit # mstar to code units +R_cut_out = float(arguments["--R_cut_out"]) / length_unit filename = arguments["--filename"] localdir = arguments["--localdir"] +# set gravitational constant +G = 4325.69 / length_unit / (v_unit**2) * (mass_unit) -#set gravitational constant -G = 4325.69 /length_unit / (v_unit**2) * (mass_unit) - if localdir: turb_path = "turb" glass_path = "glass_256.npy" -#get rsink radius from sonic radius -csound=200/v_unit #200 m/s -rsonic_Bondi = G*M_BH/2.0/(csound**2)/length_unit -msonic=4.0*np.pi/3.0*rho_gas*(rsonic_Bondi**3.0) -Rsink = Rs_over_Rsonic*rsonic_Bondi -R = R_over_Rsonic*rsonic_Bondi -print( "Radius set as %g pc"%(R*length_unit)) +# get rsink radius from sonic radius +csound = 200 / v_unit # 200 m/s +rsonic_Bondi = G * M_BH / 2.0 / (csound**2) / length_unit +msonic = 4.0 * np.pi / 3.0 * rho_gas * (rsonic_Bondi**3.0) +Rsink = Rs_over_Rsonic * rsonic_Bondi +R = R_over_Rsonic * rsonic_Bondi +print("Radius set as %g pc" % (R * length_unit)) if arguments["--boxsize"] is not None: -# print((arguments["--boxsize"])) - boxsize = float(arguments["--boxsize"])*length_unit + # print((arguments["--boxsize"])) + boxsize = float(arguments["--boxsize"]) * length_unit else: - boxsize = 10*R -center = np.ones(3)*boxsize/2.0 -res_effective = int(N_gas**(1.0/3.0)+0.5) - -#Load Bondi Spherical solution -IC_data=np.loadtxt("BONDI_SOLUTION.dat") #load IC data from Hubber -data_r=IC_data[:,0]*rsonic_Bondi #coordinate originally in Rsonic -data_within_R_index=(data_r<=R);data_r=data_r[data_within_R_index] #restrict to within R radius -data_mass_in_r=IC_data[data_within_R_index,1]*msonic#gas mass in units of 4*pi*rho0*Rsonic^3/3 -data_density=IC_data[data_within_R_index,2]*rho_gas #density originally in rho0 -data_vr=IC_data[data_within_R_index,3]*csound #velocity originally in csound -M_gas=np.max(data_mass_in_r) #normalize to total mass -#Get glass -mgas = np.repeat(M_gas/N_gas, N_gas) #gas mass init -x = 2*(np.load(glass_path)-0.5) #load glass to have basic structure -Nx=len(x[:,0]); #original glass size -r = np.sum(x**2,axis=1)**0.5 #calculate radius -x = x[r.argsort()][:N_gas] #sort the particles by radius, it should be spherically symmetric, now take the right number of particles -x *= (float(Nx) / N_gas * 4*np.pi/3 / 8)**(1./3) #scale up the sphere to have unit radius -r = np.sum(x**2,axis=1)**0.5 #recalculate radius -#Strech the particles based on how much mass is in a given radius -int_mass=np.cumsum(mgas) #integrated mass -newr=np.interp(int_mass, data_mass_in_r, data_r) #calculate the radius the particles should be at to have the sane mass within the radius a sthe solution -stretch=newr/r #stretch factor -x[:,0] *= stretch;x[:,1] *= stretch;x[:,2] *= stretch; -u = np.ones_like(mgas)*0.101*((1000/v_unit)**2)/2.0 #/2 needed because it is molecular -rho = np.interp(int_mass, data_mass_in_r, data_density) #interpolate the density -h = (32*mgas/rho)**(1./3) -vr = np.interp(int_mass, data_mass_in_r, data_vr) #interpolate the velocity -v = -x; v[:,0] *= vr/newr; v[:,1] *= vr/newr; v[:,2] *= vr/newr; -#print chek data -print( 'mdot avg: %g std %g'%(np.mean(data_density*data_r*data_r*data_vr*np.pi*4.0),np.std(data_density*data_r*data_r*data_vr*np.pi*4.0))) -print( 'density') -print( data_density) -print( 'density_alt') -altdens=np.diff(data_mass_in_r)/np.diff(data_r)/(4.0*np.pi*data_r[1:]*data_r[1:]) -print( altdens) -print( 'relative') -print( data_density[1:]/altdens) -#Calculate NEWSNK parameters -print( 'Calculating t_radial...') -for z in [0.125,0.5,1.0,2.0,4.0]: - print( '\t At %g Rsonic = %g pc'%(z,(rsonic_Bondi*z*length_unit))) - ind=(newr<=(z*rsonic_Bondi)) - print( '\t mass enclosed %g in msun'%(np.sum(mgas[ind])*mass_unit)) - #tot_weight=np.sum(mgas[ind]/rho[ind]) - tot_weight=4.0/3.0*np.pi*np.max(newr[ind])**3 - x_dot_v=np.sum(x*v,axis=1) - t_rad=-np.sum(mgas[ind])*tot_weight/np.sum(4.0*np.pi*(x_dot_v*newr*mgas)[ind]) - print( '\t t_rad=%g in code units'%(t_rad)) - print( '\t m/t_rad=%g in code units'%(np.sum(mgas[ind])/t_rad)) - #from IC - ind=np.arange(len(data_r))[(data_r<=(z*rsonic_Bondi))] - ind2=np.argmax(data_r[ind]) - dm_ic=np.diff(data_mass_in_r) - t_rad_IC=(data_mass_in_r[ind2]*(4.0/3.0*np.pi)*data_r[ind2]**3)/(4.0*np.pi*np.sum(data_r[ind]*data_r[ind]*data_vr[ind]*dm_ic[ind])) - print( '\t t_rad_IC=%g in code units'%(t_rad_IC)) - print( '\t m_IC/t_rad_IC=%g in code units'%(data_mass_in_r[ind2]/t_rad_IC)) -#Calculate flux -print( 'Calculating density ...') -dr=np.diff(newr) -rho_alt=(mgas[1:]/dr)/(4.0*np.pi*(newr[1:]**2)) -print( rho) -print( rho_alt) -print( rho[1:]/rho_alt) -#keep the one sthat are not too close -ind=newr>R_cut_out -print("Removing %d particles that are inside R_cut_out of %g"%((N_gas-np.sum(ind)),R_cut_out)) -N_gas=np.sum(ind) -M_gas=np.sum(mgas[ind]) - -NGBvals=[1,2,5,10,20,50,100,200,400,800,1600] + boxsize = 10 * R +center = np.ones(3) * boxsize / 2.0 +res_effective = int(N_gas ** (1.0 / 3.0) + 0.5) + +# Load Bondi Spherical solution +IC_data = np.loadtxt("BONDI_SOLUTION.dat") # load IC data from Hubber +data_r = IC_data[:, 0] * rsonic_Bondi # coordinate originally in Rsonic +data_within_R_index = data_r <= R +data_r = data_r[data_within_R_index] # restrict to within R radius +data_mass_in_r = ( + IC_data[data_within_R_index, 1] * msonic +) # gas mass in units of 4*pi*rho0*Rsonic^3/3 +data_density = ( + IC_data[data_within_R_index, 2] * rho_gas +) # density originally in rho0 +data_vr = ( + IC_data[data_within_R_index, 3] * csound +) # velocity originally in csound +M_gas = np.max(data_mass_in_r) # normalize to total mass +# Get glass +mgas = np.repeat(M_gas / N_gas, N_gas) # gas mass init +x = 2 * (np.load(glass_path) - 0.5) # load glass to have basic structure +Nx = len(x[:, 0]) +# original glass size +r = np.sum(x**2, axis=1) ** 0.5 # calculate radius +x = x[r.argsort()][ + :N_gas +] # sort the particles by radius, it should be spherically symmetric, now take the right number of particles +x *= (float(Nx) / N_gas * 4 * np.pi / 3 / 8) ** ( + 1.0 / 3 +) # scale up the sphere to have unit radius +r = np.sum(x**2, axis=1) ** 0.5 # recalculate radius +# Strech the particles based on how much mass is in a given radius +int_mass = np.cumsum(mgas) # integrated mass +newr = np.interp( + int_mass, data_mass_in_r, data_r +) # calculate the radius the particles should be at to have the sane mass within the radius a sthe solution +stretch = newr / r # stretch factor +x[:, 0] *= stretch +x[:, 1] *= stretch +x[:, 2] *= stretch +u = ( + np.ones_like(mgas) * 0.101 * ((1000 / v_unit) ** 2) / 2.0 +) # /2 needed because it is molecular +rho = np.interp( + int_mass, data_mass_in_r, data_density +) # interpolate the density +h = (32 * mgas / rho) ** (1.0 / 3) +vr = np.interp(int_mass, data_mass_in_r, data_vr) # interpolate the velocity +v = -x +v[:, 0] *= vr / newr +v[:, 1] *= vr / newr +v[:, 2] *= vr / newr +# print chek data +print( + "mdot avg: %g std %g" + % ( + np.mean(data_density * data_r * data_r * data_vr * np.pi * 4.0), + np.std(data_density * data_r * data_r * data_vr * np.pi * 4.0), + ) +) +print("density") +print(data_density) +print("density_alt") +altdens = ( + np.diff(data_mass_in_r) + / np.diff(data_r) + / (4.0 * np.pi * data_r[1:] * data_r[1:]) +) +print(altdens) +print("relative") +print(data_density[1:] / altdens) +# Calculate NEWSNK parameters +print("Calculating t_radial...") +for z in [0.125, 0.5, 1.0, 2.0, 4.0]: + print("\t At %g Rsonic = %g pc" % (z, (rsonic_Bondi * z * length_unit))) + ind = newr <= (z * rsonic_Bondi) + print("\t mass enclosed %g in msun" % (np.sum(mgas[ind]) * mass_unit)) + # tot_weight=np.sum(mgas[ind]/rho[ind]) + tot_weight = 4.0 / 3.0 * np.pi * np.max(newr[ind]) ** 3 + x_dot_v = np.sum(x * v, axis=1) + t_rad = ( + -np.sum(mgas[ind]) + * tot_weight + / np.sum(4.0 * np.pi * (x_dot_v * newr * mgas)[ind]) + ) + print("\t t_rad=%g in code units" % (t_rad)) + print("\t m/t_rad=%g in code units" % (np.sum(mgas[ind]) / t_rad)) + # from IC + ind = np.arange(len(data_r))[(data_r <= (z * rsonic_Bondi))] + ind2 = np.argmax(data_r[ind]) + dm_ic = np.diff(data_mass_in_r) + t_rad_IC = ( + data_mass_in_r[ind2] * (4.0 / 3.0 * np.pi) * data_r[ind2] ** 3 + ) / ( + 4.0 + * np.pi + * np.sum(data_r[ind] * data_r[ind] * data_vr[ind] * dm_ic[ind]) + ) + print("\t t_rad_IC=%g in code units" % (t_rad_IC)) + print( + "\t m_IC/t_rad_IC=%g in code units" % (data_mass_in_r[ind2] / t_rad_IC) + ) +# Calculate flux +print("Calculating density ...") +dr = np.diff(newr) +rho_alt = (mgas[1:] / dr) / (4.0 * np.pi * (newr[1:] ** 2)) +print(rho) +print(rho_alt) +print(rho[1:] / rho_alt) +# keep the one sthat are not too close +ind = newr > R_cut_out +print( + "Removing %d particles that are inside R_cut_out of %g" + % ((N_gas - np.sum(ind)), R_cut_out) +) +N_gas = np.sum(ind) +M_gas = np.sum(mgas[ind]) + +NGBvals = [1, 2, 5, 10, 20, 50, 100, 200, 400, 800, 1600] for ngb in NGBvals: - print("Radius at Ngbfactor of %g is %g pc"%(ngb,newr[32*ngb])) + print("Radius at Ngbfactor of %g is %g pc" % (ngb, newr[32 * ngb])) -#center coordinates -x = x + boxsize/2.0 +# center coordinates +x = x + boxsize / 2.0 print("Writing snapshot...") if filename is None: - filename = "rho%3.2g_"%(rho_gas*(mass_unit/(length_unit**3))) + ("MBH%g_"%(mass_unit*M_BH) if M_BH>0 else "") + "R_over_Rsonic%g_Res%d_Rs_over_Rsonic%g"%(R_over_Rsonic,res_effective,Rs_over_Rsonic) + ".hdf5" - filename = filename.replace("+","").replace('e0','e') + filename = ( + "rho%3.2g_" % (rho_gas * (mass_unit / (length_unit**3))) + + ("MBH%g_" % (mass_unit * M_BH) if M_BH > 0 else "") + + "R_over_Rsonic%g_Res%d_Rs_over_Rsonic%g" + % (R_over_Rsonic, res_effective, Rs_over_Rsonic) + + ".hdf5" + ) + filename = filename.replace("+", "").replace("e0", "e") filename = "".join(filename.split()) - -F=h5py.File(filename, 'w') + +F = h5py.File(filename, "w") F.create_group("PartType0") F.create_group("Header") -F["Header"].attrs["NumPart_ThisFile"] = [N_gas,0,0,0,0,(1 if M_BH>0 else 0)] -F["Header"].attrs["NumPart_Total"] = [N_gas,0,0,0,0,(1 if M_BH>0 else 0)] -F["Header"].attrs["MassTable"] = [M_gas/N_gas,0,0,0,0, M_BH] +F["Header"].attrs["NumPart_ThisFile"] = [ + N_gas, + 0, + 0, + 0, + 0, + (1 if M_BH > 0 else 0), +] +F["Header"].attrs["NumPart_Total"] = [ + N_gas, + 0, + 0, + 0, + 0, + (1 if M_BH > 0 else 0), +] +F["Header"].attrs["MassTable"] = [M_gas / N_gas, 0, 0, 0, 0, M_BH] F["Header"].attrs["BoxSize"] = boxsize F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas[ind]) -F["PartType0"].create_dataset("Coordinates", data=x[ind,:]) -F["PartType0"].create_dataset("Velocities", data=v[ind,:]) -F["PartType0"].create_dataset("ParticleIDs", data=np.arange(N_gas)+1) +F["PartType0"].create_dataset("Coordinates", data=x[ind, :]) +F["PartType0"].create_dataset("Velocities", data=v[ind, :]) +F["PartType0"].create_dataset("ParticleIDs", data=np.arange(N_gas) + 1) F["PartType0"].create_dataset("InternalEnergy", data=u[ind]) F["PartType0"].create_dataset("Density", data=rho[ind]) F["PartType0"].create_dataset("SmoothingLength", data=h[ind]) -if M_BH>0: +if M_BH > 0: F.create_group("PartType5") F["PartType5"].create_dataset("Masses", data=np.array([M_BH])) F["PartType5"].create_dataset("BH_Mass", data=np.array([M_BH])) F["PartType5"].create_dataset("Coordinates", data=center) - F["PartType5"].create_dataset("Velocities", data=v[0,:]*0.0) - F["PartType5"].create_dataset("ParticleIDs", data=np.array([N_gas+1])) + F["PartType5"].create_dataset("Velocities", data=v[0, :] * 0.0) + F["PartType5"].create_dataset("ParticleIDs", data=np.array([N_gas + 1])) F["PartType5"].create_dataset("SinkRadius", data=np.array([Rsink])) F.close() -if GMC_units: - delta_m = M_gas/N_gas - rhocrit = 421/ delta_m**2 - rho_avg = 3*M_gas/(R**3)/(4*np.pi) - softening = 0.000173148 # 100AU/2.8 #(delta_m/rhocrit)**(1./3) - ncrit = 1.0e11 #8920 / delta_m**2 +if GMC_units: + delta_m = M_gas / N_gas + rhocrit = 421 / delta_m**2 + rho_avg = 3 * M_gas / (R**3) / (4 * np.pi) + softening = 0.000173148 # 100AU/2.8 #(delta_m/rhocrit)**(1./3) + ncrit = 1.0e11 # 8920 / delta_m**2 tff = 8.275e-3 * rho_avg**-0.5 - mdot_Bondi = np.exp(1.5)*np.pi*(G**2)*(M_BH**2)*rho_gas/(csound**3) - tend_Bondi = 2.0*(2.0*G*M_BH/(csound**3)) - print( 'Bondi accretion parameters: \n \t Mgas:\t\t', M_gas, '\n \t Mdot:\t\t',mdot_Bondi, '\n \t Rsonic:\t',rsonic_Bondi,'\n \t Rsink:\t\t',Rsink, '\n \t t_end:\t\t',tend_Bondi, '\n \t m_acc:\t\t',tend_Bondi*mdot_Bondi, '\n \t f_acc:\t\t',tend_Bondi*mdot_Bondi) - paramsfile = str(open(os.path.realpath(__file__).replace("MakeCloud_Bondi.py","params.txt"), 'r').read()) + mdot_Bondi = ( + np.exp(1.5) * np.pi * (G**2) * (M_BH**2) * rho_gas / (csound**3) + ) + tend_Bondi = 2.0 * (2.0 * G * M_BH / (csound**3)) + print( + "Bondi accretion parameters: \n \t Mgas:\t\t", + M_gas, + "\n \t Mdot:\t\t", + mdot_Bondi, + "\n \t Rsonic:\t", + rsonic_Bondi, + "\n \t Rsink:\t\t", + Rsink, + "\n \t t_end:\t\t", + tend_Bondi, + "\n \t m_acc:\t\t", + tend_Bondi * mdot_Bondi, + "\n \t f_acc:\t\t", + tend_Bondi * mdot_Bondi, + ) + paramsfile = str( + open( + os.path.realpath(__file__).replace( + "MakeCloud_Bondi.py", "params.txt" + ), + "r", + ).read() + ) - replacements = {"NAME": "../ICs/"+filename.replace(".hdf5",""), "DTSNAP": tend_Bondi/200, "SOFTENING": softening, "GASSOFT": 2.0e-8, "TMAX": tend_Bondi, "RHOMAX": ncrit, "BOXSIZE": boxsize, "OUTFOLDER": "output_%g"%(Rs_over_Rsonic)} + replacements = { + "NAME": "../ICs/" + filename.replace(".hdf5", ""), + "DTSNAP": tend_Bondi / 200, + "SOFTENING": softening, + "GASSOFT": 2.0e-8, + "TMAX": tend_Bondi, + "RHOMAX": ncrit, + "BOXSIZE": boxsize, + "OUTFOLDER": "output_%g" % (Rs_over_Rsonic), + } print(replacements["NAME"]) -# print(paramsfile) + # print(paramsfile) for k in replacements.keys(): - paramsfile = paramsfile.replace(k, str(replacements[k])) - open("params_"+filename.replace(".hdf5","")+".txt", "w").write(paramsfile) - - + paramsfile = paramsfile.replace(k, str(replacements[k])) + open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write( + paramsfile + ) diff --git a/calcGMCProps.py b/calcGMCProps.py index 2bc6ba8..9e5a012 100644 --- a/calcGMCProps.py +++ b/calcGMCProps.py @@ -1,105 +1,111 @@ -''' +""" Simple routines to calculate the cloud's: radius, average density, surface density, free-fall time (in years), and number of particles (N) given a cloud mass and cell mass resolution to choose Written by: Anna Rosen, 11/2/22 -''' - +""" import math import numpy -MSUN = 1.989e33 #g -PCCM = 3.086e18 #cm -G = 6.6743e-8 #cm^3 g^-1 s -SECYR = 3600*24*365.25 #s - -def calcRhoAve(Mcl, Rcl = None, Sigma = 1): - ''' - Inputs: - #cloud mass = Mcl [Msun] - #cloud radius = Rcl [pc] - #Sigma = Cloud surface density, default = 1 g/cm^2 - - Output: - rho_ave [g/cm^3] - ''' - if Rcl is None and Sigma > 0.0: - print('No Rcl supplied, calc. Rcl assuming Sigma = %.2e g/cm^2' % Sigma) - Rcl = calcRcl(Mcl, Sigma) - else: - print('Error: must supple Rcl (in pc) or Sigma > 0') - M = Mcl*MSUN - R = Rcl*PCCM - rho_ave = 3*M/(4*math.pi * R**3.) - print('Average cloud density = %.2e g/cm^3' % rho_ave) - return(rho_ave) +MSUN = 1.989e33 # g +PCCM = 3.086e18 # cm +G = 6.6743e-8 # cm^3 g^-1 s +SECYR = 3600 * 24 * 365.25 # s + + +def calcRhoAve(Mcl, Rcl=None, Sigma=1): + """ + Inputs: + #cloud mass = Mcl [Msun] + #cloud radius = Rcl [pc] + #Sigma = Cloud surface density, default = 1 g/cm^2 + + Output: + rho_ave [g/cm^3] + """ + if Rcl is None and Sigma > 0.0: + print( + "No Rcl supplied, calc. Rcl assuming Sigma = %.2e g/cm^2" % Sigma + ) + Rcl = calcRcl(Mcl, Sigma) + else: + print("Error: must supple Rcl (in pc) or Sigma > 0") + M = Mcl * MSUN + R = Rcl * PCCM + rho_ave = 3 * M / (4 * math.pi * R**3.0) + print("Average cloud density = %.2e g/cm^3" % rho_ave) + return rho_ave def calcRcl(Mcl, Sigma): - ''' - Inputs: - #cloud mass = Mcl [Msun] - #cloud surface density = Sigma [g/cm^2] - - Output: - Rcl [pc] - ''' - M = Mcl * MSUN - Rcl = (M/(math.pi *Sigma))**0.5/PCCM - print('Cloud Radius = %.2f pc' % Rcl) - return(Rcl) + """ + Inputs: + #cloud mass = Mcl [Msun] + #cloud surface density = Sigma [g/cm^2] -def calc_SurfaceDensity(Mcl, Rcl): - ''' - Inputs: - #cloud mass = Mcl [Msun] - #cloud radius = Sigma [g/cm^2] - - Output: - Sigma [g/cm^2] - ''' - M = Mcl*MSUN - R = Rcl*PCCM - Sigma = M/(math.pi * R**3.) - print('Cloud surface density = %.2e g/cm^2' % Sigma) - return(rho_ave) - -def calc_tff(Mcl, Rcl = None, Sigma = 1): - ''' - Inputs: - #cloud mass = Mcl [Msun] - #cloud surface density = Sigma [g/cm^2] - - Output: - tff [yr] - ''' - if Rcl is None and Sigma > 0.0: - print('No Rcl supplied, calc. Rcl assuming Sigma = %.2e g/cm^2' % Sigma) - Rcl = calcRcl(Mcl, Sigma) - else: - print('Error: must supple Rcl (in pc) or Sigma > 0') - - rho_ave = calcRhoAve(Mcl, Rcl=Rcl, Sigma = 1) - - tff = (3*math.pi/(32*G*rho_ave))**0.5/SECYR - print('tff = %.2f yr' %tff) - return(tff) - -def calc_Npart(Mcl, massRes = 1e-3): - ''' - Inputs: - #cloud mass = Mcl [Msun] - #mass/cell [Msun] - - Output: - N [particles] - ''' - - N = Mcl/massRes - print('For Mcl = %.2e Msun, Mcell = %.2e' % (Mcl, massRes)) - print('N = %.2e particles' % N) - return N + Output: + Rcl [pc] + """ + M = Mcl * MSUN + Rcl = (M / (math.pi * Sigma)) ** 0.5 / PCCM + print("Cloud Radius = %.2f pc" % Rcl) + return Rcl + +def calc_SurfaceDensity(Mcl, Rcl): + """ + Inputs: + #cloud mass = Mcl [Msun] + #cloud radius = Sigma [g/cm^2] + + Output: + Sigma [g/cm^2] + """ + M = Mcl * MSUN + R = Rcl * PCCM + Sigma = M / (math.pi * R**3.0) + print("Cloud surface density = %.2e g/cm^2" % Sigma) + return rho_ave + + +def calc_tff(Mcl, Rcl=None, Sigma=1): + """ + Inputs: + #cloud mass = Mcl [Msun] + #cloud surface density = Sigma [g/cm^2] + + Output: + tff [yr] + """ + if Rcl is None and Sigma > 0.0: + print( + "No Rcl supplied, calc. Rcl assuming Sigma = %.2e g/cm^2" % Sigma + ) + Rcl = calcRcl(Mcl, Sigma) + else: + print("Error: must supple Rcl (in pc) or Sigma > 0") + + rho_ave = calcRhoAve(Mcl, Rcl=Rcl, Sigma=1) + + tff = (3 * math.pi / (32 * G * rho_ave)) ** 0.5 / SECYR + print("tff = %.2f yr" % tff) + return tff + + +def calc_Npart(Mcl, massRes=1e-3): + """ + Inputs: + #cloud mass = Mcl [Msun] + #mass/cell [Msun] + + Output: + N [particles] + """ + + N = Mcl / massRes + print("For Mcl = %.2e Msun, Mcell = %.2e" % (Mcl, massRes)) + print("N = %.2e particles" % N) + return N