From 266f5fb66dac86f4deae6ec746c54f74bd5bfc12 Mon Sep 17 00:00:00 2001 From: Marta Reina Campos Date: Tue, 24 Sep 2024 16:35:43 -0400 Subject: [PATCH 1/9] Added an input parameter to change the slope of the turbulent power spectrum --- MakeCloud.py | 699 ++++++++++++++++++++------------------------------- 1 file changed, 268 insertions(+), 431 deletions(-) diff --git a/MakeCloud.py b/MakeCloud.py index 92f1045..cdad59b 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -6,15 +6,15 @@ Options: -h --help Show this screen. - --R= Outer radius of the cloud in pc [default: 10.0] - --M= Mass of the cloud in msun [default: 2e4] + --R= Outer radius of the cloud in pc [default: 20.0] + --M= Mass of the cloud in msun [default: 1e5] --filename= Name of the IC file to be generated - --N= Number of gas particles [default: 2000000] + --N= Number of gas particles [default: 125000] --density_exponent= Power law exponent of the density profile [default: 0.0] --spin= Spin parameter: fraction of binding energy in solid-body rotation [default: 0.0] --omega_exponent= Powerlaw exponent of rotational frequency as a function of cylindrical radius [default: 0.0] - --turb_type= Type of initial turbulent velocity (and possibly density field): 'gaussian' or 'full' [default: gaussian] - --turb_sol= Fraction of turbulence in solenoidal modes, used when turb_type is 'gaussian' [default: 0.5] + --turb_slope= Slope of the turbulent power spectra [default: 2.0] + --turb_sol= Fraction of turbulence in solenoidal modes [default: 0.5] --alpha_turb= Turbulent virial parameter (BM92 convention: 2Eturb/|Egrav|) [default: 2.] --bturb= Magnetic energy as a fraction of the binding energy [default: 0.1] --bfixed= Magnetic field in magnitude in code units, used instead of bturb if not set to zero [default: 0] @@ -49,98 +49,84 @@ --Z= Metallicity of the cloud in Solar units (just for params file) [default: 1.0] --ISRF= Interstellar radiation background of the cloud in Solar neighborhood units (just for params file) [default: 1.0] """ -# Example: python MakeCloud.py --M=1000 --N=1e7 --R=1.0 --localdir --param_only +#Example: python MakeCloud.py --M=1000 --N=1e7 --R=1.0 --localdir --param_only + -import os import numpy as np -from scipy import fftpack, interpolate +from scipy import fftpack, interpolate, ndimage +from scipy.integrate import quad, odeint, solve_bvp +from scipy.spatial import cKDTree from scipy.spatial.distance import cdist +from scipy.optimize import minimize +import sys import h5py +import os from docopt import docopt - def get_glass_coords(N_gas, glass_path): x = np.load(glass_path) Nx = len(x) - while len(x) * np.pi * 4 / 3 / 8 < N_gas: - print( - "Need %d particles, have %d. Tessellating 8 copies of the glass file to get required particle number" - % (N_gas * 8 / (4 * np.pi / 3), len(x)) - ) - x = np.concatenate( - [ - x / 2 - + i * np.array([0.5, 0, 0]) - + j * np.array([0, 0.5, 0]) - + k * np.array([0, 0, 0.5]) - for i in range(2) - for j in range(2) - for k in range(2) - ] - ) + while len(x)*np.pi*4/3 / 8 < N_gas: + print("Need %d particles, have %d. Tessellating 8 copies of the glass file to get required particle number"%(N_gas * 8 /(4*np.pi/3), len(x))) + x = np.concatenate([x/2 + i * np.array([0.5,0,0]) + j * np.array([0,0.5,0]) + k * np.array([0, 0, 0.5]) for i in range(2) for j in range(2) for k in range(2)]) Nx = len(x) print("Glass loaded!") return x - -def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): +def TurbField(res=256, minmode=2, maxmode = 64, slope = 2.0, sol_weight=1., seed=42): freqs = fftpack.fftfreq(res) - freq3d = np.array(np.meshgrid(freqs, freqs, freqs, indexing="ij")) - intfreq = np.around(freq3d * res) - kSqr = np.sum(np.abs(freq3d) ** 2, axis=0) - intkSqr = np.sum(np.abs(intfreq) ** 2, axis=0) + freq3d = np.array(np.meshgrid(freqs,freqs,freqs,indexing='ij')) + intfreq = np.around(freq3d*res) + kSqr = np.sum(np.abs(freq3d)**slope,axis=0) + intkSqr = np.sum(np.abs(intfreq)**2, axis=0) VK = [] + vSqr = 0.0 # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field for i in range(3): - np.random.seed(seed + i) - rand_phase = fftpack.fftn( - np.random.normal(size=kSqr.shape) - ) # fourier transform of white noise - vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300) - vk[intkSqr == 0] = 0.0 - vk[intkSqr < minmode**2] *= ( - intkSqr[intkSqr < minmode**2] ** 2 / minmode**4 - ) # smoother filter than mode-freezing; should give less "ringing" artifacts - vk *= np.exp(-intkSqr / maxmode**2) + np.random.seed(seed+i) + rand_phase = fftpack.fftn(np.random.normal(size=kSqr.shape)) # fourier transform of white noise + vk = rand_phase * (float(minmode)/res)**2 / (kSqr+1e-300) + #vk[intkSqr < minmode**2] = 0.0 # freeze out modes lower than minmode + vk[intkSqr==0] = 0.0 +# vk[intkSqr>0] *= np.exp(-minmode**2/intkSqr) + vk[intkSqr < minmode**2] *= intkSqr[intkSqr < minmode**2]**2/minmode**4 # smoother filter than mode-freezing; should give less "ringing" artifacts + vk *= np.exp(-intkSqr/maxmode**2) VK.append(vk) VK = np.array(VK) - + vk_new = np.zeros_like(VK) - + # do projection operator to get the correct mix of compressive and solenoidal for i in range(3): for j in range(3): - if i == j: + 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] - ) - vk_new[:, kSqr == 0] = 0.0 + vk_new[i] += (1 - 2 * sol_weight) * freq3d[i]*freq3d[j]/(kSqr+1e-300) * VK[j] + vk_new[:,kSqr==0] = 0.0 VK = vk_new - - 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.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 return np.array(vel) arguments = docopt(__doc__) R = float(arguments["--R"]) M_gas = float(arguments["--M"]) -N_gas = int(float(arguments["--N"]) + 0.5) +N_gas = int(float(arguments["--N"])+0.5) M_star = float(arguments["--Mstar"]) -v_star = np.array([float(v) for v in arguments["--v_star"].split(",")]) +v_star = np.array([float(v) for v in arguments["--v_star"].split(',')]) spin = float(arguments["--spin"]) omega_exponent = float(arguments["--omega_exponent"]) -turbulence = float(arguments["--alpha_turb"]) / 2 -seed = int(float(arguments["--turb_seed"]) + 0.5) +turbulence = float(arguments["--alpha_turb"])/2 +seed = int(float(arguments["--turb_seed"])+0.5) tmax = int(float(arguments["--tmax"])) nsnap = int(float(arguments["--nsnap"])) +turb_slope = float(arguments["--turb_slope"]) turb_sol = float(arguments["--turb_sol"]) magnetic_field = float(arguments["--bturb"]) bfixed = float(arguments["--bfixed"]) @@ -153,7 +139,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): length_unit = float(arguments["--length_unit"]) mass_unit = float(arguments["--mass_unit"]) v_unit = float(arguments["--v_unit"]) -t_unit = length_unit / v_unit +t_unit = length_unit/v_unit G = 4300.71 * v_unit**-2 * mass_unit / length_unit makebox = arguments["--makebox"] impact_param = float(arguments["--impact_param"]) @@ -161,27 +147,20 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): v_impact = float(arguments["--v_impact"]) impact_axis = arguments["--impact_axis"] makecylinder = arguments["--makecylinder"] -cyl_aspect_ratio = float(arguments["--cyl_aspect_ratio"]) -fixed_ncrit = float(arguments["--fixed_ncrit"]) -density_exponent = float(arguments["--density_exponent"]) -metallicity = float(arguments["--Z"]) +cyl_aspect_ratio=float(arguments["--cyl_aspect_ratio"]) +fixed_ncrit=float(arguments["--fixed_ncrit"]) +density_exponent=float(arguments["--density_exponent"]) +metallicity=float(arguments["--Z"]) ISRF = float(arguments["--ISRF"]) -if arguments["--turb_path"]: - turb_path = arguments["--turb_path"] -else: - turb_path = os.path.expanduser("~") + "/turb" -if arguments["--glass_path"]: - glass_path = arguments["--glass_path"] +if arguments["--turb_path"]: turb_path = arguments["--turb_path"] +else: turb_path = os.path.expanduser('~') + "/turb" +if arguments["--glass_path"]: glass_path = arguments["--glass_path"] else: - glass_path = os.path.expanduser("~") + "/glass_orig.npy" + glass_path = os.path.expanduser('~') + "/glass_orig.npy" if not os.path.exists(glass_path): import urllib.request - 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 - ) + urllib.request.urlretrieve("http://www.tapir.caltech.edu/~mgrudich/glass_orig.npy",glass_path) if localdir: turb_path = "turb" @@ -190,490 +169,348 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): if arguments["--boxsize"] is not None: boxsize = float(arguments["--boxsize"]) else: - boxsize = 10 * R + boxsize = 10*R if arguments["--x_star"]: - x_star = np.array([float(x) for x in arguments["--x_star"].split(",")]) -else: # default to center of box - x_star = np.repeat(0.5 * boxsize, 3) + x_star = np.array([float(x) for x in arguments["--x_star"].split(',')]) +else: # default to center of box + x_star = np.repeat(0.5*boxsize,3) derefinement = arguments["--derefinement"] -res_effective = int(N_gas ** (1.0 / 3.0) + 0.5) -phimode = float(arguments["--phimode"]) - -filename = ( - "M%3.2g_" % (M_gas) - + ("Mstar%g_" % (M_star) if M_star > 0 else "") - + ("rho_exp%g_" % (-density_exponent) if density_exponent < 0 else "") - + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_sol%g" - % ( - R, - metallicity, - spin, - 2 * turbulence, - magnetic_field, - ISRF, - res_effective, - minmode, - turb_sol, - ) - + ("_%d" % seed) - + ( - "_collision_%g_%g_%g_%s" % (impact_dist, impact_param, v_impact, impact_axis) - if impact_dist > 0 - else "" - ) - + ".hdf5" -) -filename = filename.replace("+", "").replace("e0", "e") -filename = "".join(filename.split()) +res_effective = int(N_gas**(1.0/3.0)+0.5) +phimode=float(arguments["--phimode"]) -delta_m = M_gas / N_gas +filename = "M%3.2g_"%(M_gas) + ("Mstar%g_"%(M_star) if M_star>0 else "") + ("rho_exp%g_"%(-density_exponent) if density_exponent<0 else "") \ + + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_beta%g_sol%g"%(R,metallicity,spin,2*turbulence,magnetic_field,ISRF,res_effective,minmode,turb_slope,turb_sol) \ + + ("_%d"%seed) + ("_collision_%g_%g_%g_%s"%(impact_dist,impact_param,v_impact,impact_axis) if impact_dist>0 else "") + ".hdf5" +filename = filename.replace("+","").replace('e0','e') +filename = "".join(filename.split()) + +delta_m = M_gas/N_gas 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) - ncrit = 1e13 # ~100x the opacity limit -else: # something more FIRE-like, where we rely on a sub-grid prescription turning gas into star particles +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) + 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 ncrit = 100 if fixed_ncrit: - ncrit = fixed_ncrit - -tff = (3 * np.pi / (32 * G * rho_avg)) ** 0.5 -L = (4 * np.pi * R**3 / 3) ** (1.0 / 3) # volume-equivalent box size -vrms = (6 / 5 * G * M_gas / R) ** 0.5 * turbulence**0.5 + ncrit=fixed_ncrit + +tff = (3*np.pi/(32*G*rho_avg))**0.5 +L = (4*np.pi*R**3/3)**(1./3) # volume-equivalent box size +vrms = (6/5 * G * M_gas / R)**0.5 * turbulence**0.5 if turbulence: - tcross = L / vrms + tcross = L/vrms else: tcross = tff -turbenergy = ( - 0.019111097819633344 * vrms**3 / L -) # 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() -) - -jet_particle_mass = min(delta_m, max(1e-4, delta_m / 10.0)) -MS_wind_particle_mass = ( - jet_particle_mass / 10 -) # MS winds have lower mdot than jets, so we should be able to better resolve them this way - -replacements = { - "NAME": filename.replace(".hdf5", ""), - "DTSNAP": tff / nsnap, - "MAXTIMESTEP": tff / (nsnap), - "SOFTENING": softening, - "GASSOFT": 2.0e-8, - "TMAX": tff * tmax, - "RHOMAX": ncrit, - "BOXSIZE": boxsize, - "OUTFOLDER": "output", - "JET_PART_MASS": jet_particle_mass, - "MS_WIND_PART_MASS": MS_wind_particle_mass, - "BH_SEED_MASS": delta_m / 2.0, - "TURBDECAY": tcross / 2, - "TURBENERGY": turbenergy, - "TURBFREQ": tcross / 20, - "TURB_KMIN": int(100 * 2 * np.pi / L) / 100.0, - "TURB_KMAX": int(100 * 4 * np.pi / (L) + 1) / 100.0, - "TURB_SIGMA": vrms, - "TURB_MINLAMBDA": int(100 * R / 2) / 100, - "TURB_MAXLAMBDA": int(100 * R * 2) / 100, - "TURB_COHERENCE_TIME": tcross / 2, - "UNIT_L": 3.085678e18 * length_unit, - "UNIT_M": 1.989e33 * mass_unit, - "UNIT_V": v_unit * 1e2, - "UNIT_B": B_unit, - "ZINIT": metallicity, - "ISRF": ISRF, -} - -for k, r in replacements.items(): - paramsfile = paramsfile.replace(k,(r if isinstance(r,str) else '{:.2e}'.format(r))) - -open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write(paramsfile) +turbenergy = 0.019111097819633344*vrms**3/L # 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()) + +jet_particle_mass = min(delta_m,max(1e-4, delta_m/10.0)) +MS_wind_particle_mass = jet_particle_mass / 10 #MS winds have lower mdot than jets, so we should be able to better resolve them this way + +replacements = {"NAME": filename.replace(".hdf5",""), "DTSNAP": tff/nsnap, "MAXTIMESTEP": tff/(nsnap), + "SOFTENING": softening, "GASSOFT": 2.0e-8, "TMAX": tff*tmax, "RHOMAX": ncrit, + "BOXSIZE": boxsize, "OUTFOLDER": "output", "JET_PART_MASS": jet_particle_mass, + "MS_WIND_PART_MASS": MS_wind_particle_mass, "BH_SEED_MASS": delta_m/2.0 , "TURBDECAY": tcross/2, + "TURBENERGY": turbenergy, "TURBFREQ": tcross/20, "TURB_KMIN": int(100 * 2*np.pi/L)/100., + "TURB_KMAX": int(100*4*np.pi/(L)+1)/100., "TURB_SIGMA": vrms, "TURB_MINLAMBDA": int(100*R/2)/100, + "TURB_MAXLAMBDA": int(100*R*2)/100, "TURB_COHERENCE_TIME": tcross/2, "UNIT_L": 3.085678e18*length_unit, + "UNIT_M": 1.989e33*mass_unit, "UNIT_V": v_unit*1e2, "UNIT_B": B_unit, "ZINIT": metallicity, "ISRF": ISRF} + +for k in replacements.keys(): + paramsfile = paramsfile.replace(k, str(replacements[k])) +open("params_"+filename.replace(".hdf5","")+".txt", "w").write(paramsfile) if makebox: replacements_box = replacements.copy() - replacements_box["NAME"] = filename.replace(".hdf5", "_BOX") + replacements_box["NAME"] = filename.replace(".hdf5","_BOX") replacements_box["BOXSIZE"] = L - replacements_box["TURB_MINLAMBDA"] = int(100 * L / 2) / 100 - replacements_box["TURB_MAXLAMBDA"] = int(100 * L * 2) / 100 - paramsfile = str( - open( - os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r" - ).read() - ) + replacements_box["TURB_MINLAMBDA"] = int(100*L/2)/100; replacements_box["TURB_MAXLAMBDA"] = int(100*L*2)/100; + paramsfile = str(open(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) + paramsfile = paramsfile.replace(k, str(replacements_box[k])) + open("params_"+filename.replace(".hdf5","")+"_BOX.txt", "w").write(paramsfile) if makecylinder: - # Get cylinder params - R_cyl = R * np.sqrt( - np.pi / (4 * cyl_aspect_ratio) - ) # surface density equivalent cylinder - L_cyl = R_cyl * 2 * cyl_aspect_ratio - vrms_cyl = ( - 2 * G * M_gas / L_cyl - ) ** 0.5 * turbulence**0.5 # the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L - vrms_cyl *= 0.71 # additional scaling found numerically to make the stirring run reproduce the right alpha and filament length (similarly determined numerical factor added to GIZMO) - tcross_cyl = 2 * R_cyl / vrms_cyl - boxsize_cyl = ( - L_cyl * 1.5 + R_cyl * 5 - ) # the box should fit the cylinder and be many times bigger than its width - print( - "Cylinder params: L=%g R=%g boxsize=%g vrms=%g" - % (L_cyl, R_cyl, boxsize_cyl, vrms_cyl) - ) + #Get cylinder params + R_cyl = R * np.sqrt( np.pi/(4*cyl_aspect_ratio) ) #surface density equivalent cylinder + L_cyl = R_cyl*2*cyl_aspect_ratio + vrms_cyl = (2 * G * M_gas / L_cyl)**0.5 * turbulence**0.5 #the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L + vrms_cyl *= 0.71 #additional scaling found numerically to make the stirring run reproduce the right alpha and filament length (similarly determined numerical factor added to GIZMO) + tcross_cyl = 2*R_cyl/vrms_cyl + boxsize_cyl = L_cyl*1.5+R_cyl*5 #the box should fit the cylinder and be many times bigger than its width + print("Cylinder params: L=%g R=%g boxsize=%g vrms=%g"%(L_cyl,R_cyl,boxsize_cyl,vrms_cyl)) replacements_cyl = replacements.copy() - replacements_cyl["NAME"] = filename.replace(".hdf5", "_CYL") + replacements_cyl["NAME"] = filename.replace(".hdf5","_CYL") replacements_cyl["BOXSIZE"] = boxsize_cyl - # New driving params - replacements_cyl["TURB_MINLAMBDA"] = int(100 * R_cyl) / 100 - replacements_cyl["TURB_MAXLAMBDA"] = int(100 * R_cyl * 4) / 100 - replacements_cyl["TURB_SIGMA"] = vrms_cyl - replacements_cyl["TURB_COHERENCE_TIME"] = tcross_cyl / 2 - # Legacy driving params, probably needs tuning - replacements_cyl["TURBDECAY"] = tcross_cyl / 2 - replacements_cyl["TURBENERGY"] = 0.019111097819633344 * vrms_cyl**3 / R_cyl - replacements_cyl["TURBFREQ"] = tcross_cyl / 20 - replacements_cyl["TURB_KMIN"] = int(100 * 2 * np.pi / R_cyl) / 100.0 - 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" - ).read() - ) + #New driving params + replacements_cyl["TURB_MINLAMBDA"] = int(100*R_cyl)/100; replacements_cyl["TURB_MAXLAMBDA"] = int(100*R_cyl*4)/100; + replacements_cyl["TURB_SIGMA"] = vrms_cyl; replacements_cyl["TURB_COHERENCE_TIME"] = tcross_cyl/2; + #Legacy driving params, probably needs tuning + replacements_cyl["TURBDECAY"] = tcross_cyl/2; replacements_cyl["TURBENERGY"] = 0.019111097819633344*vrms_cyl**3/R_cyl; replacements_cyl["TURBFREQ"] = tcross_cyl/20; + replacements_cyl["TURB_KMIN"] = int(100 * 2*np.pi/R_cyl)/100.; replacements_cyl["TURB_KMAX"] = int(100*4*np.pi/(R_cyl)+1)/100.; + paramsfile = str(open(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) - + paramsfile = paramsfile.replace(k, str(replacements_cyl[k])) + open("params_"+filename.replace(".hdf5","")+"_CYL.txt", "w").write(paramsfile) + if param_only: - print("Parameters only run, exiting...") + print('Parameters only run, exiting...') exit() -dm = M_gas / N_gas +dm = M_gas/N_gas mgas = np.repeat(dm, N_gas) x = get_glass_coords(N_gas, glass_path) Nx = len(x) -x = 2 * (x - 0.5) +x = 2*(x-0.5) print("Computing radii...") -r = cdist(x, [np.zeros(3)])[:, 0] +r = cdist(x, [np.zeros(3)])[:,0] print("Done! Sorting coordinates...") x = x[r.argsort()][:N_gas] print("Done! Rescaling...") -x *= (float(Nx) / N_gas * 4 * np.pi / 3 / 8) ** (1.0 / 3) * R +x *= (float(Nx) / N_gas * 4*np.pi/3 / 8)**(1./3)*R print("Done! Recomupting radii...") -r = cdist(x, [np.zeros(3)])[:, 0] -x, r = x / r.max(), r / r.max() +r = cdist(x, [np.zeros(3)])[:,0] +x, r = x/r.max(), r/r.max() print("Doing density profile...") -rnew = r ** (3.0 / (3 + density_exponent)) * R -x = x * (rnew / r)[:, None] -r = np.sum(x**2, axis=1) ** 0.5 +rnew = r**(3. / (3+density_exponent)) * R +x=x * (rnew/r)[:,None] +r = np.sum(x**2, axis=1)**0.5 r_order = r.argsort() -x, r = np.take(x, r_order, axis=0), r[r_order] +x, r = np.take(x,r_order,axis=0), r[r_order] -if not os.path.exists(turb_path): - os.makedirs(turb_path) -fname = turb_path + "/vturb%d_sol%g_seed%d.npy" % (minmode, turb_sol, seed) +if not os.path.exists(turb_path): os.makedirs(turb_path) +fname = turb_path + "/vturb%d_beta%g_sol%g_seed%d.npy"%(minmode,turb_slope,turb_sol, seed) if not os.path.isfile(fname): - vt = TurbField(minmode=minmode, sol_weight=turb_sol, seed=seed) - nmin, nmax = vt.shape[-1] // 4, 3 * vt.shape[-1] // 4 - vt = vt[ - :, nmin:nmax, nmin:nmax, nmin:nmax - ] # we take the central cube of size L/2 so that opposide sides of the cloud are not correlated + vt = TurbField(minmode=minmode, slope = turb_slope, sol_weight=turb_sol, seed=seed) + nmin, nmax = vt.shape[-1]// 4, 3*vt.shape[-1]//4 + vt = vt[:,nmin:nmax, nmin:nmax, nmin:nmax] # we take the central cube of size L/2 so that opposide sides of the cloud are not correlated np.save(fname, vt) else: vt = np.load(fname) - -xgrid = np.linspace(-R, R, vt.shape[-1]) + +xgrid = np.linspace(-R,R,vt.shape[-1]) v = [] for i in range(3): - v.append(interpolate.interpn((xgrid, xgrid, xgrid), vt[i, :, :, :], x)) + v.append(interpolate.interpn((xgrid,xgrid,xgrid),vt[i,:,:,:],x)) v = np.array(v).T print("Coordinates obtained!") - + Mr = mgas.cumsum() -ugrav = G * np.sum(Mr / r * mgas) -v -= np.average(v, axis=0) -Eturb = 0.5 * M_gas / N_gas * np.sum(v**2) -v *= np.sqrt(turbulence * ugrav / Eturb) +ugrav = G * np.sum(Mr/ r * mgas) +v -= np.average(v,axis=0) +Eturb = 0.5*M_gas/N_gas*np.sum(v**2) +v *= np.sqrt(turbulence*ugrav/Eturb) E_rot_target = spin * ugrav -Rcyl = np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2) +Rcyl = np.sqrt(x[:,0]**2+x[:,1]**2) omega = Rcyl**omega_exponent -vrot = np.cross(np.c_[np.zeros_like(omega), np.zeros_like(omega), omega], x) -Erot_actual = np.sum(0.5 * mgas[:, None] * vrot**2) -vrot *= np.sqrt(E_rot_target / Erot_actual) +vrot = np.cross(np.c_[np.zeros_like(omega),np.zeros_like(omega),omega], x) +Erot_actual = np.sum(0.5*mgas[:,None]*vrot**2) +vrot *= np.sqrt(E_rot_target/Erot_actual) v += vrot 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 -) # alfven speed for unit magnetic field -uB = 0.5 * M_gas * vA_unit**2 # magnetic energy we would have for unit magnetic field -if bfixed > 0: +vA_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 +if (bfixed>0): B = B * bfixed else: - B = B * np.sqrt( - magnetic_field * ugrav / uB - ) # renormalize to desired magnetic energy + B = B * np.sqrt(magnetic_field*ugrav/uB) # renormalize to desired magnetic energy v = v - np.average(v, axis=0) x = x - np.average(x, axis=0) -r, phi = np.sum(x**2, axis=1) ** 0.5, np.arctan2(x[:, 1], x[:, 0]) -theta = np.arccos(x[:, 2] / r) -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)] -) +r, phi = np.sum(x**2,axis=1)**0.5, np.arctan2(x[:,1],x[:,0]) +theta = np.arccos(x[:,2]/r) +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)] if makecylinder: - - def ind_in_cylinder(x, L_cyl, R_cyl): - return (np.abs(x[:, 0]) < L_cyl / 2) & ( - np.sum(x[:, 1:] ** 2, axis=1) < R_cyl**2 - ) - - # Just get a roughly homogeneous cylinder along the x axis, we will stir it anyway - N_cyl = 0 - while ( - N_cyl <= N_gas - ): # should be very unlikely that we need to repeat, but let's check to be sure - x_cyl = np.random.rand(2 * N_gas, 3) * 2 - 1 - x_cyl[:, 0] *= L_cyl / 2 - x_cyl[:, 1] *= R_cyl - x_cyl[:, 2] *= R_cyl - x_cyl = x_cyl[ind_in_cylinder(x_cyl, L_cyl, R_cyl)] + def ind_in_cylinder(x,L_cyl,R_cyl): + return ( (np.abs(x[:,0]) < L_cyl/2) & ( np.sum(x[:,1:]**2,axis=1) 0: - x = np.concatenate([x, x]) - impact_dir = { - "x": np.array([1.0, 0, 0]), - "y": np.array([0, 1, 0]), - "z": np.array([0, 0, 1]), - }[impact_axis] - impact_param_dir = { - "x": np.array([0, 1, 0]), - "y": np.array([0, 0, 1]), - "z": np.array([1, 0, 0]), - }[impact_axis] +u = np.ones_like(mgas)*0.101/2.0 #/2 needed because it is molecular + +if impact_dist > 0: + x = np.concatenate([x,x]) + impact_dir = {"x": np.array([1.,0,0]), "y": np.array([0,1,0]), 'z': np.array([0,0,1])}[impact_axis] + impact_param_dir = {"x": np.array([0,1,0]), "y": np.array([0,0,1]), 'z': np.array([1,0,0])}[impact_axis] x[:N_gas] += impact_dist * R * impact_dir x[N_gas:] -= impact_dist * R * impact_dir x[:N_gas] += 0.5 * impact_param * R * impact_param_dir x[N_gas:] -= 0.5 * impact_param * R * impact_param_dir - v = np.concatenate([v, v]) - vrms = np.sum(v**2, axis=1).mean() ** 0.5 + v = np.concatenate([v,v]) + vrms = np.sum(v**2,axis=1).mean()**0.5 v[:N_gas] -= v_impact * vrms * impact_dir v[N_gas:] += v_impact * vrms * impact_dir - B = np.concatenate([B, B]) - u = np.concatenate([u, u]) - mgas = np.concatenate([mgas, mgas]) + B = np.concatenate([B,B]) + u = np.concatenate([u,u]) + mgas = np.concatenate([mgas,mgas]) -u = ( - np.ones_like(mgas) * (200 / v_unit) ** 2 -) # start with specific internal energy of (200m/s)^2, this is overwritten unless starting with restart flag 2###### #0.101/2.0 #/2 needed because it is molecular +u = np.ones_like(mgas)*(200/v_unit)**2 #start with specific internal energy of (200m/s)^2, this is overwritten unless starting with restart flag 2###### #0.101/2.0 #/2 needed because it is molecular if diffuse_gas: # assuming 10K vs 10^4K gas: factor of ~10^3 density contrast - rho_warm = M_gas * 3 / (4 * np.pi * R**3) / 1000 - M_warm = ( - boxsize**3 - (4 * np.pi * R**3 / 3) - ) * rho_warm # mass of diffuse box-filling medium - N_warm = int(M_warm / (M_gas / N_gas)) + rho_warm = M_gas*3/(4*np.pi*R**3) / 1000 + M_warm = (boxsize**3 - (4*np.pi*R**3 / 3)) * rho_warm # mass of diffuse box-filling medium + N_warm = int(M_warm/(M_gas/N_gas)) if derefinement: x0 = get_glass_coords(N_gas, glass_path) Nx = len(x0) - x0 = 2 * (x0 - 0.5) - r0 = (x0 * x0).sum(1) ** 0.5 + x0 = 2*(x0-0.5) + r0 = (x0*x0).sum(1)**0.5 x0, r0 = x0[r0.argsort()], r0[r0.argsort()] # first lay down the stuff within 3*R - N_warm = int( - 4 * np.pi * rho_warm * (3 * R) ** 3 / 3 / dm - ) # number of cells within 3R - x_warm = ( - x0[:N_warm] * 3 * R / r0[N_warm - 1] - ) # uniform density of cells within 3R - x0 = x0[ - N_warm: - ] # now we take the ones outside the initial sphere and map them to a n(R) ~ R^-3 profile so that we get constant number of cells per log radius interval + N_warm = int(4*np.pi*rho_warm*(3*R)**3/3/dm) # number of cells within 3R + print(N_warm) + x_warm = x0[:N_warm] * 3 * R / r0[N_warm-1] # uniform density of cells within 3R + x0 = x0[N_warm:] # now we take the ones outside the initial sphere and map them to a n(R) ~ R^-3 profile so that we get constant number of cells per log radius interval r0 = r0[N_warm:] - rnew = 3 * R * np.exp(np.arange(len(x0)) / N_warm / 3) - x_warm = np.concatenate([x_warm, (rnew / r0)[:, None] * x0], axis=0) - x_warm = x_warm[np.max(np.abs(x_warm), axis=1) < boxsize / 2] + rnew = 3*R*np.exp(np.arange(len(x0))/N_warm/3) + x_warm = np.concatenate([x_warm, (rnew/r0)[:,None] * x0],axis=0) + 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)]) + 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)]) 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] + 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 - B = np.concatenate( - [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([0, 0, 1])] - ) - u = np.concatenate([u, np.repeat(101.0, N_warm)]) + v = np.concatenate([v, np.zeros((N_warm,3))]) + Bmag = np.average(np.sum(B**2,axis=1))**0.5 + B = np.concatenate([B, np.repeat(Bmag,N_warm)[:,np.newaxis] * np.array([0,0,1])]) + u = np.concatenate([u, np.repeat(101.,N_warm)]) if makecylinder: - # The magnetic field is paralell to the cylinder (true at low densities, so probably fine for IC) - B_cyl = np.concatenate( - [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([1, 0, 0])] - ) - # Add diffuse medium - M_warm_cyl = (boxsize_cyl**3 - (4 * np.pi * R**3 / 3)) * rho_warm - N_warm_cyl = int(M_warm_cyl / (M_gas / N_gas)) - x_warm = ( - boxsize_cyl * np.random.rand(N_warm_cyl, 3) - boxsize_cyl / 2 - ) # will be recentered later - x_warm = x_warm[ - ~ind_in_cylinder(x_warm, L_cyl, R_cyl) - ] # keep only warm gas outside the cylinder - # print("N_warm_cyl: %g N_warm_cyl_kept %g "%(N_warm_cyl,len(x_warm))) + #The magnetic field is paralell to the cylinder (true at low densities, so probably fine for IC) + B_cyl = np.concatenate([B, np.repeat(Bmag,N_warm)[:,np.newaxis] * np.array([1,0,0])]) + #Add diffuse medium + M_warm_cyl = (boxsize_cyl**3 - (4*np.pi*R**3 / 3)) * rho_warm + N_warm_cyl = int(M_warm_cyl/(M_gas/N_gas)) + x_warm = boxsize_cyl*np.random.rand(N_warm_cyl, 3) - boxsize_cyl/2 #will be recentered later + x_warm = x_warm[ ~ind_in_cylinder(x_warm,L_cyl,R_cyl) ] #keep only warm gas outside the cylinder + #print("N_warm_cyl: %g N_warm_cyl_kept %g "%(N_warm_cyl,len(x_warm))) N_warm_cyl = len(x_warm) x_cyl = np.concatenate([x_cyl, x_warm]) - v_cyl = np.concatenate([v_cyl, np.zeros((N_warm, 3))]) + v_cyl = np.concatenate([v_cyl, np.zeros((N_warm,3))]) else: N_warm = 0 -rho = np.repeat(3 * M_gas / (4 * np.pi * R**3), len(mgas)) -if diffuse_gas: - rho[-N_warm:] /= 1000 -h = (32 * mgas / rho) ** (1.0 / 3) +rho = np.repeat(3*M_gas/(4*np.pi*R**3), len(mgas)) +if diffuse_gas: rho[-N_warm:] /= 1000 +h = (32*mgas/rho)**(1./3) -x += boxsize / 2 # cloud is always centered at (boxsize/2,boxsize/2,boxsize/2) +x += boxsize/2 # cloud is always centered at (boxsize/2,boxsize/2,boxsize/2) if makecylinder: - x_cyl += boxsize_cyl / 2 + x_cyl += boxsize_cyl/2 print("Writing snapshot...") -F = h5py.File(filename, "w") +F=h5py.File(filename, 'w') F.create_group("PartType0") F.create_group("Header") -F["Header"].attrs["NumPart_ThisFile"] = [ - 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["NumPart_ThisFile"] = [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) F["PartType0"].create_dataset("Coordinates", data=x) F["PartType0"].create_dataset("Velocities", data=v) -F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(len(mgas))) +F["PartType0"].create_dataset("ParticleIDs", data=1+np.arange(len(mgas))) F["PartType0"].create_dataset("InternalEnergy", data=u) -if M_star > 0: +if M_star>0: F.create_group("PartType5") - # Let's add the sink at the center + #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("Velocities", data=[v_star]) # at rest - F["PartType5"].create_dataset( - "ParticleIDs", data=np.array([F["PartType0/ParticleIDs"][:].max() + 1]) - ) - # Advanced properties for sinks - F["PartType5"].create_dataset( - "BH_Mass", data=M_star - ) # all the mass in the sink/protostar/star - F["PartType5"].create_dataset( - "BH_Mass_AlphaDisk", data=np.array([0.0]) - ) # starts with no disk - F["PartType5"].create_dataset( - "BH_Mdot", data=np.array([0.0]) - ) # starts with no mdot - F["PartType5"].create_dataset( - "BH_Specific_AngMom", data=np.array([0.0]) - ) # starts with no angular momentum - F["PartType5"].create_dataset( - "SinkRadius", data=np.array([softening]) - ) # Sinkradius set to softening - F["PartType5"].create_dataset("StellarFormationTime", data=np.array([0.0])) - F["PartType5"].create_dataset("ProtoStellarAge", data=np.array([0.0])) - F["PartType5"].create_dataset( - "ProtoStellarStage", data=np.array([5], dtype=np.int32), dtype=np.int32 - ) - # Stellar properties - # if (central_star or central_SN): - # if central_star: - # print("Assuming central sink is a ZAMS star") - # starts as ZAMS star - # else: - # print("Assuming central sink is a ZAMS star about to go supernova") - # F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([6],dtype=np.int32), dtype=np.int32) #starts as ZAMS star going SN - # Set guess for ZAMS stellar radius, will be overwritten - if (M_star) > 1.0: - R_ZAMS = (M_star) ** 0.57 + 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])) + #Advanced properties for sinks + F["PartType5"].create_dataset("BH_Mass", data=M_star) #all the mass in the sink/protostar/star + F["PartType5"].create_dataset("BH_Mass_AlphaDisk", data=np.array([0.])) #starts with no disk + F["PartType5"].create_dataset("BH_Mdot", data=np.array([0.])) #starts with no mdot + F["PartType5"].create_dataset("BH_Specific_AngMom", data=np.array([0.])) #starts with no angular momentum + F["PartType5"].create_dataset("SinkRadius", data=np.array([softening])) #Sinkradius set to softening + F["PartType5"].create_dataset("StellarFormationTime", data=np.array([0.])) + F["PartType5"].create_dataset("ProtoStellarAge", data=np.array([0.])) + F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([5],dtype=np.int32), dtype=np.int32) + #Stellar properties + #if (central_star or central_SN): + #if central_star: +# print("Assuming central sink is a ZAMS star") + #starts as ZAMS star + #else: + #print("Assuming central sink is a ZAMS star about to go supernova") + #F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([6],dtype=np.int32), dtype=np.int32) #starts as ZAMS star going SN + #Set guess for ZAMS stellar radius, will be overwritten + if ((M_star)>1.0): + R_ZAMS = (M_star)**0.57 else: - R_ZAMS = (M_star) ** 0.8 - 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("Mass_D", data=np.array([0.0])) # No D left - + R_ZAMS = (M_star)**0.8 + 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.])) #dummy + F["PartType5"].create_dataset("Mass_D", data=np.array([0.])) #No D left + if magnetic_field > 0.0: F["PartType0"].create_dataset("MagneticField", data=B) F.close() if makebox: - F = h5py.File(filename.replace(".hdf5", "_BOX.hdf5"), "w") + F=h5py.File(filename.replace(".hdf5","_BOX.hdf5"), 'w') F.create_group("PartType0") F.create_group("Header") - F["Header"].attrs["NumPart_ThisFile"] = [len(mgas), 0, 0, 0, 0, 0] - F["Header"].attrs["NumPart_Total"] = [len(mgas), 0, 0, 0, 0, 0] - F["Header"].attrs["MassTable"] = [M_gas / len(mgas), 0, 0, 0, 0, 0] - F["Header"].attrs["BoxSize"] = (4 * np.pi * R**3 / 3) ** (1.0 / 3) + F["Header"].attrs["NumPart_ThisFile"] = [len(mgas),0,0,0,0,0] + F["Header"].attrs["NumPart_Total"] = [len(mgas),0,0,0,0,0] + F["Header"].attrs["MassTable"] = [M_gas/len(mgas),0,0,0,0,0] + F["Header"].attrs["BoxSize"] = (4 * np.pi * R**3 / 3)**(1./3) 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"] - ) - F["PartType0"].create_dataset("Velocities", data=np.zeros((len(mgas), 3))) - F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(len(mgas))) + 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"]) + F["PartType0"].create_dataset("Velocities", data=np.zeros((len(mgas),3))) + F["PartType0"].create_dataset("ParticleIDs", data=1+np.arange(len(mgas))) F["PartType0"].create_dataset("InternalEnergy", data=u) if magnetic_field > 0.0: - F["PartType0"].create_dataset("MagneticField", data=B[: len(mgas)]) + F["PartType0"].create_dataset("MagneticField", data=B[:len(mgas)]) F.close() - + if makecylinder: - F = h5py.File(filename.replace(".hdf5", "_CYL.hdf5"), "w") + F=h5py.File(filename.replace(".hdf5","_CYL.hdf5"), 'w') F.create_group("PartType0") F.create_group("Header") - F["Header"].attrs["NumPart_ThisFile"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] - F["Header"].attrs["NumPart_Total"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] - F["Header"].attrs["MassTable"] = [M_gas / N_gas, 0, 0, 0, 0, 0] + F["Header"].attrs["NumPart_ThisFile"] = [N_gas+N_warm_cyl,0,0,0,0,0] + F["Header"].attrs["NumPart_Total"] = [N_gas+N_warm_cyl,0,0,0,0,0] + F["Header"].attrs["MassTable"] = [M_gas/N_gas,0,0,0,0,0] F["Header"].attrs["BoxSize"] = boxsize_cyl F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas) - F["PartType0"].create_dataset("Coordinates", data=x_cyl) + 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) F.close() + From 9952b0ba745fec6cb4ae66f1c0508f5278b75797 Mon Sep 17 00:00:00 2001 From: Marta Reina Campos Date: Tue, 24 Sep 2024 21:46:03 -0400 Subject: [PATCH 2/9] Fixed the change to the slope of the turbulent power spectra --- MakeCloud.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/MakeCloud.py b/MakeCloud.py index cdad59b..8b0684c 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -78,7 +78,7 @@ def TurbField(res=256, minmode=2, maxmode = 64, slope = 2.0, sol_weight=1., seed freqs = fftpack.fftfreq(res) freq3d = np.array(np.meshgrid(freqs,freqs,freqs,indexing='ij')) intfreq = np.around(freq3d*res) - kSqr = np.sum(np.abs(freq3d)**slope,axis=0) + kSqr = np.sum(np.abs(freq3d)**2,axis=0) intkSqr = np.sum(np.abs(intfreq)**2, axis=0) VK = [] @@ -87,7 +87,7 @@ def TurbField(res=256, minmode=2, maxmode = 64, slope = 2.0, sol_weight=1., seed for i in range(3): np.random.seed(seed+i) rand_phase = fftpack.fftn(np.random.normal(size=kSqr.shape)) # fourier transform of white noise - vk = rand_phase * (float(minmode)/res)**2 / (kSqr+1e-300) + vk = rand_phase * (float(minmode)/res)**2 / (np.power(kSqr, slope/2.0)+1e-300) #vk[intkSqr < minmode**2] = 0.0 # freeze out modes lower than minmode vk[intkSqr==0] = 0.0 # vk[intkSqr>0] *= np.exp(-minmode**2/intkSqr) From 4369f54bc300a200656753e01501eee48eba372f Mon Sep 17 00:00:00 2001 From: Marta Reina Campos Date: Tue, 24 Sep 2024 21:56:57 -0400 Subject: [PATCH 3/9] Changes had been made to an older version of MakeCloud, now this is the latest --- MakeCloud.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/MakeCloud.py b/MakeCloud.py index 92f1045..4c8f94a 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -13,7 +13,7 @@ --density_exponent= Power law exponent of the density profile [default: 0.0] --spin= Spin parameter: fraction of binding energy in solid-body rotation [default: 0.0] --omega_exponent= Powerlaw exponent of rotational frequency as a function of cylindrical radius [default: 0.0] - --turb_type= Type of initial turbulent velocity (and possibly density field): 'gaussian' or 'full' [default: gaussian] + --turb_slope= Slope of the turbulent power spectra [default: 2.0] --turb_sol= Fraction of turbulence in solenoidal modes, used when turb_type is 'gaussian' [default: 0.5] --alpha_turb= Turbulent virial parameter (BM92 convention: 2Eturb/|Egrav|) [default: 2.] --bturb= Magnetic energy as a fraction of the binding energy [default: 0.1] @@ -84,7 +84,7 @@ def get_glass_coords(N_gas, glass_path): return x -def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): +def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed=42): freqs = fftpack.fftfreq(res) freq3d = np.array(np.meshgrid(freqs, freqs, freqs, indexing="ij")) intfreq = np.around(freq3d * res) @@ -98,7 +98,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): rand_phase = fftpack.fftn( np.random.normal(size=kSqr.shape) ) # fourier transform of white noise - vk = rand_phase * (float(minmode) / res) ** 2 / (kSqr + 1e-300) + vk = rand_phase * (float(minmode) / res) ** 2 / (np.power(kSqr, slope/2.0) + 1e-300) vk[intkSqr == 0] = 0.0 vk[intkSqr < minmode**2] *= ( intkSqr[intkSqr < minmode**2] ** 2 / minmode**4 @@ -141,6 +141,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): seed = int(float(arguments["--turb_seed"]) + 0.5) tmax = int(float(arguments["--tmax"])) nsnap = int(float(arguments["--nsnap"])) +turb_slope = float(arguments["--turb_slope"]) turb_sol = float(arguments["--turb_sol"]) magnetic_field = float(arguments["--bturb"]) bfixed = float(arguments["--bfixed"]) @@ -206,7 +207,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): "M%3.2g_" % (M_gas) + ("Mstar%g_" % (M_star) if M_star > 0 else "") + ("rho_exp%g_" % (-density_exponent) if density_exponent < 0 else "") - + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_sol%g" + + "R%g_Z%g_S%g_A%g_B%g_I%g_Res%d_n%d_beta%g_sol%g" % ( R, metallicity, @@ -216,6 +217,7 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): ISRF, res_effective, minmode, + turb_slope, turb_sol, ) + ("_%d" % seed) @@ -381,9 +383,9 @@ def TurbField(res=256, minmode=2, maxmode=64, sol_weight=1.0, seed=42): if not os.path.exists(turb_path): os.makedirs(turb_path) -fname = turb_path + "/vturb%d_sol%g_seed%d.npy" % (minmode, turb_sol, seed) +fname = turb_path + "/vturb%d_beta%g_sol%g_seed%d.npy" % (minmode, turb_slope, turb_sol, seed) if not os.path.isfile(fname): - vt = TurbField(minmode=minmode, sol_weight=turb_sol, seed=seed) + vt = TurbField(minmode=minmode, slope = turb_slope, sol_weight=turb_sol, seed=seed) nmin, nmax = vt.shape[-1] // 4, 3 * vt.shape[-1] // 4 vt = vt[ :, nmin:nmax, nmin:nmax, nmin:nmax From 55eb65acd15e9ca47698c1317b3108e321e4a23e Mon Sep 17 00:00:00 2001 From: Marta Reina Campos Date: Tue, 24 Sep 2024 22:14:00 -0400 Subject: [PATCH 4/9] Now it should be the latest version --- MakeCloud.py | 624 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 381 insertions(+), 243 deletions(-) diff --git a/MakeCloud.py b/MakeCloud.py index 9cbdbb4..4fecc10 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -6,15 +6,15 @@ Options: -h --help Show this screen. - --R= Outer radius of the cloud in pc [default: 20.0] - --M= Mass of the cloud in msun [default: 1e5] + --R= Outer radius of the cloud in pc [default: 10.0] + --M= Mass of the cloud in msun [default: 2e4] --filename= Name of the IC file to be generated - --N= Number of gas particles [default: 125000] + --N= Number of gas particles [default: 2000000] --density_exponent= Power law exponent of the density profile [default: 0.0] --spin= Spin parameter: fraction of binding energy in solid-body rotation [default: 0.0] --omega_exponent= Powerlaw exponent of rotational frequency as a function of cylindrical radius [default: 0.0] --turb_slope= Slope of the turbulent power spectra [default: 2.0] - --turb_sol= Fraction of turbulence in solenoidal modes, used when turb_type is 'gaussian' [default: 0.5] + --turb_sol= Fraction of turbulence in solenoidal modes [default: 0.5] --alpha_turb= Turbulent virial parameter (BM92 convention: 2Eturb/|Egrav|) [default: 2.] --bturb= Magnetic energy as a fraction of the binding energy [default: 0.1] --bfixed= Magnetic field in magnitude in code units, used instead of bturb if not set to zero [default: 0] @@ -49,40 +49,49 @@ --Z= Metallicity of the cloud in Solar units (just for params file) [default: 1.0] --ISRF= Interstellar radiation background of the cloud in Solar neighborhood units (just for params file) [default: 1.0] """ -#Example: python MakeCloud.py --M=1000 --N=1e7 --R=1.0 --localdir --param_only - +# Example: python MakeCloud.py --M=1000 --N=1e7 --R=1.0 --localdir --param_only +import os import numpy as np -from scipy import fftpack, interpolate, ndimage -from scipy.integrate import quad, odeint, solve_bvp -from scipy.spatial import cKDTree +from scipy import fftpack, interpolate from scipy.spatial.distance import cdist -from scipy.optimize import minimize -import sys import h5py -import os from docopt import docopt + def get_glass_coords(N_gas, glass_path): x = np.load(glass_path) Nx = len(x) - while len(x)*np.pi*4/3 / 8 < N_gas: - print("Need %d particles, have %d. Tessellating 8 copies of the glass file to get required particle number"%(N_gas * 8 /(4*np.pi/3), len(x))) - x = np.concatenate([x/2 + i * np.array([0.5,0,0]) + j * np.array([0,0.5,0]) + k * np.array([0, 0, 0.5]) for i in range(2) for j in range(2) for k in range(2)]) + while len(x) * np.pi * 4 / 3 / 8 < N_gas: + print( + "Need %d particles, have %d. Tessellating 8 copies of the glass file to get required particle number" + % (N_gas * 8 / (4 * np.pi / 3), len(x)) + ) + x = np.concatenate( + [ + x / 2 + + i * np.array([0.5, 0, 0]) + + j * np.array([0, 0.5, 0]) + + k * np.array([0, 0, 0.5]) + for i in range(2) + for j in range(2) + for k in range(2) + ] + ) Nx = len(x) print("Glass loaded!") return x -def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed=42): + +def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42): freqs = fftpack.fftfreq(res) - freq3d = np.array(np.meshgrid(freqs,freqs,freqs,indexing='ij')) - intfreq = np.around(freq3d*res) - kSqr = np.sum(np.abs(freq3d)**2,axis=0) - intkSqr = np.sum(np.abs(intfreq)**2, axis=0) + freq3d = np.array(np.meshgrid(freqs, freqs, freqs, indexing="ij")) + intfreq = np.around(freq3d * res) + kSqr = np.sum(np.abs(freq3d) ** 2, axis=0) + intkSqr = np.sum(np.abs(intfreq) ** 2, axis=0) VK = [] - vSqr = 0.0 # apply ~k^-2 exp(-k^2/kmax^2) filter to white noise to get x, y, and z components of velocity field for i in range(3): np.random.seed(seed + i) @@ -98,34 +107,38 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= VK.append(vk) VK = np.array(VK) - + vk_new = np.zeros_like(VK) - + # do projection operator to get the correct mix of compressive and solenoidal for i in range(3): for j in range(3): - if i==j: + 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] - vk_new[:,kSqr==0] = 0.0 + vk_new[i] += ( + (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j] + ) + vk_new[:, kSqr == 0] = 0.0 VK = vk_new - - 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.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 return np.array(vel) arguments = docopt(__doc__) R = float(arguments["--R"]) M_gas = float(arguments["--M"]) -N_gas = int(float(arguments["--N"])+0.5) +N_gas = int(float(arguments["--N"]) + 0.5) M_star = float(arguments["--Mstar"]) -v_star = np.array([float(v) for v in arguments["--v_star"].split(',')]) +v_star = np.array([float(v) for v in arguments["--v_star"].split(",")]) spin = float(arguments["--spin"]) omega_exponent = float(arguments["--omega_exponent"]) -turbulence = float(arguments["--alpha_turb"])/2 -seed = int(float(arguments["--turb_seed"])+0.5) +turbulence = float(arguments["--alpha_turb"]) / 2 +seed = int(float(arguments["--turb_seed"]) + 0.5) tmax = int(float(arguments["--tmax"])) nsnap = int(float(arguments["--nsnap"])) turb_slope = float(arguments["--turb_slope"]) @@ -141,7 +154,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= length_unit = float(arguments["--length_unit"]) mass_unit = float(arguments["--mass_unit"]) v_unit = float(arguments["--v_unit"]) -t_unit = length_unit/v_unit +t_unit = length_unit / v_unit G = 4300.71 * v_unit**-2 * mass_unit / length_unit makebox = arguments["--makebox"] impact_param = float(arguments["--impact_param"]) @@ -149,20 +162,27 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= v_impact = float(arguments["--v_impact"]) impact_axis = arguments["--impact_axis"] makecylinder = arguments["--makecylinder"] -cyl_aspect_ratio=float(arguments["--cyl_aspect_ratio"]) -fixed_ncrit=float(arguments["--fixed_ncrit"]) -density_exponent=float(arguments["--density_exponent"]) -metallicity=float(arguments["--Z"]) +cyl_aspect_ratio = float(arguments["--cyl_aspect_ratio"]) +fixed_ncrit = float(arguments["--fixed_ncrit"]) +density_exponent = float(arguments["--density_exponent"]) +metallicity = float(arguments["--Z"]) ISRF = float(arguments["--ISRF"]) -if arguments["--turb_path"]: turb_path = arguments["--turb_path"] -else: turb_path = os.path.expanduser('~') + "/turb" -if arguments["--glass_path"]: glass_path = arguments["--glass_path"] +if arguments["--turb_path"]: + turb_path = arguments["--turb_path"] else: - glass_path = os.path.expanduser('~') + "/glass_orig.npy" + turb_path = os.path.expanduser("~") + "/turb" +if arguments["--glass_path"]: + glass_path = arguments["--glass_path"] +else: + glass_path = os.path.expanduser("~") + "/glass_orig.npy" if not os.path.exists(glass_path): import urllib.request + print("Downloading glass file...") - urllib.request.urlretrieve("http://www.tapir.caltech.edu/~mgrudich/glass_orig.npy",glass_path) + 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 + ) if localdir: turb_path = "turb" @@ -171,12 +191,12 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= if arguments["--boxsize"] is not None: boxsize = float(arguments["--boxsize"]) else: - boxsize = 10*R + boxsize = 10 * R if arguments["--x_star"]: - x_star = np.array([float(x) for x in arguments["--x_star"].split(',')]) -else: # default to center of box - x_star = np.repeat(0.5*boxsize,3) + x_star = np.array([float(x) for x in arguments["--x_star"].split(",")]) +else: # default to center of box + x_star = np.repeat(0.5 * boxsize, 3) derefinement = arguments["--derefinement"] @@ -211,104 +231,155 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= filename = filename.replace("+", "").replace("e0", "e") filename = "".join(filename.split()) -delta_m = M_gas/N_gas +delta_m = M_gas / N_gas 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) - ncrit = 1e13 # ~100x the opacity limit -else: # something more FIRE-like, where we rely on a sub-grid prescription turning gas into star particles +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) + 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 ncrit = 100 if fixed_ncrit: - ncrit=fixed_ncrit - -tff = (3*np.pi/(32*G*rho_avg))**0.5 -L = (4*np.pi*R**3/3)**(1./3) # volume-equivalent box size -vrms = (6/5 * G * M_gas / R)**0.5 * turbulence**0.5 + ncrit = fixed_ncrit + +tff = (3 * np.pi / (32 * G * rho_avg)) ** 0.5 +L = (4 * np.pi * R**3 / 3) ** (1.0 / 3) # volume-equivalent box size +vrms = (6 / 5 * G * M_gas / R) ** 0.5 * turbulence**0.5 if turbulence: - tcross = L/vrms + tcross = L / vrms else: tcross = tff -turbenergy = 0.019111097819633344*vrms**3/L # 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()) +turbenergy = ( + 0.019111097819633344 * vrms**3 / L +) # ST_Energy sets the dissipation rate of SPECIFIC energy ~ v^2 / (L/v) ~ v^3/L -jet_particle_mass = min(delta_m,max(1e-4, delta_m/10.0)) -MS_wind_particle_mass = jet_particle_mass / 10 #MS winds have lower mdot than jets, so we should be able to better resolve them this way - -replacements = {"NAME": filename.replace(".hdf5",""), "DTSNAP": tff/nsnap, "MAXTIMESTEP": tff/(nsnap), - "SOFTENING": softening, "GASSOFT": 2.0e-8, "TMAX": tff*tmax, "RHOMAX": ncrit, - "BOXSIZE": boxsize, "OUTFOLDER": "output", "JET_PART_MASS": jet_particle_mass, - "MS_WIND_PART_MASS": MS_wind_particle_mass, "BH_SEED_MASS": delta_m/2.0 , "TURBDECAY": tcross/2, - "TURBENERGY": turbenergy, "TURBFREQ": tcross/20, "TURB_KMIN": int(100 * 2*np.pi/L)/100., - "TURB_KMAX": int(100*4*np.pi/(L)+1)/100., "TURB_SIGMA": vrms, "TURB_MINLAMBDA": int(100*R/2)/100, - "TURB_MAXLAMBDA": int(100*R*2)/100, "TURB_COHERENCE_TIME": tcross/2, "UNIT_L": 3.085678e18*length_unit, - "UNIT_M": 1.989e33*mass_unit, "UNIT_V": v_unit*1e2, "UNIT_B": B_unit, "ZINIT": metallicity, "ISRF": ISRF} +paramsfile = str( + open(os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r").read() +) -for k in replacements.keys(): - paramsfile = paramsfile.replace(k, str(replacements[k])) -open("params_"+filename.replace(".hdf5","")+".txt", "w").write(paramsfile) +jet_particle_mass = min(delta_m, max(1e-4, delta_m / 10.0)) +MS_wind_particle_mass = ( + jet_particle_mass / 10 +) # MS winds have lower mdot than jets, so we should be able to better resolve them this way + +replacements = { + "NAME": filename.replace(".hdf5", ""), + "DTSNAP": tff / nsnap, + "MAXTIMESTEP": tff / (nsnap), + "SOFTENING": softening, + "GASSOFT": 2.0e-8, + "TMAX": tff * tmax, + "RHOMAX": ncrit, + "BOXSIZE": boxsize, + "OUTFOLDER": "output", + "JET_PART_MASS": jet_particle_mass, + "MS_WIND_PART_MASS": MS_wind_particle_mass, + "BH_SEED_MASS": delta_m / 2.0, + "TURBDECAY": tcross / 2, + "TURBENERGY": turbenergy, + "TURBFREQ": tcross / 20, + "TURB_KMIN": int(100 * 2 * np.pi / L) / 100.0, + "TURB_KMAX": int(100 * 4 * np.pi / (L) + 1) / 100.0, + "TURB_SIGMA": vrms, + "TURB_MINLAMBDA": int(100 * R / 2) / 100, + "TURB_MAXLAMBDA": int(100 * R * 2) / 100, + "TURB_COHERENCE_TIME": tcross / 2, + "UNIT_L": 3.085678e18 * length_unit, + "UNIT_M": 1.989e33 * mass_unit, + "UNIT_V": v_unit * 1e2, + "UNIT_B": B_unit, + "ZINIT": metallicity, + "ISRF": ISRF, +} + +for k, r in replacements.items(): + paramsfile = paramsfile.replace(k,(r if isinstance(r,str) else '{:.2e}'.format(r))) + +open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write(paramsfile) if makebox: replacements_box = replacements.copy() - replacements_box["NAME"] = filename.replace(".hdf5","_BOX") + replacements_box["NAME"] = filename.replace(".hdf5", "_BOX") replacements_box["BOXSIZE"] = L - replacements_box["TURB_MINLAMBDA"] = int(100*L/2)/100; replacements_box["TURB_MAXLAMBDA"] = int(100*L*2)/100; - paramsfile = str(open(os.path.realpath(__file__).replace("MakeCloud.py","params.txt"), 'r').read()) + replacements_box["TURB_MINLAMBDA"] = int(100 * L / 2) / 100 + replacements_box["TURB_MAXLAMBDA"] = int(100 * L * 2) / 100 + paramsfile = str( + open( + 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) + paramsfile = paramsfile.replace(k, str(replacements_box[k])) + open("params_" + filename.replace(".hdf5", "") + "_BOX.txt", "w").write(paramsfile) if makecylinder: - #Get cylinder params - R_cyl = R * np.sqrt( np.pi/(4*cyl_aspect_ratio) ) #surface density equivalent cylinder - L_cyl = R_cyl*2*cyl_aspect_ratio - vrms_cyl = (2 * G * M_gas / L_cyl)**0.5 * turbulence**0.5 #the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L - vrms_cyl *= 0.71 #additional scaling found numerically to make the stirring run reproduce the right alpha and filament length (similarly determined numerical factor added to GIZMO) - tcross_cyl = 2*R_cyl/vrms_cyl - boxsize_cyl = L_cyl*1.5+R_cyl*5 #the box should fit the cylinder and be many times bigger than its width - print("Cylinder params: L=%g R=%g boxsize=%g vrms=%g"%(L_cyl,R_cyl,boxsize_cyl,vrms_cyl)) + # Get cylinder params + R_cyl = R * np.sqrt( + np.pi / (4 * cyl_aspect_ratio) + ) # surface density equivalent cylinder + L_cyl = R_cyl * 2 * cyl_aspect_ratio + vrms_cyl = ( + 2 * G * M_gas / L_cyl + ) ** 0.5 * turbulence**0.5 # the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L + vrms_cyl *= 0.71 # additional scaling found numerically to make the stirring run reproduce the right alpha and filament length (similarly determined numerical factor added to GIZMO) + tcross_cyl = 2 * R_cyl / vrms_cyl + boxsize_cyl = ( + L_cyl * 1.5 + R_cyl * 5 + ) # the box should fit the cylinder and be many times bigger than its width + print( + "Cylinder params: L=%g R=%g boxsize=%g vrms=%g" + % (L_cyl, R_cyl, boxsize_cyl, vrms_cyl) + ) replacements_cyl = replacements.copy() - replacements_cyl["NAME"] = filename.replace(".hdf5","_CYL") + replacements_cyl["NAME"] = filename.replace(".hdf5", "_CYL") replacements_cyl["BOXSIZE"] = boxsize_cyl - #New driving params - replacements_cyl["TURB_MINLAMBDA"] = int(100*R_cyl)/100; replacements_cyl["TURB_MAXLAMBDA"] = int(100*R_cyl*4)/100; - replacements_cyl["TURB_SIGMA"] = vrms_cyl; replacements_cyl["TURB_COHERENCE_TIME"] = tcross_cyl/2; - #Legacy driving params, probably needs tuning - replacements_cyl["TURBDECAY"] = tcross_cyl/2; replacements_cyl["TURBENERGY"] = 0.019111097819633344*vrms_cyl**3/R_cyl; replacements_cyl["TURBFREQ"] = tcross_cyl/20; - replacements_cyl["TURB_KMIN"] = int(100 * 2*np.pi/R_cyl)/100.; replacements_cyl["TURB_KMAX"] = int(100*4*np.pi/(R_cyl)+1)/100.; - paramsfile = str(open(os.path.realpath(__file__).replace("MakeCloud.py","params.txt"), 'r').read()) + # New driving params + replacements_cyl["TURB_MINLAMBDA"] = int(100 * R_cyl) / 100 + replacements_cyl["TURB_MAXLAMBDA"] = int(100 * R_cyl * 4) / 100 + replacements_cyl["TURB_SIGMA"] = vrms_cyl + replacements_cyl["TURB_COHERENCE_TIME"] = tcross_cyl / 2 + # Legacy driving params, probably needs tuning + replacements_cyl["TURBDECAY"] = tcross_cyl / 2 + replacements_cyl["TURBENERGY"] = 0.019111097819633344 * vrms_cyl**3 / R_cyl + replacements_cyl["TURBFREQ"] = tcross_cyl / 20 + replacements_cyl["TURB_KMIN"] = int(100 * 2 * np.pi / R_cyl) / 100.0 + 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" + ).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) - + paramsfile = paramsfile.replace(k, str(replacements_cyl[k])) + open("params_" + filename.replace(".hdf5", "") + "_CYL.txt", "w").write(paramsfile) + if param_only: - print('Parameters only run, exiting...') + print("Parameters only run, exiting...") exit() -dm = M_gas/N_gas +dm = M_gas / N_gas mgas = np.repeat(dm, N_gas) x = get_glass_coords(N_gas, glass_path) Nx = len(x) -x = 2*(x-0.5) +x = 2 * (x - 0.5) print("Computing radii...") -r = cdist(x, [np.zeros(3)])[:,0] +r = cdist(x, [np.zeros(3)])[:, 0] print("Done! Sorting coordinates...") x = x[r.argsort()][:N_gas] print("Done! Rescaling...") -x *= (float(Nx) / N_gas * 4*np.pi/3 / 8)**(1./3)*R +x *= (float(Nx) / N_gas * 4 * np.pi / 3 / 8) ** (1.0 / 3) * R print("Done! Recomupting radii...") -r = cdist(x, [np.zeros(3)])[:,0] -x, r = x/r.max(), r/r.max() +r = cdist(x, [np.zeros(3)])[:, 0] +x, r = x / r.max(), r / r.max() print("Doing density profile...") -rnew = r**(3. / (3+density_exponent)) * R -x=x * (rnew/r)[:,None] -r = np.sum(x**2, axis=1)**0.5 +rnew = r ** (3.0 / (3 + density_exponent)) * R +x = x * (rnew / r)[:, None] +r = np.sum(x**2, axis=1) ** 0.5 r_order = r.argsort() -x, r = np.take(x,r_order,axis=0), r[r_order] +x, r = np.take(x, r_order, axis=0), r[r_order] if not os.path.exists(turb_path): os.makedirs(turb_path) @@ -322,222 +393,289 @@ def TurbField(res=256, minmode=2, maxmode=64, slope = 2.0, sol_weight=1.0, seed= np.save(fname, vt) else: vt = np.load(fname) - -xgrid = np.linspace(-R,R,vt.shape[-1]) + +xgrid = np.linspace(-R, R, vt.shape[-1]) v = [] for i in range(3): - v.append(interpolate.interpn((xgrid,xgrid,xgrid),vt[i,:,:,:],x)) + v.append(interpolate.interpn((xgrid, xgrid, xgrid), vt[i, :, :, :], x)) v = np.array(v).T print("Coordinates obtained!") - + Mr = mgas.cumsum() -ugrav = G * np.sum(Mr/ r * mgas) -v -= np.average(v,axis=0) -Eturb = 0.5*M_gas/N_gas*np.sum(v**2) -v *= np.sqrt(turbulence*ugrav/Eturb) +ugrav = G * np.sum(Mr / r * mgas) +v -= np.average(v, axis=0) +Eturb = 0.5 * M_gas / N_gas * np.sum(v**2) +v *= np.sqrt(turbulence * ugrav / Eturb) E_rot_target = spin * ugrav -Rcyl = np.sqrt(x[:,0]**2+x[:,1]**2) +Rcyl = np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2) omega = Rcyl**omega_exponent -vrot = np.cross(np.c_[np.zeros_like(omega),np.zeros_like(omega),omega], x) -Erot_actual = np.sum(0.5*mgas[:,None]*vrot**2) -vrot *= np.sqrt(E_rot_target/Erot_actual) +vrot = np.cross(np.c_[np.zeros_like(omega), np.zeros_like(omega), omega], x) +Erot_actual = np.sum(0.5 * mgas[:, None] * vrot**2) +vrot *= np.sqrt(E_rot_target / Erot_actual) v += vrot 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 # alfven speed for unit magnetic field -uB = 0.5*M_gas*vA_unit**2 # magnetic energy we would have for unit magnetic field -if (bfixed>0): +vA_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 +if bfixed > 0: B = B * bfixed else: - B = B * np.sqrt(magnetic_field*ugrav/uB) # renormalize to desired magnetic energy + B = B * np.sqrt( + magnetic_field * ugrav / uB + ) # renormalize to desired magnetic energy v = v - np.average(v, axis=0) x = x - np.average(x, axis=0) -r, phi = np.sum(x**2,axis=1)**0.5, np.arctan2(x[:,1],x[:,0]) -theta = np.arccos(x[:,2]/r) -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)] +r, phi = np.sum(x**2, axis=1) ** 0.5, np.arctan2(x[:, 1], x[:, 0]) +theta = np.arccos(x[:, 2] / r) +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)] +) if makecylinder: - def ind_in_cylinder(x,L_cyl,R_cyl): - return ( (np.abs(x[:,0]) < L_cyl/2) & ( np.sum(x[:,1:]**2,axis=1) 0: - x = np.concatenate([x,x]) - impact_dir = {"x": np.array([1.,0,0]), "y": np.array([0,1,0]), 'z': np.array([0,0,1])}[impact_axis] - impact_param_dir = {"x": np.array([0,1,0]), "y": np.array([0,0,1]), 'z': np.array([1,0,0])}[impact_axis] +u = np.ones_like(mgas) * 0.101 / 2.0 # /2 needed because it is molecular + +if impact_dist > 0: + x = np.concatenate([x, x]) + impact_dir = { + "x": np.array([1.0, 0, 0]), + "y": np.array([0, 1, 0]), + "z": np.array([0, 0, 1]), + }[impact_axis] + impact_param_dir = { + "x": np.array([0, 1, 0]), + "y": np.array([0, 0, 1]), + "z": np.array([1, 0, 0]), + }[impact_axis] x[:N_gas] += impact_dist * R * impact_dir x[N_gas:] -= impact_dist * R * impact_dir x[:N_gas] += 0.5 * impact_param * R * impact_param_dir x[N_gas:] -= 0.5 * impact_param * R * impact_param_dir - v = np.concatenate([v,v]) - vrms = np.sum(v**2,axis=1).mean()**0.5 + v = np.concatenate([v, v]) + vrms = np.sum(v**2, axis=1).mean() ** 0.5 v[:N_gas] -= v_impact * vrms * impact_dir v[N_gas:] += v_impact * vrms * impact_dir - B = np.concatenate([B,B]) - u = np.concatenate([u,u]) - mgas = np.concatenate([mgas,mgas]) + B = np.concatenate([B, B]) + u = np.concatenate([u, u]) + mgas = np.concatenate([mgas, mgas]) -u = np.ones_like(mgas)*(200/v_unit)**2 #start with specific internal energy of (200m/s)^2, this is overwritten unless starting with restart flag 2###### #0.101/2.0 #/2 needed because it is molecular +u = ( + np.ones_like(mgas) * (200 / v_unit) ** 2 +) # start with specific internal energy of (200m/s)^2, this is overwritten unless starting with restart flag 2###### #0.101/2.0 #/2 needed because it is molecular if diffuse_gas: # assuming 10K vs 10^4K gas: factor of ~10^3 density contrast - rho_warm = M_gas*3/(4*np.pi*R**3) / 1000 - M_warm = (boxsize**3 - (4*np.pi*R**3 / 3)) * rho_warm # mass of diffuse box-filling medium - N_warm = int(M_warm/(M_gas/N_gas)) + rho_warm = M_gas * 3 / (4 * np.pi * R**3) / 1000 + M_warm = ( + boxsize**3 - (4 * np.pi * R**3 / 3) + ) * rho_warm # mass of diffuse box-filling medium + N_warm = int(M_warm / (M_gas / N_gas)) if derefinement: x0 = get_glass_coords(N_gas, glass_path) Nx = len(x0) - x0 = 2*(x0-0.5) - r0 = (x0*x0).sum(1)**0.5 + x0 = 2 * (x0 - 0.5) + r0 = (x0 * x0).sum(1) ** 0.5 x0, r0 = x0[r0.argsort()], r0[r0.argsort()] # first lay down the stuff within 3*R - N_warm = int(4*np.pi*rho_warm*(3*R)**3/3/dm) # number of cells within 3R - print(N_warm) - x_warm = x0[:N_warm] * 3 * R / r0[N_warm-1] # uniform density of cells within 3R - x0 = x0[N_warm:] # now we take the ones outside the initial sphere and map them to a n(R) ~ R^-3 profile so that we get constant number of cells per log radius interval + N_warm = int( + 4 * np.pi * rho_warm * (3 * R) ** 3 / 3 / dm + ) # number of cells within 3R + x_warm = ( + x0[:N_warm] * 3 * R / r0[N_warm - 1] + ) # uniform density of cells within 3R + x0 = x0[ + N_warm: + ] # now we take the ones outside the initial sphere and map them to a n(R) ~ R^-3 profile so that we get constant number of cells per log radius interval r0 = r0[N_warm:] - rnew = 3*R*np.exp(np.arange(len(x0))/N_warm/3) - x_warm = np.concatenate([x_warm, (rnew/r0)[:,None] * x0],axis=0) - x_warm = x_warm[np.max(np.abs(x_warm),axis=1) < boxsize/2] + rnew = 3 * R * np.exp(np.arange(len(x0)) / N_warm / 3) + x_warm = np.concatenate([x_warm, (rnew / r0)[:, None] * x0], axis=0) + 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)]) + 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)]) 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] + 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 - B = np.concatenate([B, np.repeat(Bmag,N_warm)[:,np.newaxis] * np.array([0,0,1])]) - u = np.concatenate([u, np.repeat(101.,N_warm)]) + v = np.concatenate([v, np.zeros((N_warm, 3))]) + Bmag = np.average(np.sum(B**2, axis=1)) ** 0.5 + B = np.concatenate( + [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([0, 0, 1])] + ) + u = np.concatenate([u, np.repeat(101.0, N_warm)]) if makecylinder: - #The magnetic field is paralell to the cylinder (true at low densities, so probably fine for IC) - B_cyl = np.concatenate([B, np.repeat(Bmag,N_warm)[:,np.newaxis] * np.array([1,0,0])]) - #Add diffuse medium - M_warm_cyl = (boxsize_cyl**3 - (4*np.pi*R**3 / 3)) * rho_warm - N_warm_cyl = int(M_warm_cyl/(M_gas/N_gas)) - x_warm = boxsize_cyl*np.random.rand(N_warm_cyl, 3) - boxsize_cyl/2 #will be recentered later - x_warm = x_warm[ ~ind_in_cylinder(x_warm,L_cyl,R_cyl) ] #keep only warm gas outside the cylinder - #print("N_warm_cyl: %g N_warm_cyl_kept %g "%(N_warm_cyl,len(x_warm))) + # The magnetic field is paralell to the cylinder (true at low densities, so probably fine for IC) + B_cyl = np.concatenate( + [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([1, 0, 0])] + ) + # Add diffuse medium + M_warm_cyl = (boxsize_cyl**3 - (4 * np.pi * R**3 / 3)) * rho_warm + N_warm_cyl = int(M_warm_cyl / (M_gas / N_gas)) + x_warm = ( + boxsize_cyl * np.random.rand(N_warm_cyl, 3) - boxsize_cyl / 2 + ) # will be recentered later + x_warm = x_warm[ + ~ind_in_cylinder(x_warm, L_cyl, R_cyl) + ] # keep only warm gas outside the cylinder + # print("N_warm_cyl: %g N_warm_cyl_kept %g "%(N_warm_cyl,len(x_warm))) N_warm_cyl = len(x_warm) x_cyl = np.concatenate([x_cyl, x_warm]) - v_cyl = np.concatenate([v_cyl, np.zeros((N_warm,3))]) + v_cyl = np.concatenate([v_cyl, np.zeros((N_warm, 3))]) else: N_warm = 0 -rho = np.repeat(3*M_gas/(4*np.pi*R**3), len(mgas)) -if diffuse_gas: rho[-N_warm:] /= 1000 -h = (32*mgas/rho)**(1./3) +rho = np.repeat(3 * M_gas / (4 * np.pi * R**3), len(mgas)) +if diffuse_gas: + rho[-N_warm:] /= 1000 +h = (32 * mgas / rho) ** (1.0 / 3) -x += boxsize/2 # cloud is always centered at (boxsize/2,boxsize/2,boxsize/2) +x += boxsize / 2 # cloud is always centered at (boxsize/2,boxsize/2,boxsize/2) if makecylinder: - x_cyl += boxsize_cyl/2 + x_cyl += boxsize_cyl / 2 print("Writing snapshot...") -F=h5py.File(filename, 'w') +F = h5py.File(filename, "w") F.create_group("PartType0") F.create_group("Header") -F["Header"].attrs["NumPart_ThisFile"] = [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["NumPart_ThisFile"] = [ + 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) F["PartType0"].create_dataset("Coordinates", data=x) F["PartType0"].create_dataset("Velocities", data=v) -F["PartType0"].create_dataset("ParticleIDs", data=1+np.arange(len(mgas))) +F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(len(mgas))) F["PartType0"].create_dataset("InternalEnergy", data=u) -if M_star>0: +if M_star > 0: F.create_group("PartType5") - #Let's add the sink at the center + # 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("Velocities", data=[v_star]) #at rest - F["PartType5"].create_dataset("ParticleIDs", data=np.array([F["PartType0/ParticleIDs"][:].max() + 1])) - #Advanced properties for sinks - F["PartType5"].create_dataset("BH_Mass", data=M_star) #all the mass in the sink/protostar/star - F["PartType5"].create_dataset("BH_Mass_AlphaDisk", data=np.array([0.])) #starts with no disk - F["PartType5"].create_dataset("BH_Mdot", data=np.array([0.])) #starts with no mdot - F["PartType5"].create_dataset("BH_Specific_AngMom", data=np.array([0.])) #starts with no angular momentum - F["PartType5"].create_dataset("SinkRadius", data=np.array([softening])) #Sinkradius set to softening - F["PartType5"].create_dataset("StellarFormationTime", data=np.array([0.])) - F["PartType5"].create_dataset("ProtoStellarAge", data=np.array([0.])) - F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([5],dtype=np.int32), dtype=np.int32) - #Stellar properties - #if (central_star or central_SN): - #if central_star: -# print("Assuming central sink is a ZAMS star") - #starts as ZAMS star - #else: - #print("Assuming central sink is a ZAMS star about to go supernova") - #F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([6],dtype=np.int32), dtype=np.int32) #starts as ZAMS star going SN - #Set guess for ZAMS stellar radius, will be overwritten - if ((M_star)>1.0): - R_ZAMS = (M_star)**0.57 + 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]) + ) + # Advanced properties for sinks + F["PartType5"].create_dataset( + "BH_Mass", data=M_star + ) # all the mass in the sink/protostar/star + F["PartType5"].create_dataset( + "BH_Mass_AlphaDisk", data=np.array([0.0]) + ) # starts with no disk + F["PartType5"].create_dataset( + "BH_Mdot", data=np.array([0.0]) + ) # starts with no mdot + F["PartType5"].create_dataset( + "BH_Specific_AngMom", data=np.array([0.0]) + ) # starts with no angular momentum + F["PartType5"].create_dataset( + "SinkRadius", data=np.array([softening]) + ) # Sinkradius set to softening + F["PartType5"].create_dataset("StellarFormationTime", data=np.array([0.0])) + F["PartType5"].create_dataset("ProtoStellarAge", data=np.array([0.0])) + F["PartType5"].create_dataset( + "ProtoStellarStage", data=np.array([5], dtype=np.int32), dtype=np.int32 + ) + # Stellar properties + # if (central_star or central_SN): + # if central_star: + # print("Assuming central sink is a ZAMS star") + # starts as ZAMS star + # else: + # print("Assuming central sink is a ZAMS star about to go supernova") + # F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([6],dtype=np.int32), dtype=np.int32) #starts as ZAMS star going SN + # Set guess for ZAMS stellar radius, will be overwritten + if (M_star) > 1.0: + R_ZAMS = (M_star) ** 0.57 else: - R_ZAMS = (M_star)**0.8 - 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.])) #dummy - F["PartType5"].create_dataset("Mass_D", data=np.array([0.])) #No D left - + R_ZAMS = (M_star) ** 0.8 + 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("Mass_D", data=np.array([0.0])) # No D left + if magnetic_field > 0.0: F["PartType0"].create_dataset("MagneticField", data=B) F.close() if makebox: - F=h5py.File(filename.replace(".hdf5","_BOX.hdf5"), 'w') + F = h5py.File(filename.replace(".hdf5", "_BOX.hdf5"), "w") F.create_group("PartType0") F.create_group("Header") - F["Header"].attrs["NumPart_ThisFile"] = [len(mgas),0,0,0,0,0] - F["Header"].attrs["NumPart_Total"] = [len(mgas),0,0,0,0,0] - F["Header"].attrs["MassTable"] = [M_gas/len(mgas),0,0,0,0,0] - F["Header"].attrs["BoxSize"] = (4 * np.pi * R**3 / 3)**(1./3) + F["Header"].attrs["NumPart_ThisFile"] = [len(mgas), 0, 0, 0, 0, 0] + F["Header"].attrs["NumPart_Total"] = [len(mgas), 0, 0, 0, 0, 0] + F["Header"].attrs["MassTable"] = [M_gas / len(mgas), 0, 0, 0, 0, 0] + F["Header"].attrs["BoxSize"] = (4 * np.pi * R**3 / 3) ** (1.0 / 3) 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"]) - F["PartType0"].create_dataset("Velocities", data=np.zeros((len(mgas),3))) - F["PartType0"].create_dataset("ParticleIDs", data=1+np.arange(len(mgas))) + 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"] + ) + F["PartType0"].create_dataset("Velocities", data=np.zeros((len(mgas), 3))) + F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(len(mgas))) F["PartType0"].create_dataset("InternalEnergy", data=u) if magnetic_field > 0.0: - F["PartType0"].create_dataset("MagneticField", data=B[:len(mgas)]) + F["PartType0"].create_dataset("MagneticField", data=B[: len(mgas)]) F.close() - + if makecylinder: - F=h5py.File(filename.replace(".hdf5","_CYL.hdf5"), 'w') + F = h5py.File(filename.replace(".hdf5", "_CYL.hdf5"), "w") F.create_group("PartType0") F.create_group("Header") - F["Header"].attrs["NumPart_ThisFile"] = [N_gas+N_warm_cyl,0,0,0,0,0] - F["Header"].attrs["NumPart_Total"] = [N_gas+N_warm_cyl,0,0,0,0,0] - F["Header"].attrs["MassTable"] = [M_gas/N_gas,0,0,0,0,0] + F["Header"].attrs["NumPart_ThisFile"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] + F["Header"].attrs["NumPart_Total"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] + F["Header"].attrs["MassTable"] = [M_gas / N_gas, 0, 0, 0, 0, 0] F["Header"].attrs["BoxSize"] = boxsize_cyl F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas) - F["PartType0"].create_dataset("Coordinates", data = x_cyl) + 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) F.close() - From bbe5497de5c51d9a83356aeafbb48321ee0449b6 Mon Sep 17 00:00:00 2001 From: Mike Grudic Date: Wed, 25 Sep 2024 11:44:53 -0400 Subject: [PATCH 5/9] formatting workflow --- .github/workflows/python-package.yml | 11 +++++++++++ .github/workflows/python-package.yml~ | 11 +++++++++++ 2 files changed, 22 insertions(+) create mode 100644 .github/workflows/python-package.yml create mode 100644 .github/workflows/python-package.yml~ diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml new file mode 100644 index 0000000..df503e4 --- /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: ". --format" \ No newline at end of file diff --git a/.github/workflows/python-package.yml~ b/.github/workflows/python-package.yml~ new file mode 100644 index 0000000..bd2c714 --- /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: ". --check" \ No newline at end of file From 2d4df3f5c1ecbca9048e13e0eaa94ea0db8cc43a Mon Sep 17 00:00:00 2001 From: Mike Grudic Date: Wed, 25 Sep 2024 11:46:58 -0400 Subject: [PATCH 6/9] format workflow --- .github/workflows/python-package.yml~ | 11 ----------- 1 file changed, 11 deletions(-) delete mode 100644 .github/workflows/python-package.yml~ diff --git a/.github/workflows/python-package.yml~ b/.github/workflows/python-package.yml~ deleted file mode 100644 index bd2c714..0000000 --- a/.github/workflows/python-package.yml~ +++ /dev/null @@ -1,11 +0,0 @@ -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: ". --check" \ No newline at end of file From 467ff5660825baca34450360d83f88f942905acf Mon Sep 17 00:00:00 2001 From: Mike Grudic Date: Wed, 25 Sep 2024 11:50:14 -0400 Subject: [PATCH 7/9] black args fix --- .github/workflows/python-package.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index df503e4..f2f094e 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -8,4 +8,4 @@ jobs: - uses: actions/checkout@v4 - uses: rickstaa/action-black@v1 with: - black_args: ". --format" \ No newline at end of file + black_args: ". -l 119" \ No newline at end of file From d36e8bfbfa2dbfd3df02495bf042d2473659577c Mon Sep 17 00:00:00 2001 From: Mike Grudic Date: Wed, 25 Sep 2024 11:54:06 -0400 Subject: [PATCH 8/9] black args fix --- MakeCloud.py | 95 ++++++++++--- MakeCloud_Bondi.py | 347 +++++++++++++++++++++++++++++---------------- calcGMCProps.py | 186 ++++++++++++------------ 3 files changed, 393 insertions(+), 235 deletions(-) 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 From 538ca8bc1745d71f185ee6f70377759e2b5d3491 Mon Sep 17 00:00:00 2001 From: Mike Grudic Date: Wed, 25 Sep 2024 11:56:57 -0400 Subject: [PATCH 9/9] black args fix --- .github/workflows/python-package.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index f2f094e..3855b91 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -8,4 +8,4 @@ jobs: - uses: actions/checkout@v4 - uses: rickstaa/action-black@v1 with: - black_args: ". -l 119" \ No newline at end of file + black_args: "-l 119 ." \ No newline at end of file