3. Advanced Customization

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

from tripodpy import Simulation

sim = Simulation()

Customizing the Grids

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

In that case, we 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
    Nr           : NoneType
    OmegaK       : NoneType
    r            : NoneType
    ri           : NoneType
  -----

sim.makegrids()

sim.grid

Group (Grid quantities)
-----------------------
    A            : Field (Radial grid annulus area [cm²]), 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: The Keplerian frequency has not been initialized at this point.

The Radial Grid

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

radius_ratio = sim.grid.r[1:]/sim.grid.r[:-1]
print('maximum ratio=', max(radius_ratio),'; minimum ratio=', min(radius_ratio))

maximum ratio= 1.0715193052376093 ; minimum ratio= 1.0715193052376046

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, we 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. Note that this and many othr useful functions for setting up your simulations ad adding additional physics can be found in the Dustpylib repository.

import numpy as np

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

    Parameters
    ----------
    ri : array
        Radial interfaces
    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
    for _ in range(num):
        ind = np.argmin(r0 > ri) - 1
        ind_left  = ind-num
        ind_right = ind+num+1
        ri_left  = ri[:ind_left]
        ri_right = ri[ind_right:]
        number_refined_cells = (2*num+1)*2
        ri_mid = np.empty(number_refined_cells)
        for i in range(0, number_refined_cells, 2):
            j = ind-num+int(i/2)
            ri_mid[i] = ri[j]
            ri_mid[i+1] = 0.5*(ri[j]+ri[j+1])
        ri = np.concatenate((ri_left, ri_mid, ri_right))
    return ri

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_unrefined = np.logspace(0., 3., num=100, base=10.) * c.au

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

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize=(8,1))
ax.plot(ri/c.au, ri*0, 'rx', label='local refined')
ax.plot(ri_unrefined/c.au, ri_unrefined*0, 'b+', label='even log-spaced')
ax.set_xscale('log')
ax.set_xlabel('cell faces')
ax.legend(ncol=2)
ax.set_xlim(3,8)
ax.set_yticks([])
plt.show()

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
    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
    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 we 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

120

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

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  True  True  True  True  True  True]

Customizing the Physics of a Field

In this example, we want to have a fragmentation velocity that depends on the temperature in the disk. If the temperature is 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 that particles coated in water ice are stickier than pure silicate particles and can consequently withstand higher collision velocities (see e.g. Pinilla et al. (2017)). Keep in mind, however, 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()

now lets look at the fields we want to modify, sim.dust.v.frag. Each field has an updater that is called to update the values of the field in each timestep:

sim.dust.v.frag.updater

Heartbeat
---------

Systole:  None
Updater:  None
Diastole: None

Docstrings
----------

Systole:
The type of the None singleton.

Updater:
The type of the None singleton.

Diastole:
The type of the None singleton.

We can see that the updater is empty (None) for the details of updaters and fields work beyond what we explain in the documentation i refere the reader to Simframe the underlying simulation framework we use.

The fragmentation velocity is of shape (Nr,), i.e. there is one value at every location in the grid.

sim.dust.v.frag.shape

(120,)

At the moment it is 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. To this end, 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

(120,)

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.Heartbeat object. When calling the update() function, the heartbeat which consists of a systole, the actual updater, and a diastole will be executed. 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, we can directly assign the function only to the updater leaving the systoles/diastoles as they are.

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 and needs to be instructed to update itself. We can either update the whole simulation frame with Simulation.update(), or we merely 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: Upon customizing a quantity on which other quantities depend, one also has to update these quantities. In our case, this would be the sticking/fragmentation probabilites. Therefore, it is always better to update the entire 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]")

0.0

The grid group now has a new member.

sim.grid

Group (Grid quantities)
-----------------------
    A            : Field (Radial grid annulus area [cm²]), 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}\). As a consequence, the disk temperature decreases. 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 TwoPopPy instructions how to update the snowline location, but we have not yet told it to actually update it regularily.

TwoPopPy 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', 'components', 'gas', 'dust']

This means that if we call Simulation.update(), we 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.grid.updateorder

['OmegaK']

sim.gas.updateorder

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

