Running a simulation with dust growth and fragmentation

The dust growth and fragmentation algorithm is described in Vericel et al. (in prep).

Important

Dust growth is not restricted to single star simulations! You can use dustgrowth around binaries, you can put planets in your disc, you can tilt and warp everything if you want and many more. Be creative :).

If you find a bug, please report it either on the phantom slack channel or by sending me an email at arnaud.vericel@univ-lyon1.fr

The algorithm is tested automatically by the nightly test suite, but you can check the test manually in your phantom repository:

cd ~/phantom; make testgrowth

Dust growth is a compile time option and you have two ways of setting it up.

You can import the Makefile using the dustydisc setup and set DUSTGROWTH to yes before compiling:

~/phantom/scripts/writemake.sh dustydisc > Makefile
export DUSTGROWTH=yes
make; make setup

You can directly import the Makefile using the “growingdisc” setup:

~/phantom/scripts/writemake.sh growingdisc > Makefile
make; make setup

When you make the setup file ./phantomsetup ilovegrowth, dust growth adds an option block after the dust section:

# options for growth and fragmentation of dust
               ifrag =           1    ! fragmentation of dust (0=off,1=on,2=Kobayashi)
               isnow =           0    ! snow line (0=off,1=position based,2=temperature based)
               rsnow =        100.    ! snow line position in AU
               Tsnow =        150.    ! snow line condensation temperature in K
               vfrag =         15.    ! uniform fragmentation threshold in m/s
             vfragin =       5.000    ! internal fragmentation threshold in m/s
            vfragout =         15.    ! external fragmentation threshold in m/s
        grainsizemin =       0.005    ! minimum allowed grain size in cm

Here’s a brief description of each of them (remember that they are described in more details in Vericel et al. (in prep))

ifrag =           1    ! fragmentation of dust (0=off,1=on,2=Kobayashi)

switch choosing between pure growth (= 0), growth and harsh fragmentation (= 1) or growth and smoother fragmentation (= 2)

isnow =           0    ! snow line (0=off,1=position based,2=temperature based)

switch setting up a snow line by distance to the star (= 1) or by its corresponding temperature (= 2)

rsnow =        100.    ! snow line position in AU

snow line distance to the star if isnow = 1

Tsnow =        150.    ! snow line condensation temperature in K

snow line corresponding sublimation temperature if isnow = 2

vfrag =         15.    ! uniform fragmentation threshold in m/s

uniform fragmentation threshold velocity above which fragmentation occurs (if ifrag > 0 and isnow = 0)

 vfragin =       5.000    ! inward fragmentation threshold in m/s
vfragout =         15.    ! inward fragmentation threshold in m/s

internal and external fragmentation threshold velocities above which fragmentation occurs (if ifrag > 0 and isnow > 0)

grainsizemin =       0.005    ! minimum allowed grain size in cm

minimum size you allow dust grains to reach if ifrag > 0

After re-executing the setup file ./phantomsetup ilovegrowth.setup, ilovegrowth.in contains an extra dust growth block only with the relevant informations, e.g here:

# options controlling growth
               ifrag =           1    ! dust fragmentation (0=off,1=on,2=Kobayashi)
        grainsizemin =       0.005    ! minimum grain size in cm
               isnow =           0    ! snow line (0=off,1=position based,2=temperature based)
               vfrag =         15.    ! uniform fragmentation threshold in m/s

The relative velocity between dust particles in computed using the viscosity parameter α which is stored in the shearparam parameter in the physical viscosity block of the infile.

# options controlling physical viscosity
           irealvisc =           0    ! physical viscosity type (0=none,1=const,2=Shakura/Sunyaev)
          shearparam =       0.010    ! magnitude of shear viscosity (irealvisc=1) or alpha_SS (irealvisc=2)

