<jupyter_start><jupyter_text>Accretion discIn this tutorial we set up a protoplanetary disc around a star represented by a sink particle, and we add a planet. This notebook generates a Phantom "temporary" dump file that can be read by Phantom as an initial condition. It also generates a Phantom "in" file. Together, these files can start a Phantom simulation.  InitializationFirst we import the required modules.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import phantomsetup
from phantomsetup import constants
from astropy import units as u<jupyter_output><empty_output><jupyter_text>Here we set some constants for convenience.<jupyter_code>igas = phantomsetup.defaults.PARTICLE_TYPE['igas']
put_planet=False<jupyter_output><empty_output><jupyter_text>ParametersNow we set the parameters for the problem.First is the `prefix` which sets the file name for the dump file and Phantom in file.<jupyter_code>prefix = 'disc'<jupyter_output><empty_output><jupyter_text>ResolutionWe choose the resolution to be $10^6$ gas particles.<jupyter_code>number_of_particles = 10000
particle_type = igas<jupyter_output><empty_output><jupyter_text>UnitsWe set the length and mass units to be au and solar masses, respectively. We will also set the time unit such that the gravitational constant is unity.<jupyter_code>length_unit = phantomsetup.units.unit_string_to_cgs('au')
mass_unit = phantomsetup.units.unit_string_to_cgs('solarm')
gravitational_constant = 1.0

#find unit time

time_unit = np.sqrt(length_unit ** 3 / (constants.G * mass_unit))*u.second
time_unit=time_unit.to_value(u.yr)

print("Time Unit:",time_unit,"yr")<jupyter_output>Time Unit: 0.15916067096618036 yr<jupyter_text>Runtime options<jupyter_code># end time (years)
tmax=5000
# time between dumps (years)
dtmax=5
#full dump every n dumps
nfulldump=10
#total number of dumps
ndumps=-1
#tree accuracy
tree_accuracy=0.2

#convert to code units
tmax=tmax/time_unit
dtmax=dtmax/time_unit<jupyter_output><empty_output><jupyter_text>ViscosityThe SPH $\alpha$ viscosity parameter is its minimal value of 0.1.<jupyter_code>alpha_artificial = 0.99
alphamax=1<jupyter_output><empty_output><jupyter_text>StarThe star is of solar mass, at the origin, with a 5 au accretion radius.<jupyter_code>stellar_mass = 1.0
stellar_accretion_radius = 1.0
stellar_position = (0.0, 0.0, 0.0)
stellar_velocity = (0.0, 0.0, 0.0)<jupyter_output><empty_output><jupyter_text>DiscThe disc has mass 0.01 solar masses, it extends from 10 au to 200 au.<jupyter_code>radius_min = 1.0
radius_max = 100.0
disc_mass = 0.01<jupyter_output><empty_output><jupyter_text>Equation of stateThe equation of state is locally isothermal. We set the aspect ratio H/R at a reference radius.<jupyter_code>ieos = 6
isink=1

q_index = 0.5  #0.75
aspect_ratio = 0.5 #0.05
reference_radius = 10.0<jupyter_output><empty_output><jupyter_text>PlanetWe add a planet at 100 au.<jupyter_code>planet_mass = 0.001
planet_position = (5.2, 0.0, 0.0)
orbital_radius = np.linalg.norm(planet_position)
planet_velocity = (0,np.sqrt(gravitational_constant * stellar_mass / orbital_radius),0)<jupyter_output><empty_output><jupyter_text>We set the planet accretion radius as a fraction of the Hill sphere radius.<jupyter_code>planet_accretion_radius_fraction_hill_radius = 0.25
planet_hill_radius = phantomsetup.orbits.hill_sphere_radius(
    orbital_radius, planet_mass, stellar_mass
)
planet_accretion_radius = (
    planet_accretion_radius_fraction_hill_radius * planet_hill_radius
)<jupyter_output><empty_output><jupyter_text>Surface density distributionFor the surface density distribution we use the Lynden-Bell and Pringle (1974) self-similar solution, i.e. a power law with an exponential taper.<jupyter_code>def density_distribution(radius, radius_critical, gamma):
    """Self-similar disc surface density distribution.

    This is the Lyden-Bell and Pringle (1974) solution, i.e. a power law
    with an exponential taper.
    """
    return phantomsetup.disc.self_similar_accretion_disc(radius, radius_critical, gamma)