sim.gas.S.updateorder

['ext', 'tot']

sim.gas.v.updateorder

['visc', 'rad']

sim.dust.updateorder

['delta',
 'rhos',
 'fill',
 'backreaction',
 'f',
 'qrec',
 'a',
 'm',
 'St',
 'H',
 'rho',
 'D',
 'eps',
 'v',
 'p',
 'q',
 'SigmaFloor',
 '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', 'fragtrans', 'driftfrag']

sim.dust.S.updateorder

['compo', 'ext', 'tot', 'smax_hyd']

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, 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 we can see, the grid updater is only updating the Keplerian frequency, but not our snowline location. So we can simply add it to the list.

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

If we 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 we previously saw from the update order the grid is updated before the gas. The snowline location, however, needs the gas temperature and, therefore, has to be updated after the gas. On the other hand, 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. These are functions that are called before and after the actual updater respectively. 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. Alternatively, 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()

We then 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 has been 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}\).

sim.t.snapshots = np.logspace(3., 4., num=21, base=10.) * 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()

tripodpy v1.0.0

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
Execution time: 0:00:02

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

from tripodpy import plot

sim.dust.v.rel.tot[:,-2,-1]

[9.89573147e+01 9.99700269e+01 9.99816380e+01 9.99795476e+01
 9.99793393e+01 9.99784526e+01 9.99755832e+01 9.99729235e+01
 1.00062050e+02 1.02801513e+02 1.66129999e+02 9.94209033e+02
 9.99064410e+02 9.98508517e+02 9.97316291e+02 9.95912211e+02
 9.94586883e+02 9.94388105e+02 9.93159217e+02 9.91990969e+02
 9.91810650e+02 9.92189259e+02 9.91472052e+02 9.90759690e+02
 9.91007207e+02 9.90579368e+02 9.90026175e+02 9.89886066e+02
 9.89765992e+02 9.89620537e+02 9.89471329e+02 9.89316184e+02
 9.89153082e+02 9.88979790e+02 9.88768627e+02 9.88588472e+02
 9.88396302e+02 9.88148220e+02 9.87500185e+02 9.87598115e+02
 9.87589396e+02 9.87167335e+02 9.85846687e+02 9.86131466e+02
 9.86006149e+02 9.83087789e+02 9.83956664e+02 9.83020641e+02
 9.80871611e+02 9.78475302e+02 9.75475351e+02 9.71412769e+02
 9.65484095e+02 9.56281980e+02 9.41445651e+02 9.17181625e+02
 8.77656420e+02 8.14612987e+02 7.18125446e+02 5.82398472e+02
 4.26173835e+02 2.89905167e+02 1.89920778e+02 1.21401378e+02
 7.00025016e+01 3.06967233e+01 8.22527489e+00 6.17477792e+00
 4.80836892e+00 3.87983125e+00 3.43761907e+00 2.90657122e+00
 2.55028645e+00 2.30944635e+00 2.14906082e+00 2.04775325e+00
 1.99226007e+00 1.97442200e+00 1.98948036e+00 2.03510348e+00
 2.11084532e+00 2.21787855e+00 2.35891932e+00 2.53830480e+00
 2.76221229e+00 3.03902974e+00 3.37990768e+00 3.79954576e+00
 4.31730106e+00 4.95876021e+00 5.75801536e+00 6.76106858e+00
 8.03113020e+00 9.65655090e+00 1.12968943e+01 1.39588282e+01
 1.75675389e+01 2.25791417e+01 2.90941472e+01 3.92961547e+01
 5.52673998e+01 8.35954126e+01 1.33788306e+02 2.10997169e+02
 3.35328346e+02 5.48781762e+02 9.32635389e+02 1.64510219e+03
 2.97594296e+03 5.28541001e+03 8.00115238e+03 8.00982059e+03
 9.02852515e+03 8.59801065e+03 4.90500343e+03 2.27115662e+03
 9.62479116e+02 3.81478456e+02 1.31810488e+02 1.32071602e+01]

plot.panel(sim,show_limits=True)

The change in fragmentation velocity has an obvious effect on the particle 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()