-
Notifications
You must be signed in to change notification settings - Fork 3
Home
Make sure you have the following packages installed:
- Numpy
- Scipy
- Pynbody
- matplotlib
-
custom_python_packages
Contains additional utilities used by ICgen. To install, clone the repository at https://github.com/ibackus/custom_python_packages.git and add to your PYTHONPATH environmental variable -
changa_uw
Currently, the UW version of ChaNGa is required.
You can clone the repo using https://github.com/ibackus/ICgen.git
Make sure to add the created ICgen directory to your PYTHONPATH environmental variable
ChaNGa is required by ICgen. To allow ICgen to properly execute ChaNGa, you'll have to edit or create a preset for running ChaNGa. This can be done in python and is pretty simple. This only needs to be done once, the settings will be saved.
ChaNGa presets are store as a 4 element list (or tuple). These basically supply the command-line arguments which are used for executing ChaNGa and they are all strings [runner, runner_arguments, ChaNGa, ChaNGa_arguments]
- runner e.g., 'charmrun' or 'mpirun', etc...
- runner_arguments Additional arguments used by the runner. For no extra arguments, set to an empty string ''
- ChaNGa Command for running ChaNGa, e.g. 'ChaNGa' or '~/bin/ChaNGa_uw'
- ChaNGa_arguments Additional arguments supplied to ChaNGa. For no extra arguments, set to an empty string ''
Here's an example of setting up the 'local' preset.
import ICgen
# Echo the default global settings if you want
ICgen.global_settings()
# Edit the 'local' changa preset
ICgen.global_settings['changa_presets']['local'] = ['charmrun', '++local', 'ChaNGa', '-D 3']
# Change which preset is used by default
ICgen.global_settings['changa_presets']['default'] = 'local'
# Save settings
ICgen.global_settings.save()
If needed, the default settings can be restored:
ICgen.global_settings.restore_defaults()
ICgen.global_settings.save()
The following is an example script for generating a set of initial conditions which can be found here: IC_example.py
import ICgen
import pynbody
SimArray = pynbody.array.SimArray
Initialize a blank initial conditions (IC) object:
IC = ICgen.IC()
Echo the default settings to get an idea of what parameters can be changed
IC.settings()
Set up the surface density profile settings. Notice that right now the only setting under 'Sigma profile parameters' is kind: None
. Lets use a simple powerlaw with cut-offs at the interior and edge of the disk. The available kinds are 'powerlaw', 'MQWS', and 'exponential'
IC.settings.sigma.kind = 'powerlaw'
Now echo the sigma settings (defaults)
IC.settings.sigma()
The important parameters to set are:
- Rd : Disk radius, should be a pynbody SimArray
- Qmin : Minimum estimated Toomre Q in the disk. This will in general differ from the actual Q min.
- n_points : number of radial points to calculate sigma at
The other parameters include:
- rin : Dimensionless radius (multiple of disk radius) where the interior cut-off is applied.
- rmax : Maximum dimensionless radius (multiple of disk radius) to evaluate surface density at.
- cutlength : Dimensionless length scale of the exterior exponential cutoff (multiple of disk radius).
Let's generate a disk with powerlaw up to 10 au followed by a cutoff with a Q of approximately 1.0. The actual Toomre Q will differ from this number significantly. This approximation is just a simple way of setting the surface density normalization. To save time, lets do it at low resolution
IC.settings.sigma.Qmin = 1.0
IC.settings.sigma.n_points = 100
IC.settings.sigma.Rd = SimArray(10.0, 'au')
Lets be careful and save what we've done. This will save the ICs to IC.p in the current directory
IC.save()
Let's change the settings used for numerically calculating the gas density. For this we'll use a low resolution of 100 to save time. The default of 1000 is more robust
IC.settings.rho_calc() # Echo defaults
IC.settings.rho_calc.nr = 100 # Number of radial points to calculate on
IC.settings.rho_calc.nz = 100 # Number of vertical points to calculate on
Set the number of particles
IC.settings.pos_gen.nParticles = 50000
Set up the temperature profile to use. Available kinds are 'powerlaw' and 'MQWS' We'll use something of the form T = T0(r/r0)^Tpower
IC.settings.physical.kind = 'powerlaw'
IC.settings.physical.Tpower = -0.5 # exponent
IC.settings.physical.T0 = SimArray(150.0, 'K') # temperature at r0
IC.settings.physical.Tmin = SimArray(10.0, 'K') # Minimum temperature
IC.settings.physical.r0 = SimArray(1.0, 'au')
Let's set the star mass and gas mass assuming H2
IC.settings.physical.M = SimArray(1.0, 'Msol') # star mass in solar masses
IC.settings.physical.m = SimArray(2.0, 'm_p') # mass of H2
Have changa run on the local preset
IC.settings.changa_run.preset = 'local'
Save our work to IC.p
IC.save()
IC.settings()
All that needs to be done now is to tell the ICs to generate. This can be done step-by-step or all at once.
To do it all at once:
IC.generate()
IC.save()
To do it step-by-step:
IC.maker.sigma_gen() # Generate surface density profile and CDF
IC.maker.rho_gen() # Calculate density according to hydrodynamic equilibrium
IC.maker.pos_gen() # Generate particle positions
IC.maker.snapshot_gen() # Generate the final tipsy snapshot with velocities etc
IC.save()
This will run through the whole procedure and save a tipsy snapshot to snapshot.std with a basic .param file saved to snapshot.param