radius_critical = 100.0
gamma = 1.0
args = (radius_critical, gamma)<jupyter_output><empty_output><jupyter_text>Instantiate the `Setup` objectThe following instantiates the `phantomsetup.Setup` object.<jupyter_code>setup = phantomsetup.Setup()<jupyter_output><empty_output><jupyter_text>Set attributes and add particles PrefixSet the prefix.<jupyter_code>setup.prefix = prefix<jupyter_output><empty_output><jupyter_text>UnitsSet units.<jupyter_code>setup.set_units(
    length=length_unit, mass=mass_unit, gravitational_constant_is_unity=True
)<jupyter_output><empty_output><jupyter_text>Equation of stateSet the equation of state. We get `polyk` from the aspect ratio parametrization.<jupyter_code>polyk = phantomsetup.eos.polyk_for_locally_isothermal_disc(
    q_index, reference_radius, aspect_ratio, stellar_mass, gravitational_constant
)
setup.set_run_option('isink', isink)

setup.set_equation_of_state(ieos=ieos, polyk=polyk)<jupyter_output><empty_output><jupyter_text>ViscositySet the numerical viscosity to Phantom disc viscosity.<jupyter_code>setup.set_dissipation(disc_viscosity=True, alpha=alpha_artificial)<jupyter_output><empty_output><jupyter_text>GravityInclude disc self-gravity and define tree accuracy<jupyter_code>setup.set_compile_option('GRAVITY', True)
setup.set_run_option('tree_accuracy', tree_accuracy)
if ndumps<0 and dtmax>0:
    setup.set_output(tmax=float(tmax),dtmax=float(dtmax), nfulldump=nfulldump)
elif dtmax<0 and ndumps>0:
    setup.set_output(tmax=float(tmax),ndumps=ndumps, nfulldump=nfulldump)
else:
    sys.exit("Issue with ndumps, dtmax (only one of the two have to be set)")<jupyter_output><empty_output><jupyter_text>Options for injecting/removing particles<jupyter_code>setup.set_compile_option('INJECT_PARTICLES',False)
setup.set_run_option('rkill',-1) #off<jupyter_output><empty_output><jupyter_text>StarAdd a star at the origin.<jupyter_code>setup.add_sink(
    mass=stellar_mass,
    accretion_radius=stellar_accretion_radius,
    position=stellar_position,
    velocity=stellar_velocity,
)<jupyter_output><empty_output><jupyter_text>DiscAdd the disc around the star.<jupyter_code>disc = phantomsetup.Disc(
    particle_type=particle_type,
    number_of_particles=number_of_particles,
    disc_mass=disc_mass,
    density_distribution=density_distribution,
    radius_range=(radius_min, radius_max),
    q_index=q_index,
    aspect_ratio=aspect_ratio,
    reference_radius=reference_radius,
    stellar_mass=stellar_mass,
    gravitational_constant=gravitational_constant,
    extra_args=(radius_critical, gamma),
)

setup.add_container(disc)<jupyter_output><empty_output><jupyter_text>PlanetAdd a planet in orbit around the star.<jupyter_code>if (put_planet==True):
    setup.add_sink(
        mass=planet_mass,
        accretion_radius=planet_accretion_radius,
        position=planet_position,
        velocity=planet_velocity,
    )<jupyter_output><empty_output><jupyter_text>PlotNow we plot some quantities to see what we have set up.First is the particles in the xy-plane. The sink particles are marked in red.<jupyter_code>x, y, z = disc.arrays['position'][:, 0], disc.arrays['position'][:, 1], disc.arrays['position'][:, 2]

fig, ax = plt.subplots()
ax.plot(x[::10], y[::10], 'k.', ms=0.5)
for sink in setup.sinks:
    ax.plot(sink.position[0], sink.position[1], 'ro')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')<jupyter_output><empty_output><jupyter_text>Next we plot the particles in the rz-plane.<jupyter_code>position_cylindrical, velocity_cylindrical = phantomsetup.geometry.coordinate_transform(
    position=disc.arrays['position'],
    velocity=disc.arrays['velocity'],
    geometry_from='cartesian',
    geometry_to='cylindrical'
)