Independently of the use of physical viscosity or not (irealvisc = 0 or 2), the algorithm will use shearparam as the value for α. However, if you set irealvisc to 1 and set a constant viscosity, dustgrowth will be unavailable. In the case where you have no physical viscosity, make sure shearparam is equal to the alphaSS that you specified during the setup. This is important for consistency.

After launching the simulation ./phantom ilovegrowth.in, the full dumpfiles will contain the grain sizes, grain densities, the ratio Vrel/Vfrag (if ifrag > 0, else gives Vrel/(1 m/s)), the Stokes number St and other usefull informations. Meanwhile, small dumps will only contain the grainsizes and graindensities.

Tips to perform great dust growth simulations

Set more gas particles than dust particles. I typically recommend a ratio of 10 gas particles per dust particle (e.g. 10^6 gas and 10^5 dust).

Set the dust disc slighlty inside the gas disc to avoid numerical artefacts at the outer edge of the disc (e.g, rout_gas = 300 and rout_dust = 250).

I also suggest to eventually use another moddump to remove the few particles that got ejected at large distances from the star. This will accelerate your simulation. You can make that utility with make moddump MODFILE=moddump_removeparticles_radius.f90.

If you consider fragmentation, grainsizemin should not be too small or else small grains will dictate the timestepping and make your simulation ridiculously slow. I typically recommend to set that minimum between 10 and 50 μm.

Make synthetic images out of dustgrowth dumps with mcfost

Dustgrowth dumps are not naturally ready to be interpreted by mcfost. However, a Python module is available for transforming a dump into a friendlier format for mcfost.

To do so, you can open your favorite Python code editor and import the module:

# you need to copy the script to your current directory before importing it
import os
os.system(f"cp path_to_phantom/scripts/growthtomcfost.py .")
import growthtomcfost as gtmcf

And simply call the function pimp_my_sim, which takes a few arguments. You can learn more about those by printing the docstring attached to the function by typing gtmc.pimp_my_sim? in iPython or in a Jupyter Notebook. This should give you the following:

Signature:
g.pimp_my_sim(
    gdump_name,
    outdump_name,
    path_to_phantom,
    bins_per_dex=5,
    force_smax=False,
    smax_user=0.2,
    maxdustlarge=25,
    compil_logfile='compil.log',
    logfile_moddump='moddump.log',
    logfile_phantom='phantom.log',
    save_plot=False,
    scale='log',
    color='black',
    show=True)
Docstring:
pimp_my_sim transforms a dustgrowth dump into a multi large grains dump.
it runs phantom on that dump for 1 time step to recompute the densities.
the output dump is ready to be processed by mcfost through command line or pymcfost.

input parameters are:

gdump_name      : (str)   - name of dustgrowth dump (input)
outdump_name    : (str)   - prefix of the name of desired output dump (output) - final name will be outdump_name_00000
path_to_phantom : (str)   - path to phantom's directory
bins_per_dex    : (int)   - number of bins per magnitude of size
force_smax      : (bool)  - wether or not to force a maximum size for binning, else find it automatically
smax_user       : (float) - size of forced maximum in cm
maxdustlarge    : (int)   - maximum number of large dust species for memory allocation in phantom
compil_logfile  : (str)   - name of logfile to store compilation output
logfile_moddump : (str)   - name of logfile to store moddump run output
logfile_phantom : (str)   - name of logfile to store phantom run output
save_plot       : (bool)  - wether or not to save plot of size distribution in pdf file, else shows it interactively
scale           : (str)   - npart axis scale, options are "linear" or "log"
color           : (str)   - histogram color
show            : (bool)  - wether or not to show the distribution, only applies if save_plot=True
File:      /Users/Arnaud/Documents/Codes/phantom/scripts/growthtomcfost.py
Type:      function

The resulting dump is ready to be processed by mcfost by typing

mcfost <paramfile> -phantom <output_dump> -<options>

The advantage of doing this in a Jupyter Notebook or a python file is that the transition is made easy with pymcfost which can be called directly after.

Have fun :)