3. Advanced Customization

The core principle of DustPy is that you can change anything easily. Not only the initial conditions as shown in the previous chapter, but also the physics behind the simulation.

from dustpy import Simulation

sim = Simulation()

Customizing the Grids

By default the radial and the mass grid will be created when calling Simulation.initialize(). But there can be situations where you need to know the grid sizes before completely initializing the Simulation object. For example if you want to create custom fields and you need to initialize them with the correct shape.

In that case you can call Simulation.makegrids() to only create the grids without initializing the simulation objects. In fact, Simulation.makegrids() is by default called within Simulation.initialize().

sim.grid

Group (Grid quantities)
-----------------------
    A            : NoneType
    m            : NoneType
    Nm           : NoneType
    Nr           : NoneType
    OmegaK       : NoneType
    r            : NoneType
    ri           : NoneType
  -----

sim.makegrids()

sim.grid

Group (Grid quantities)
-----------------------
    A            : Field (Radial grid annulus area [cm²]), constant
    m            : Field (Mass grid [g]), constant
    Nm           : Field (# of mass bins), constant
    Nr           : Field (# of radial grid cells), constant
    r            : Field (Radial grid cell centers [cm]), constant
    ri           : Field (Radial grid cell interfaces [cm]), constant
  -----
    OmegaK       : NoneType
  -----

Note that the Keplerian frequency has not been initialized at this point.

The Radial Grid

By default the radial grid is a regular logarithmic grid. Meaning, the ratio of adjacent grid cells is constant.

sim.grid.r[1:]/sim.grid.r[:-1]

[1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931 1.07151931
 1.07151931 1.07151931 1.07151931]

As explained in the previous chapter, the location of the grid boundaries and the number of grid cells can be controlled via Simulation.ini.grid.rmin, Simulation.ini.grid.rmax, and Simulation.ini.grid.Nr.

Simulation.makegrids() will use these parameters to create the radial grid.

But it is also possible to completely customize the radial grid. To do so you have to set the locations of the radial grid cell interfaces Simulation.grid.ri before calling either Simulation.makegrids() or Simulation.initialize().

In this example we simply want to refine the grid at a given location. We use this helper function, which takes an existing grid ri and doubles the number of grid cells in a region num grid cells on both sides around location r0. We also recursively call this function with reduced num to even further refine the grid and to have a smooth transition between the high and low resolution regions.

import numpy as np

def refinegrid(ri, r0, num=3):
    """Function to refine the radial grid

    Parameters
    ----------
    ri : array
        Radial grid
    r0 : float
        Radial location around which grid should be refined
    num : int, option, default : 3
        Number of refinement iterations

    Returns
    -------
    ri : array
        New refined radial grid"""
    if num == 0:
        return ri
    ind = np.argmin(r0 > ri) - 1
    indl = ind-num
    indr = ind+num+1
    ril = ri[:indl]
    rir = ri[indr:]
    N = (2*num+1)*2
    rim = np.empty(N)
    for i in range(0, N, 2):
        j = ind-num+int(i/2)
        rim[i] = ri[j]
        rim[i+1] = 0.5*(ri[j]+ri[j+1])
    ri = np.concatenate((ril, rim, rir))
    return refinegrid(ri, r0, num=num-1)

We now create a regular logarithmic grid and feed it to our function. We want to refine the grid in a location around \(4.5\,\mathrm{AU}\).

import dustpy.constants as c

ri = np.logspace(0., 3., num=100, base=10.) * c.au

ri = refinegrid(ri, 4.5*c.au, num=3)

We can now create a new empty Simulation object, assign the grid cell interfaces and initialize the grids.

sim = Simulation()

sim.grid.ri = ri

sim.grid

Group (Grid quantities)
-----------------------
    A            : NoneType
    m            : NoneType
    Nm           : NoneType
    Nr           : NoneType
    OmegaK       : NoneType
    r            : NoneType
    ri           : ndarray
  -----

Note: it is sufficient to assign a numpy.ndarray to Simulation.grid.ri and not a simframe.Field.

We can now make the grids.

sim.makegrids()

sim.grid

Group (Grid quantities)
-----------------------
    A            : Field (Radial grid annulus area [cm²]), constant
    m            : Field (Mass grid [g]), constant
    Nm           : Field (# of mass bins), constant
    Nr           : Field (# of radial grid cells), constant
    r            : Field (Radial grid cell centers [cm]), constant
    ri           : Field (Radial grid cell interfaces [cm]), constant
  -----
    OmegaK       : NoneType
  -----

As you can see, Simulation.grid.ri was automatically converted to a simframe.Field and the other fields were created. The number of radial grid cells is greater than \(100\) as we added more grid cells.

sim.grid.Nr

114

To see that we actually refined the grid at the correct location, we can plot the location of the radial grid cells.

import matplotlib.pyplot as plt

fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogy(sim.grid.r/c.au)
ax.axhline(4.5, c="gray", lw=1)
ax.set_xlabel("# of radial grid cell")
ax.set_ylabel("Location of radial grid cell [AU]")
fig.tight_layout()
plt.show()

The position of the radial grid cells have to be exactly in the center between their grid cell interfaces and are automatically calculated by Simulation.makegrids().

sim.grid.r == 0.5 * (sim.grid.ri[1:] + sim.grid.ri[:-1])

[ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True]

The Mass Grid

You should NEVER set the mass grid manually! The mass grid has to be strictly logarithmic. Only customize the mass grid by setting Simulation.ini.grid.mmin, Simulation.ini.grid.mmax, and Simulation.ini.grid.Nmbpd.

If you have to create your own non-logarithmic mass grid for some reason, be aware that you have to re-write the entire coagulation algorithm as well, since it only conserves mass on a logarithmic grid.

Customizing the Physics of a Field

In this example we want to have a fragmentation velocity that depends on the temperature in the disk. Is the temperature below 150 K, we want to have a fragmentation velocity of 10 m/s, otherwise it shall be 1 m/s. The idea behind this approach is than particles coated in water ice are stickier that pure silicate particles and can widthstand higher collision velocities. See for example Pinilla et al. (2017). However, keep in mind that newer experiments suggest that particles covered in water ice do not have a beneficial collisional behavior, see e.g. Musiolik & Wurm (2019).

First, we initialize our simulation object.

sim.initialize()

The fragmentation velocity has the shape (Nr,), meaning there is one value at every location in the grid.

sim.dust.v.frag.shape

(114,)

But right now it’s constant.

fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogx(sim.grid.r/c.au, sim.dust.v.frag)
ax.set_xlabel("Distance from star [AU]")
ax.set_ylabel("Fragmentation velocity [cm/s]")
fig.tight_layout()
plt.show()

We have to write a function that takes the simulation object as input parameter and returns our desired fragmentation velocities. We can use the fact that the gas temperature has the same shape. Keep in mind that everything has to be in cgs units.

sim.gas.T.shape

(114,)

def v_frag(sim):
    return np.where(sim.gas.T<150., 1000., 100)

We can now assign this function to the updater of the dust fragmentation velocities. For details of this process, please have a look at the Simframe documentation.

sim.dust.v.frag.updater = v_frag

The updater of a group/field stores a simframe.Heatbeat object. When calling the update() function the heartbeat will be executed which consists of a systole, the actual updater, and a diastole. The systole is executed before the actual update functions, the diastole afterwards.

When assigning a function (or None) to the updater of a group/field a new Heartbeat object will be created with empty systoles and diastoles only executing the update function. If the existing updater already has systoles/diastoles, those would be overwritten with an empty function.

To prevent this you can directly assign the function only to the updater leaving the systoles/diastoles as they are.

sim.dust.v.updater.updater = v_frag

The systoles/diastoles can be set with the following command. Only for demonstration here, since we assign None. Read more about this in the section about Systoles and Diastoles.

sim.dust.v.updater.systole = None
sim.dust.v.updater.diastole = None

As of now, the simulation object still holds the old data for the fragmentation velocity. We have to tell it to update itself. We can either update the whole simulation frame with Simulation.update(), or we just update the fragmentation velocities.

sim.dust.v.frag.update()

The fragmentation velocities should now show our desired behavior.

fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogx(sim.grid.r/c.au, sim.dust.v.frag)
ax.set_xlabel("Distance from star [AU]")
ax.set_ylabel("Fragmentation velocity [cm/s]")
fig.tight_layout()
plt.show()

Note: If you customized a quantity on which other quantities depend on, you also have to update these quantities. In our case this would be the sticking/fragmentation probabilites. So it is always better to update the whole simulation frame.

sim.update()

Adding Custom Fields

We can not only modify existing fields, we can also create our own fields.

In this example we want to add another field rsnow to Simulation.grid, that gives us the location of the so called snowline, i.e., the location in the disk where water ice starts to sublime.

First, we add the field and initialize it with zero.

sim.grid.addfield("rsnow", 0., description="Snowline location [cm]")

The grid group has now a new member.

sim.grid

Group (Grid quantities)
-----------------------
    A            : Field (Radial grid annulus area [cm²]), constant
    m            : Field (Mass grid [g]), constant
    Nm           : Field (# of mass bins), constant
    Nr           : Field (# of radial grid cells), constant
    OmegaK       : Field (Keplerian frequency [1/s])
    r            : Field (Radial grid cell centers [cm]), constant
    ri           : Field (Radial grid cell interfaces [cm]), constant
    rsnow        : Field (Snowline location [cm])
  -----

As a next step we have to write a function that returns us the location of the snowline. Here we simply use the first grid cell where the temperature is smaller than \(150\,\mathrm{K}\) and return the value of the inner interface of that grid cell.

def rsnow(sim):
    isnow = np.argmax(sim.gas.T<150.)
    return sim.grid.ri[isnow]

We assign this function to the updater of our snowline field.

sim.grid.rsnow.updater.updater = rsnow

And update the field.

sim.grid.rsnow.update()

print("The snowline is located at {:4.2f} AU.".format(sim.grid.rsnow/c.au))

The snowline is located at 2.15 AU.

Right now the temperature is constant throughout the simulation, because the stellar parameters do not change. To see an effect in our snowline location, we need to have a changing temperature profile.

To achieve this, we let the stellar radius decrease from a value of \(3\,M_\odot\) to \(2\,M_\odot\) within the first \(10,000\,\mathrm{yrs}\). This results in decreasing disk temperature. This is only for demonstration purposes and is not necessarily physical.

def Rstar(sim):
    dR = -1.*c.R_sun
    dt = 1.e4 * c.year
    m = dR/dt
    R = m*sim.t + 3.*c.R_sun
    R = np.maximum(R, c.R_sun)
    return R

And assign this to the updater of the stellar radius.

sim.star.R.updater.updater = Rstar

Modifying the Update Order

But we are still not done, yet. We have given DustPy instructions how to update the snowline location, but we have not yet told it to actually update it regularily.

DustPy calls Simulation.update(), the updater of the simulation object, once per timestep after the integration step and just before writing the data files. The updater of a group/field is basically a list of groups/fields, whose updater is called in that order.

For the main simulation object this is

sim.updateorder

['star', 'grid', 'gas', 'dust']

This means that if you call Simulation.update() you basically call Simulation.star.update(), Simulation.grid.update(), Simulation.gas.update(), and Simulation.dust.update() in that order.

The updaters of the sub-groups and fields look as follows

sim.star.updateorder

['M', 'R', 'T', 'L']

sim.grid.updateorder

['OmegaK']

sim.gas.updateorder

['gamma',
 'mu',
 'T',
 'alpha',
 'cs',
 'Hp',
 'nu',
 'rho',
 'n',
 'mfp',
 'P',
 'eta',
 'S']

sim.gas.S.updateorder

['ext', 'tot']

sim.gas.v.updateorder

['visc', 'rad']

sim.dust.updateorder

['delta',
 'rhos',
 'fill',
 'a',
 'St',
 'H',
 'rho',
 'backreaction',
 'v',
 'D',
 'eps',
 'kernel',
 'p',
 'S']

sim.dust.backreaction.updateorder

['A', 'B']

sim.dust.delta.updateorder

['rad', 'turb', 'vert']

sim.dust.Fi.updateorder

['adv', 'diff', 'tot']

sim.dust.p.updateorder

['frag', 'stick']

sim.dust.S.updateorder

['ext', 'tot']

sim.dust.v.updateorder

['frag', 'driftmax', 'rel']

sim.dust.v.rel.updateorder

['azi', 'brown', 'rad', 'turb', 'vert', 'tot']

Note: The gas updater does not contain the updaters of Simulation.gas.v and Simulation.gas.Fi and the updater of the gas sources does not contain the updater of Simulation.gas.S.hyd. These are quantities that are calculated in the finalization step of the integrator, since they are derived from the result of the implicit gas integration.

The same is true for Simulation.dust.v.rad, Simulation.dust.Fi, Simulation.dust.S.coag, and Simulation.dust.S.hyd in the dust updater, which are also calculated from the implicit integration in the finalization step of the integrator

As you can see, the grid updater is only updating the Keplerian frequency, but not our snowline location. So we can simply adding it to the list.

sim.grid.updater = ["OmegaK", "rsnow"]

If you assign lists to updaters, their systoles and diastoles will always be overwritten with None.

sim.grid.updateorder

['OmegaK', 'rsnow']

Systoles and Diastoles

However, the previous solution has a conceptional problem. As you can see from the update order previously the grid is updated before the gas. The snowline location, however, needs the gas temperature and, therefore, has to be updated after the gas. But we also cannot update the grid as a whole after the gas, because the gas updaters need the Keplerian frequency. We need another solution.

But first, we revert the grid updater.

sim.grid.updater = ["OmegaK"]

Every updater has a systole and a diastole. That is a function that is called before respectively after the actual updater. Since no other quantity depends on our snowline location, we can simply update it at the end and put it in the diastole of the main updater. Or we could assign it to the diastole of the gas temperature updater, since it only requires the updated gas temperature.

We therefore write a diastole function, that is updating the snowline location separately.

def diastole(sim):
    sim.grid.rsnow.update()

And assign this function to the diastole of the gas temperature updater.

sim.gas.T.updater.diastole = diastole

Now every time Simulation.gas.T.update() is called, Simulation.grid.rsnow.update() will be called at the end of it.

Customizing the Snapshots

As already explained in a previous chapter, the snapshots can be customized by simply setting Simulation.t.snapshots. In this example we only want to run the simulation for \(10,000\,\mathrm{yrs}\). The snapshots have to include the starting time to write out the initial conditions.

sim.t.snapshots = np.hstack(
    [sim.t, np.geomspace(1.e3, 1.e4, num=21) * c.year]
)

We can now change the data directory to avoid an overwrite error and start the simulation with our modifications.

sim.writer.datadir = "3_data"

sim.run()

DustPy v1.0.6

Documentation: https://stammler.github.io/dustpy/
PyPI:          https://pypi.org/project/dustpy/
GitHub:        https://github.com/stammler/dustpy/

Please cite Stammler & Birnstiel (2022).

Checking for mass conservation...

    - Sticking:
        max. rel. error:  2.75e-14
        for particle collision
            m[114] =  1.93e+04 g    with
            m[116] =  3.73e+04 g
    - Full fragmentation:
        max. rel. error:  5.55e-16
        for particle collision
            m[55] =  7.20e-05 g    with
            m[55] =  7.20e-05 g
    - Erosion:
        max. rel. error:  1.78e-15
        for particle collision
            m[110] =  5.18e+03 g    with
            m[118] =  7.20e+04 g

Creating data directory 3_data.
Writing file 3_data/data0000.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0001.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0002.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0003.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0004.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0005.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0006.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0007.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0008.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0009.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0010.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0011.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0012.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0013.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0014.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0015.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0016.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0017.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0018.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0019.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0020.hdf5
Writing dump file 3_data/frame.dmp
Writing file 3_data/data0021.hdf5
Writing dump file 3_data/frame.dmp
Execution time: 0:05:30

We can now have a look at the result of our modifications.

from dustpy import plot

plot.panel(sim)

The change in fragmentation velocity has an obvious effect on the particles sizes.

To check the time evolution of the snowline, we have to read the data. The gray lines are the positions of the radial grid cell interfaces and snapshots, which explains the discrete behavior of the snowline location.

t = sim.writer.read.sequence("t") / c.year
ri = sim.writer.read.sequence("grid.ri") / c.au
rsnow = sim.writer.read.sequence("grid.rsnow") / c.au

fig = plt.figure(dpi=150)
ax = fig.add_subplot(111)
ax.semilogx(t, rsnow)
ax.hlines(ri, 1.e3, 1e4, lw=1, color="gray", alpha=0.25)
ax.vlines(t, 2.5, 5.5, lw=1, color="gray", alpha=0.25)
ax.set_xlim(1.e3, 1.e4)
ax.set_ylim(2.5, 5.5)
ax.set_xlabel("Time [yr]")
ax.set_ylabel("Snowline location [AU]")
fig.tight_layout()
plt.show()