Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 31 additions & 17 deletions MakeCloud.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
#!/usr/bin/env python
"""
MakeCloud: "Believe me, we've got some very turbulent clouds, the best clouds. You're gonna love it."

Usage: MakeCloud.py [options]

Options:
Expand Down Expand Up @@ -30,10 +28,11 @@
--no_diffuse_gas Remove diffuse ISM envelope fills the rest of the box with uniform density.
--phimode=<f> Relative amplitude of m=2 density perturbation (e.g. for Boss-Bodenheimer test) [default: 0.0]
--localdir Changes directory defaults assuming all files are used from local directory.
--B_unit=<gauss> Unit of magnetic field in gauss [default: 1e4]
--B_unit=<gauss> Unit of magnetic field in gauss [default: 1.0]
--length_unit=<pc> Unit of length in pc [default: 1]
--mass_unit=<msun> Unit of mass in M_sun [default: 1]
--v_unit=<m/s> Unit of velocity in m/s [default: 1]
--v_unit=<m/s> Unit of velocity in m/s [default: 1e3]
--unit_system=<name> Units system to adopt (options: starforge_classic (m/s - pc - Msun - T), FIRE (km/s - kpc - 1e10Msun - uG)) [default: None]
--turb_seed=<N> Random seed for turbulence initialization [default: 42]
--tmax=<N> Maximum time to run the simulation to, in units of the freefall time [default: 5]
--nsnap=<N> Number of snapshots per freefall time [default: 150]
Expand Down Expand Up @@ -144,12 +143,27 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42
filename = arguments["--filename"]
diffuse_gas = not arguments["--no_diffuse_gas"]
param_only = arguments["--param_only"]
B_unit = float(arguments["--B_unit"])
length_unit = float(arguments["--length_unit"])
mass_unit = float(arguments["--mass_unit"])
v_unit = float(arguments["--v_unit"])
t_unit = length_unit / v_unit
G = 4300.71 * v_unit**-2 * mass_unit / length_unit
if arguments["--unit_system"] is not None:
match arguments["--unit_system"]:
case "starforge_classic":
length_unit_pc = 1
mass_unit_Msun = 1
v_unit_SI = 1
B_unit_gauss = 1e4
case "FIRE":
length_unit_pc = 1e3
mass_unit_Msun = 1e10
v_unit_SI = 1e3
B_unit_gauss = 1
else:
length_unit_pc = float(arguments["--length_unit"])
mass_unit_Msun = float(arguments["--mass_unit"])
v_unit_SI = float(arguments["--v_unit"])
B_unit_gauss = float(arguments["--B_unit"])

t_unit = length_unit_pc / v_unit_SI

G = 4300.71 * v_unit_SI**-2 * mass_unit_Msun / length_unit_pc
makebox = arguments["--makebox"]
impact_param = float(arguments["--impact_param"])
impact_dist = float(arguments["--impact_dist"])
Expand Down Expand Up @@ -221,7 +235,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42
filename = "".join(filename.split())

dm = M_gas / N_gas
dm_solar = dm / mass_unit
dm_solar = dm / mass_unit_Msun
rho_avg = 3 * M_gas / R**3 / (4 * np.pi)
if dm_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)
Expand Down Expand Up @@ -275,10 +289,10 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42
"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,
"UNIT_L": 3.085678e18 * length_unit_pc,
"UNIT_M": 1.989e33 * mass_unit_Msun,
"UNIT_V": v_unit_SI * 1e2,
"UNIT_B": B_unit_gauss,
"ZINIT": metallicity,
"ISRF": ISRF,
}
Expand Down Expand Up @@ -403,7 +417,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, 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_gauss * (M_gas) ** -0.5 * R**1.5 * np.sqrt(4 * np.pi / 3) / v_unit_SI
) # 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:
Expand Down Expand Up @@ -467,7 +481,7 @@ def ind_in_cylinder(x, L_cyl, R_cyl):
mgas = np.concatenate([mgas, mgas])

u = (
np.ones_like(mgas) * (200 / v_unit) ** 2
np.ones_like(mgas) * (200 / v_unit_SI) ** 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:
Expand Down
65 changes: 4 additions & 61 deletions params.txt
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,10 @@ MinSizeTimestep 1.0e-15 % set this very low, or get the wrong answer


%---- System of units
UnitLength_in_cm UNIT_L % should be 1 pc
UnitMass_in_g UNIT_M % should be 1 solar mass
UnitVelocity_in_cm_per_s UNIT_V % should be 1 m/sec
UnitMagneticField_in_gauss UNIT_B % should be 1 Tesla
UnitLength_in_cm UNIT_L
UnitMass_in_g UNIT_M
UnitVelocity_in_cm_per_s UNIT_V
UnitMagneticField_in_gauss UNIT_B
GravityConstantInternal 0 % calculated by code if =0

%---- Cosmological parameters
Expand Down Expand Up @@ -154,25 +154,6 @@ CritPhysDensity RHOMAX % critical physical density for star formation (c
SfEffPerFreeFall 1.0 % SFR/(Mgas/tfreefall) for gas which meets SF criteria


%---- sub-grid (Springel+Hernquist/GADGET/AREPO) "effective equation of state"
%------- star formation+feedback model (GALSF_EFFECTIVE_EQS on)
MaxSfrTimescale 4.0 % code units (SF timescale at 2-phase threshold)
TempSupernova 3.0e8 % in Kelvin (temp of hot gas in 2-phase model)
TempClouds 1000.0 % in Kelvin (temp of cold gas in 2-phase model)
FactorSN 0.1 % SNe coupling frac (frac of egy retained in hot)
FactorEVP 3000.0 % controls Kennicutt normalization
FactorForSofterEQS 1.0 % interpolate between 'stiff' and isothermal EOS
%------- the sub-grid "decoupled winds" model (GALSF_SUBGRID_WINDS on)
WindEfficiency 2.0 % mass-loading (Mdot_wind = SFR * WindEfficiency)
WindEnergyFraction 0.06 % fraction of SNe energy in winds (sets velocity)
WindFreeTravelMaxTime 0.1 % 'free-stream time' in units of t_Hubble(z)
WindFreeTravelDensFac 0.1 % 'free-stream' until density < this * CritPhysDensity
%------- alternative winds (set GALSF_SUBGRID_WIND_SCALING == 1 or 2)
%------- (scaling with local dark matter properties, as Dave/Oppenheimer/Mannucci/Illustris)
VariableWindVelFactor 1.0 % wind velocity relative to estimated halo v_escape
VariableWindSpecMomentum 5000. % wind momentum per unit stellar mass (code velocity units)



%-------------- FIRE (PFH) explicit star formation & feedback model (FIRE on)
%--- initial metallicity of gas & stars in simulation
Expand Down Expand Up @@ -234,19 +215,6 @@ GrackleDataFile CloudyData_UVB=HM2012.h5
%------------------ Driven Turbulence (Large-Eddy boxes) -----------------
%-------------------------------------------------------------------------

%-------------- Legacy Turbulent stirring parameters (TURB_DRIVING on)
%ST_decay TURBDECAY % decay time for driving-mode phase correlations
%ST_energy TURBENERGY % energy of driving-scale modes: sets norm of turb
%ST_DtFreq TURBFREQ % time interval for driving updates (set by hand)
%ST_Kmin TURB_KMIN % minimum driving-k: should be >=2.*M_PI/All.BoxSize
%ST_Kmax TURB_KMAX % maximum driving-k: set to couple times Kmin or more if cascade desired
%ST_SolWeight 1.0 % fractional wt of solenoidal modes (wt*curl + (1-wt)*div)
%ST_AmplFac 1.0 % multiplies turb amplitudes
%ST_SpectForm 3 % driving pwr-spec: 0=Ek~const; 1=sharp-peak at kc; 2=Ek~k^(-5/3); 3=Ek~k^-2
%ST_Seed 42 % random number seed for modes (so you can reproduce it)
%IsoSoundSpeed 200. % initializes gas sound speed in box to this value
%TimeBetTurbSpectrum DTSNAP % time (code units) between evaluations of turb pwrspec

%-------------- Turbulent stirring parameters (TURB_DRIVING on)
TurbDrive_ApproxRMSVturb TURB_SIGMA
TurbDrive_MinWavelength TURB_MINLAMBDA
Expand All @@ -267,31 +235,6 @@ TurbDrive_TimeBetTurbSpectrum DTSNAP
DarkEnergyConstantW -1 % time-independent DE parameter w, used only if no table
TabulatedCosmologyFile CosmoTbl % table with cosmological parameters


%-------------------------------------------------------------
%------------------ Solid bodies and Impacts -----------------
%-------------------------------------------------------------

%-------------- Parameters for custom Tillotson equation-of-state (EOS_TILLOTSON on)
%--- In ICs, set "CompositionType": 0=custom,1=granite,2=basalt,3=iron,4=ice,5=olivine/dunite,6=water;
%--- their EOS parameters will be set accordingly. If CompositionType=0, the custom parameters below
%--- are used, matched to the definitions in Table A1 of Reinhardt+Stadel 2017,MNRAS,467,4252 (below is iron)
Tillotson_EOS_params_a 0.5 % a parameter [dimensionless]
Tillotson_EOS_params_b 1.5 % b parameter [dimensionless]
Tillotson_EOS_params_u_0 9.5e10 % u_0 parameter in [erg/g]
Tillotson_EOS_params_rho_0 7.86 % rho_0 parameter in [g/cm^3]
Tillotson_EOS_params_A 1.28e12 % A parameter in [erg/cm^3]
Tillotson_EOS_params_B 1.05e12 % B parameter in [erg/cm^3]
Tillotson_EOS_params_u_s 1.42e10 % u_s parameter in [erg/g]
Tillotson_EOS_params_u_s_prime 8.45e10 % u_s^prime parameter in [erg/g]
Tillotson_EOS_params_alpha 5.0 % alpha parameter [dimensionless]
Tillotson_EOS_params_beta 5.0 % beta parameter [dimensionless]
Tillotson_EOS_params_mu 7.75e11 % elastic shear modulus in [erg/cm^3] (used if EOS_ELASTIC is on)
Tillotson_EOS_params_Y0 8.5e10 % hugoniot elastic limit in [erg/cm^3] (used if EOS_ELASTIC is on)




%--- Developer-Mode Parameters (usually hard-coded, but set manually if DEVELOPER_MODE is on) --------
ErrTolTheta 0.5 % 0.7=standard
ErrTolForceAcc 0.0025 % 0.0025=standard
Expand Down