Example Simulation

from simobject import Quantity, Updater, Simulation, DataUpdater

import numpy as np

import astropy.constants as c
import astropy.units as u

Set some constants

au = c.au.cgs.value
rc = 50 * au

Set up the simulation object

sim = Simulation()

# we define the grid, notice that `sim.r` inherits the constant flag

sim.addQuantity('nr', Quantity(100, 'nr of grid points [-]', constant=True))
sim.addQuantity('r0', Quantity(1 * au, 'inner grid radius [cm]', constant=True))
sim.addQuantity('r1', Quantity(1e3 * au, 'outer grid radius [cm]', constant=True))
sim.addQuantity('r', Quantity(np.logspace(np.log10(sim.r0), np.log10(sim.r1), sim.nr), info='radial grid [cm]'))

# time

sim.addQuantity('time', Quantity(0, 'simulation time [s]'))
sim.addQuantity('dt', Quantity(1, 'time step [s]'))

# surface density, to avoid it also being constant, we override that flag

sim.addQuantity('sigma_g', Quantity(200 * (sim.r/rc)**-1 * np.exp(-sim.r/rc), info='gas surface density [g/cm²]', constant=False))

sim.sigma_d = Quantity(sim.sigma_g/100, info='dust surface density [g/cm²]')

Here we define functions how time and densities get updated

def timeupdate(time):
    time += time.owner.dt

def densityupdate(density):
    density *= 0.99

To test those capabilities, we define also print statements to be called in the systole or diastole

def systole_printer(obj):
    print(f'systole of {obj.info}')

def diastole_printer(obj):
    print(f'diastole of {obj.info}')

Now we set/assign those update, systole, diastole functions. We also assign the DataUpdater to the general simulation diastole, so it gets called after all individual diastoles.

sim.time.updater = timeupdate
sim.sigma_g.updater = densityupdate
sim.sigma_d.updater = densityupdate

sim.sigma_d.systoler = systole_printer
sim.sigma_d.diastoler = diastole_printer

sim.diastoler = DataUpdater(['sigma_d', 'sigma_g', 'time'])

Here we define in which order things are updated. We set this to be the same for all quantities. If we didn’t assign anyting, then they would be updated in the order they were added to the object

order = ['sigma_d', 'sigma_g', 'dt', 'time']

sim.systole_order = order
sim.update_order = order
sim.diastole_order = order

Run the update

print(f'time       before update = {sim.time:g}')
print(f'sigma_d[0] before update = {sim.sigma_d[0]:g}')
print(f'sigma_g[0] before update = {sim.sigma_g[0]:g}')

time       before update = 0
sigma_d[0] before update = 98.0199
sigma_g[0] before update = 9801.99

sim.update()

systole of dust surface density [g/cm²]
diastole of dust surface density [g/cm²]
updating data

print(f'time       after update = {sim.time:g}')
print(f'sigma_d[0] after update = {sim.sigma_d[0]:g}')
print(f'sigma_g[0] after update = {sim.sigma_g[0]:g}')

time       after update = 1
sigma_d[0] after update = 97.0397
sigma_g[0] after update = 9703.97

sim.data['time']

array([1])

We also see that sim.data['sigma_d'] now contains an entry as DataUpdater was called as well.

print(f'shape of sim.data[\'sigma_d\']: {sim.data["sigma_d"].shape}')

shape of sim.data['sigma_d']: (1, 100)

If we update once more, the shape(s) increases

sim.update()
sim.update()

systole of dust surface density [g/cm²]
diastole of dust surface density [g/cm²]
updating data
systole of dust surface density [g/cm²]
diastole of dust surface density [g/cm²]
updating data

print(f'shape of sim.data[\'sigma_d\']: {sim.data["sigma_d"].shape}')
print(f'shape of sim.data[\'time\']: {sim.data["time"].shape}')
sim.data['time']

shape of sim.data['sigma_d']: (3, 100)
shape of sim.data['time']: (3, 1)

array([[1],
       [2],
       [3]])