R = position_cylindrical[:, 0]

fig, ax = plt.subplots()
ax.plot(R[::10], z[::10], 'k.', ms=0.5)
ax.set_xlabel('$R$')
ax.set_ylabel('$z$')
ax.set_aspect('equal')
ax.set_ylim(bottom=2 * z.min(), top=2 * z.max())<jupyter_output><empty_output><jupyter_text>Finally, we plot $v_{\phi}$ as a function of radius.<jupyter_code>vphi = velocity_cylindrical[:, 1]

fig, ax = plt.subplots()
ax.plot(R[::10], vphi[::10], 'k.', ms=0.5)
ax.set_xlabel('$R$')
ax.set_ylabel('$v_{\phi}$')<jupyter_output><empty_output><jupyter_text>Write to fileNow that we are happy with the setup, write the "temporary" dump file with the initial conditions and the Phantom "in" file.First we set a working directory for the simulation.<jupyter_code>working_dir = '~/phantom-runs/disc'
setup.write_dump_file(directory=working_dir)
setup.write_in_file(directory=working_dir)<jupyter_output><empty_output><jupyter_text>Viscosity options (Need to be here to ensure compilation happens with the right flag)<jupyter_code>#all needs to be False to start with
setup.set_compile_option('DISC_VISCOSITY', False)

#set the required flag here
viscosity_flag='DISC_VISCOSITY'

if viscosity_flag=='DISC_VISCOSITY':
    setup.set_compile_option('DISC_VISCOSITY', True)
elif viscosity_flag == 'CONST_AV':
    setup.set_compile_option('CONST_AV', True)
elif viscosity_flag == 'MORRIS_MONAGHAN':
    setup.set_compile_option("MORRIS_MONAGHAN", True)<jupyter_output><empty_output><jupyter_text>Compile PhantomYou can start a Phantom calculation from these two files but you must compile Phantom with the correct Makefile variables. We can use the `phantom_compile_command` method to show how Phantom would be compiled.<jupyter_code>setup.set_compile_option("STORE_TEMPERATURE",True)
#note: in Makefile the parameter is called "STORE_TEMP"
print(setup.phantom_compile_command(hdf5root='/opt/local/'))<jupyter_output>make \
  SETUP=empty \
  SYSTEM=gfortran \
  HDF5=yes \
  HDF5ROOT=/opt/local/ \
  DISC_VISCOSITY=yes \
  DRIVING=no \
  DUST=no \
  DUSTGROWTH=no \
  GRAVITY=yes \
  H2CHEM=no \
  IND_TIMESTEPS=yes \
  INJECT_PARTICLES=no \
  ISOTHERMAL=yes \
  KERNEL=cubic \
  MAXDUSTSMALL=11 \
  MAXDUSTLARGE=11 \
  MCFOST=no \
  MHD=no \
  NONIDEALMHD=no \
  PERIODIC=no \
  PHOTO=no \
  MORRIS_MONAGHAN=no \
  CONST_AV=no \
  STORE_TEMPERATURE=yes<jupyter_text>To compile Phantom we can use the `compile_phantom` method to compile Phantom. For example, the following will compile Phantom located in `phantom_dir` and copy the binary to `working_dir`.<jupyter_code>result = setup.compile_phantom(
    phantom_dir='~/phantom/',
    working_dir=working_dir,
    hdf5root='/opt/local'
)<jupyter_output>Compilation successful.
{'fraction_non_alphanumeric': 0.1026648773006135, 'fraction_numerical': 0.016008435582822087, 'mean_word_length': 6.75111441307578, 'pattern_counts': {'":': 0, '<': 120, '<?xml version=': 0, '>': 120, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '127260', 'license': 'MIT', 'n_tokens_mistral': 3886, 'n_tokens_neox': 3623, 'n_words': 793, 'path': 'examples/make_disc.ipynb', 'repo_name': 'distamio/phantom-setup', 'size': '10432'}