Custom Initial Conditions¶
This notebook will teach you how to set up a custom gas distribution as initial conditions.
First we create a new directory and change to it.
example_name = "400_custom_initial_conditions"
example_dir = f"example_dirs/{example_name}"
import os
if not os.path.basename(os.getcwd()) == example_name:
!mkdir -p $example_dir
os.chdir(example_dir)
repo_root = os.path.abspath(os.path.join(os.getcwd(), "../../../"))
print(f"Current working directory: {os.getcwd()}")
print(f"Repository root directory: {repo_root}")
Current working directory: /home/rometsch/repo/fargocpt/examples/example_dirs/400_custom_initial_conditions
Repository root directory: /home/rometsch/repo/fargocpt
Make the code¶
Make sure the code is built by running make again.
If you have not yet compiled the go, please go to the readme and follow the instructions there. You can also try to run the following cell directly, but it will only output error messages. This might make debugging harder.
%%timeit -n1 -r1
from sys import platform
if platform in ["linux", "darwin"]:
!make -j 4 -C $repo_root/src > make.log
else:
raise RuntimeError(f"Seems like you are not running MacOS or Linux but {platform}. This is unsupported. You are on your own, good luck!")
119 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
Preparing a setup file¶
We’ll take the example setup file from the examples directory and modify it in python. If you want to create setup files for a parameter study, just copy the code and make your own setup creator script.
configfile = "setup.yml"
!cp $repo_root/examples/config.yml $configfile
We’ll use the ruamel.yaml package to read and write the setup file. This can be set up to preserve comments which is very useful if you want to trace your decisions later on.
try:
import ruamel.yaml
except ImportError:
raise ImportError("Please install ruamel.yaml with `python3 -m pip install ruamel.yaml`")
yamlloader = ruamel.yaml.YAML()
with open(configfile, "r") as infile:
config = yamlloader.load(infile)
config["MonitorTimestep"] = 0.314 # monitor scalar files around every half orbit
config["Nmonitor"] = 20 # write a snapshot every orbit
config["Nsnapshots"] = 10 # wirte 100 snapshots
# use very low resolution by setting it to 2 cell per scaleheight, cps
config["cps"] = 2
with open(configfile, "w") as outfile:
yamlloader.dump(config, outfile)
Getting a data skeleton¶
The next step is to run the code for zero timesteps. This creates the zeroth and reference snapshots, which we can then modify to our liking.
from fargocpt import run
run(["start", configfile, "-N", "0"], np=2, nt=1, exe=repo_root+"/bin/fargocpt_exe", detach=False)
Running command: mpirun -np 2 --report-pid /tmp/tmp72h2aq3h -x OMP_NUM_THREADS=1 /home/rometsch/repo/fargocpt/bin/fargocpt_exe start setup.yml -N 0
fargo process pid 1410359
[0] MPI rank # 0 runs as process 1410363
[0] MPI rank # 0 OpenMP thread # 0 of 1 on cpt-kamino
[1] MPI rank # 1 runs as process 1410364
[1] MPI rank # 1 OpenMP thread # 0 of 1 on cpt-kamino
[0] fargo: This file was compiled on Nov 14 2023, 12:56:40.
[0] fargo: This version of FARGO used _GNU_SOURCE
[0] fargo: This version of FARGO used NDEBUG. So no assertion checks!
[0] Using parameter file setup.yml
[0] Computing disk quantities within 5.00000e+00 L0 from coordinate center
[0] BC: Inner composite = reflecting
[0] BC: Outer composite = reflecting
[0] BC: Sigma inner = zerogradient
[0] BC: Sigma outer = zerogradient
[0] BC: Energy inner = zerogradient
[0] BC: Energy outer = zerogradient
[0] BC: Vrad inner = reflecting
[0] BC: Vrad outer = reflecting
[0] BC: Vaz inner = keplerian
[0] BC: Vaz outer = keplerian
[0] DampingTimeFactor: 1.00000e-01 Outer damping time is computed at radius of 2.50000e+00
[0] Damping VRadial to reference value at inner boundary.
[0] Damping VRadial to reference value at outer boundary.
[0] Damping VAzimuthal to reference value at inner boundary.
[0] Damping VAzimuthal to reference value at outer boundary.
[0] Damping SurfaceDensity to reference value at inner boundary.
[0] Damping SurfaceDensity to reference value at outer boundary.
[0] Damping Energy to reference value at inner boundary.
[0] Damping Energy to reference value at outer boundary.
[0] Radiative diffusion is disabled. Using fixed omega = 1.500000 with a maximum 50000 interations.
[0] Indirect Term computed as effective Hydro center acceleratrion with shifting the Nbody system to the center.
[0] Body force on gas computed via potential.
[0] Using FARGO algorithm for azimuthal advection.
[0] Using standard forward euler scheme for source terms.
[0] Cps is set, overwriting Nrad and Naz!
[0] Grid resolution set using cps = 2.000000
[0] The grid has (Nrad, Naz) = (74, 251) cells with (1.994113, 1.997395) cps.
[0] Computing scale height with respect to primary object.
[0] Using isothermal equation of state. AdiabaticIndex = 1.400.
[0] Viscosity is of alpha type with alpha = 1.000e-03
[0] Defaulting to VanLeer flux limiter
[0] Output information:
[0] Output directory: output/out/
[0] Number of files: 80
[0] Total output size: 0.00 GB
[0] Space Available: 31.24 GB
[0] Initializing 1 RNGs per MPI process.
[0] Warning : no `radii.dat' file found. Using default.
[0] The first 1 planets are used to calculate the hydro frame center.
[0] The mass of the planets used as hydro frame center is 1.000000e+00.
[0] 2 planet(s) initialized.
[0] Planet overview:
[0]
[0] # | name | mass [m0] | x [l0] | y [l0] | vx | vy |
[0] -----+-------------------------+------------+------------+------------+------------+------------+
[0] 0 | Star | 1 | 0 | -0 | 0 | 0 |
[0] 1 | Jupiter | 0.0009546033 | 1 | 0 | -0 | 1.000477 |
[0]
[0] # | e | a | T [t0] | T [a] | accreting | Accretion Type |
[0] -----+------------+------------+------------+------------+------------+----------------+
[0] 0 | 6.368246e-17 | 1 | 6.280188 | 0.999548 | 0 | No Accretion |
[0] 1 | 6.368246e-17 | 1 | 6.280188 | 0.999548 | 0 | No Accretion |
[0]
[0] # | Temp [K] | R [l0] | irradiates | rampuptime |
[0] -----+------------+------------+------------+------------+
[0] 0 | 5778 | 0.0046505 | yes | 0 |
[0] 1 | 0 | 4.6505e-05 | no | 0 |
[0]
[0] Using Tscharnuter-Winkler (1979) artificial viscosity with C = 1.410000.
[0] Artificial viscosity is used for dissipation.
[0] Surface density factor: 2.50663
[0] Tau factor: 0.5
[0] Tau min: 0.01
[0] Kappa factor: 1
[0] Minimum temperature: 2.81162e-05 K = 3.00000e+00
[0] Maximum temperature: 9.37206e+94 K = 1.00000e+100
[0] Heating from viscous dissipation is enabled. Using a total factor of 1.
[0] Cooling (beta) is disabled and reference temperature is floor. Using beta = 10.
[0] Cooling (radiative) is enabled. Using a total factor of 1.
[0] S-curve cooling is disabled.
[0] CFL parameter: 0.5
[0] Opacity uses tables from Lin & Papaloizou, 1985
[0] Particles are disabled.
[0] Initializing Sigma(r) = 2.25093e-05 = 200 g cm^-2 * [r/(1 AU)]^(-0.5)
[0] Total disk is mass is 0.000348841 = 6.9366e+29 g.
[0] Writing output output/out/snapshots/0, Snapshot Number 0, Time 0.000000.
[0] Writing output output/out/snapshots/reference, Snapshot Number 0, Time 0.000000.
[0] -- Final: Total Hydrosteps 0, Time 0.00, Walltime 0.01 seconds, Time per Step: 0.00 milliseconds
0
# %matplotlib widget
from fargocpt import Overview
Overview("output/out",
vars=["2:Sigma",
"1:Sigma",
"1:vazi:rel",
"1:vrad"]).create()
Defining helper functions and classes¶
import numpy as np
from dataclasses import dataclass
from types import SimpleNamespace
import astropy.units as u
import astropy.constants as const
@dataclass
class FargoCPTField:
outputdir: str
snapshotid: str
name: str
def __post_init__(self):
self.filename = f"{self.outputdir}/snapshots/{self.snapshotid}/{self.name}.dat"
self.grid = get_fargo_grid(self.outputdir)
if self.name == "vrad":
self.R, self.Phi = np.meshgrid(self.grid.ri, self.grid.phic, indexing="ij")
elif self.name == "vtheta":
self.R, self.Phi = np.meshgrid(self.grid.rc, self.grid.phii[:-1], indexing="ij")
else:
self.R, self.Phi = np.meshgrid(self.grid.rc, self.grid.phic, indexing="ij")
self._load()
def _load(self):
self._data = np.fromfile(self.filename, dtype=np.float64)
self._data = self._data.reshape(self.R.shape[0], self.R.shape[1])
def save(self, altid=None):
if altid is not None:
filename = f"{self.outputdir}/snapshots/{altid}/{self.name}.dat"
self._data.tofile(self.filename)
@property
def array(self):
return self._data
@array.setter
def array(self, data):
if not self._data.shape == data.shape:
raise ValueError("Shape of data does not match shape of field")
self._data = data
def get_fargo_grid(outputdir):
Nrad, Naz = np.genfromtxt(f"{outputdir}/dimensions.dat", usecols=(4,5), dtype=int, unpack=True)
ri = np.genfromtxt(f"{outputdir}/used_rad.dat")
phii = np.linspace(0, 2*np.pi, Naz+1)
Ri, Phii = np.meshgrid(ri, phii, indexing="ij")
Xi = Ri*np.cos(Phii)
Yi = Ri*np.sin(Phii)
rc = 2/3*(ri[1:]**2/(ri[1:]+ri[:-1]) + ri[:-1]) # approx center in polar coords
phic = 0.5*(phii[1:]+phii[:-1])
Rc, Phic = np.meshgrid(rc, phic, indexing="ij")
Xc = Rc*np.cos(Phic)
Yc = Rc*np.sin(Phic)
dphi = phii[1] - phii[0]
dr = ri[1:] - ri[:-1]
A = 0.5*(Ri[1:,1:]**2 - Ri[:-1,1:]**2)*dphi
return SimpleNamespace(
Nrad=Nrad, Naz=Naz,
ri=ri, phii=phii, Ri=Ri, Phii=Phii, Xi=Xi, Yi=Yi,
rc=rc, phic=phic, Rc=Rc, Phic=Phic, Xc=Xc, Yc=Yc,
dphi=dphi, dr=dr, A=A
)
def get_fargo_powerlaw_disk(config: dict):
l0 = u.Unit(config["l0"])
m0 = u.Unit(config["m0"])
t0 = u.Unit(config["t0"])
Temp0 = u.Unit(config["Temp0"])
Mstar = u.Quantity(config["nbody"][0]["mass"]).to(m0)
GM = (const.G * Mstar).to(l0**3/t0**2).value
Rgas = (const.k_B / const.m_p).to(l0**2/t0**2/Temp0).value
mu = config["mu"]
r0 = 1
Sigma0 = u.Quantity(config["Sigma0"])
if not Sigma0.unit.is_unity():
Sigma0 = Sigma0.to_value(m0/l0**2)
betaSigma = -float(config["SigmaSlope"])
if config["EquationOfState"][:3] == "iso":
gamma = 1
else:
gamma = float(config["AdiabaticIndex"])
h0 = float(config["AspectRatio"])
betah = float(config["FlaringIndex"])
return PowerlawDisk(GM, Rgas, mu, r0, Sigma0, betaSigma, gamma, h0=h0, betah=betah)
def get_fargo_units(outdir: str) -> dict:
with open(f"{outdir}/units.yml", "r") as f:
units_info = yaml.safe_load(f)
units = {}
for key, val in units_info.items():
units[key] = u.Unit(f"{val['cgs value']} {val['cgs symbol']}")
return units
def get_fargo_powerlaw_disk_output(outdir: str):
config = get_fargo_config(f"{outdir}/snapshots/0/config.yml")
units = get_fargo_units(outdir)
l0 = units["length"]
m0 = units["mass"]
t0 = units["time"]
Temp0 = units["temperature"]
Mstar = u.Quantity(config["nbody"][0]["mass"]).to(m0)
GM = (const.G * Mstar).to(l0**3/t0**2).value
Rgas = (const.k_B / const.m_p).to(l0**2/t0**2/Temp0).value
mu = config["mu"]
r0 = 1
Sigma0 = u.Quantity(config["Sigma0"])
if not Sigma0.unit.is_unity():
Sigma0 = Sigma0.to_value(m0/l0**2)
betaSigma = -float(config["SigmaSlope"])
if config["EquationOfState"][:3] == "iso":
gamma = 1
else:
gamma = float(config["AdiabaticIndex"])
h0 = float(config["AspectRatio"])
betah = float(config["FlaringIndex"])
return PowerlawDisk(GM, Rgas, mu, r0, Sigma0, betaSigma, gamma, h0=h0, betah=betah)
import yaml
def get_fargo_config(filename):
with open(filename, "r") as f:
config = yaml.safe_load(f)
return config
@dataclass
class PowerlawDisk:
"""
Powerlaw disk model
Scaleheight is defined with the isothermal sound speed, so H = c_s/sqrt(gamma) / Omega_K
Parameters
----------
GM : float
Product of the gravitational constant and the central mass
Rgas : float
Gas constant to link c_s and T through the ideal gas law. c_s = sqrt(gamma Rg/mu T). It's Rgas = kB/mH. kB = Boltzman constant. mH = mass of a hydrogen atom.
mu : float
Mean molecular weight
r0 : float
Reference radius
Sigma0 : float
Surface density at r0
betaSigma : float
Surface density powerlaw exponent
T0 : float
Temperature at r0
betaT : float
Temperature powerlaw exponent
gamma : float
Adiabatic index
betaOmega : float (optional)
Disk angular velocity powerlaw . Defaults to -3/2 for a Keplerian disk.
"""
GM: float
Rgas: float
mu: float
r0: float
Sigma0: float
betaSigma: float
gamma: float
T0: float = None
betaT: float = None
h0: float = None
betah: float = None
betaOmega: float = -3/2
def __post_init__(self):
"""
Set derived parameters
"""
self.OmegaK0 = np.sqrt(self.GM/self.r0**3)
is_T_set = self.T0 is not None
is_h_set = self.h0 is not None
if not (is_T_set or is_h_set):
raise ValueError("Either T0 or h0 must be set.")
elif is_T_set and is_h_set:
raise ValueError("Either T0 or h0 must be set, not both.")
elif is_T_set:
self.cs0 = np.sqrt(self.gamma*self.Rgas/self.mu*self.T0)
self.h0 = self.cs0 / np.sqrt(self.gamma) / self.OmegaK0
else:
self.cs0 = self.h0 * np.sqrt(self.gamma) * self.OmegaK0
self.T0 = self.mu/self.Rgas * self.cs0**2 / self.gamma
is_betaT_set = self.betaT is not None
is_betah_set = self.betah is not None
if not (is_betaT_set or is_betah_set):
raise ValueError("Either betaT or betah must be set.")
elif is_betaT_set and is_betah_set:
raise ValueError("Either betaT or betah must be set, not both.")
elif is_betaT_set:
self.betacs = self.betaT/2
self.betah = self.betaT/2 + 1/2
else:
self.betaT = 2*self.betah - 1
self.betacs = self.betaT/2
self.betaP = self.betaSigma + self.betaT
self.P0 = self.Sigma0 * self.cs0**2
def cs(self, r: np.ndarray):
return self.cs0 * (r/self.r0)**self.betacs
def h(self, r: np.ndarray):
return self.h0 * (r/self.r0)**self.betah
def Sigma(self, r):
return self.Sigma0 * (r/self.r0)**self.betaSigma
def T(self, r: np.ndarray):
return self.T0 * (r/self.r0)**self.betaT
def P(self, r: np.ndarray):
return self.P0 * (r/self.r0)**self.betaP
def OmegaK(self, r: np.ndarray):
return self.OmegaK0 * (r/self.r0)**self.betaOmega
def Omega(self, r: np.ndarray):
"""
Compute the background disk angular velocity profile.
Use a modification of SB15 Eq. (12) which uses the disk aspect ratio h.
Omega_0(r) = Omega_K(r) * sqrt{1 + beta_p * h^2}.
Parameters
----------
r : np.ndarray
Radial positions
Returns
-------
np.ndarray
"""
return self.OmegaK(r) * np.sqrt(1 + self.betaP * self.h(r)**2)
def plot(self, r: np.ndarray):
"""
Plot an overview of the disks powerlaws.
"""
mosaic_defs = [
["Sigma", "P", "T"],
["cs", "h", "Omega"]
]
import matplotlib.pyplot as plt
fig, axd = plt.subplot_mosaic(mosaic_defs, dpi=150, figsize=(12, 8))
axd["Sigma"].plot(r, self.Sigma(r))
axd["Sigma"].set_title("Sigma")
axd["Sigma"].set_xlabel("r")
axd["Sigma"].set_ylabel("Sigma")
axd["P"].plot(r, self.P(r))
axd["P"].set_title("P")
axd["P"].set_xlabel("r")
axd["P"].set_ylabel("P")
axd["T"].plot(r, self.T(r))
axd["T"].set_title("T")
axd["T"].set_xlabel("r")
axd["T"].set_ylabel("T")
axd["T"].set_xscale("log")
axd["T"].set_yscale("log")
axd["cs"].plot(r, self.cs(r))
axd["cs"].set_title("cs")
axd["cs"].set_xlabel("r")
axd["cs"].set_ylabel("cs")
axd["h"].plot(r, self.h(r))
axd["h"].set_title("h")
axd["h"].set_xlabel("r")
axd["h"].set_ylabel("h")
axd["Omega"].plot(r, self.Omega(r))
axd["Omega"].set_title("Omega")
axd["Omega"].set_xlabel("r")
axd["Omega"].set_ylabel("Omega")
for key, ax in axd.items():
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(which="minor", alpha=0.5,)
fig.tight_layout()
return fig, axd
Modify the initial conditions¶
To modify the data, we first define a helper class to handle the data files.
config = get_fargo_config(configfile)
outdir = "output/out"
disk = get_fargo_powerlaw_disk_output(outdir)
Adding a gap¶
Let’s add a gap by reducing the density and adjusting the azimuthal velocity.
gap_depth = 0.9
gap_width = 0.1
gap_center = 1.0
field = FargoCPTField("output/out", "0", "Sigma")
R = field.R
old_density = field.array
f = -gap_depth * np.exp(-0.5*(R-gap_center)**2/gap_width**2)
field.array = old_density * (1+f)
field.save()
field = FargoCPTField("output/out", "0", "vazi")
R = field.R
# need to recompute because R is different for vrad
f = -gap_depth * np.exp(-0.5*(R-gap_center)**2/gap_width**2)
# field.array = field.array * ( 1 - 2*R*(R-gap_center)/gap_width**2 * f/(1+f) * disk.cs(R)**2/(disk.gamma*(R*disk.Omega(R))**2) )**0.5
field.array = field.array * ( 1 - 2*R*(R-gap_center)/gap_width**2 * f/(1+f) * disk.h(R)**2 )**0.5
field.save()
Overview("output/out",
vars=["2:Sigma",
"1:Sigma",
"1:vazi:rel",
"1:vrad"]).create()
Now that the initial conditions have been adjusted, we can restart the simulation from snapshot number zero, which will pick up the modified conditions.
run(["auto", configfile], np=2, nt=1, exe=repo_root+"/bin/fargocpt_exe", detach=False)
Running command: mpirun -np 2 --report-pid /tmp/tmp5b2djodn -x OMP_NUM_THREADS=1 /home/rometsch/repo/fargocpt/bin/fargocpt_exe auto setup.yml
fargo process pid 1410404
[0] MPI rank # 0 runs as process 1410408
[0] MPI rank # 0 OpenMP thread # 0 of 1 on cpt-kamino
[1] MPI rank # 1 runs as process 1410409
[1] MPI rank # 1 OpenMP thread # 0 of 1 on cpt-kamino
[0] fargo: This file was compiled on Nov 14 2023, 12:56:40.
[0] fargo: This version of FARGO used _GNU_SOURCE
[0] fargo: This version of FARGO used NDEBUG. So no assertion checks!
[0] Using parameter file setup.yml
[0] Computing disk quantities within 5.00000e+00 L0 from coordinate center
[0] BC: Inner composite = reflecting
[0] BC: Outer composite = reflecting
[0] BC: Sigma inner = zerogradient
[0] BC: Sigma outer = zerogradient
[0] BC: Energy inner = zerogradient
[0] BC: Energy outer = zerogradient
[0] BC: Vrad inner = reflecting
[0] BC: Vrad outer = reflecting
[0] BC: Vaz inner = keplerian
[0] BC: Vaz outer = keplerian
[0] DampingTimeFactor: 1.00000e-01 Outer damping time is computed at radius of 2.50000e+00
[0] Damping VRadial to reference value at inner boundary.
[0] Damping VRadial to reference value at outer boundary.
[0] Damping VAzimuthal to reference value at inner boundary.
[0] Damping VAzimuthal to reference value at outer boundary.
[0] Damping SurfaceDensity to reference value at inner boundary.
[0] Damping SurfaceDensity to reference value at outer boundary.
[0] Damping Energy to reference value at inner boundary.
[0] Damping Energy to reference value at outer boundary.
[0] Radiative diffusion is disabled. Using fixed omega = 1.500000 with a maximum 50000 interations.
[0] Indirect Term computed as effective Hydro center acceleratrion with shifting the Nbody system to the center.
[0] Body force on gas computed via potential.
[0] Getting output number of snapshot 0
[0] Using FARGO algorithm for azimuthal advection.
[0] Using standard forward euler scheme for source terms.
[0] Cps is set, overwriting Nrad and Naz!
[0] Grid resolution set using cps = 2.000000
[0] The grid has (Nrad, Naz) = (74, 251) cells with (1.994113, 1.997395) cps.
[0] Computing scale height with respect to primary object.
[0] Using isothermal equation of state. AdiabaticIndex = 1.400.
[0] Viscosity is of alpha type with alpha = 1.000e-03
[0] Defaulting to VanLeer flux limiter
[0] Output information:
[0] Output directory: output/out/
[0] Number of files: 80
[0] Total output size: 0.00 GB
[0] Space Available: 31.23 GB
[0] Initializing 1 RNGs per MPI process.
[0] Warning : no `radii.dat' file found. Using default.
[0] The first 1 planets are used to calculate the hydro frame center.
[0] The mass of the planets used as hydro frame center is 1.000000e+00.
[0] 2 planet(s) initialized.
[0] Planet overview:
[0]
[0] # | name | mass [m0] | x [l0] | y [l0] | vx | vy |
[0] -----+-------------------------+------------+------------+------------+------------+------------+
[0] 0 | Star | 1 | 0 | -0 | 0 | 0 |
[0] 1 | Jupiter | 0.0009546033 | 1 | 0 | -0 | 1.000477 |
[0]
[0] # | e | a | T [t0] | T [a] | accreting | Accretion Type |
[0] -----+------------+------------+------------+------------+------------+----------------+
[0] 0 | 6.368246e-17 | 1 | 6.280188 | 0.999548 | 0 | No Accretion |
[0] 1 | 6.368246e-17 | 1 | 6.280188 | 0.999548 | 0 | No Accretion |
[0]
[0] # | Temp [K] | R [l0] | irradiates | rampuptime |
[0] -----+------------+------------+------------+------------+
[0] 0 | 5778 | 0.0046505 | yes | 0 |
[0] 1 | 0 | 4.6505e-05 | no | 0 |
[0]
[0] Using Tscharnuter-Winkler (1979) artificial viscosity with C = 1.410000.
[0] Artificial viscosity is used for dissipation.
[0] Surface density factor: 2.50663
[0] Tau factor: 0.5
[0] Tau min: 0.01
[0] Kappa factor: 1
[0] Minimum temperature: 2.81162e-05 K = 3.00000e+00
[0] Maximum temperature: 9.37206e+94 K = 1.00000e+100
[0] Heating from viscous dissipation is enabled. Using a total factor of 1.
[0] Cooling (beta) is disabled and reference temperature is floor. Using beta = 10.
[0] Cooling (radiative) is enabled. Using a total factor of 1.
[0] S-curve cooling is disabled.
[0] CFL parameter: 0.5
[0] Opacity uses tables from Lin & Papaloizou, 1985
[0] Particles are disabled.
[0] Initializing Sigma(r) = 2.25093e-05 = 200 g cm^-2 * [r/(1 AU)]^(-0.5)
[0] Total disk is mass is 0.000348841 = 6.9366e+29 g.
[0] Loading polargrinds for damping...
[0] Reading file 'output/out/snapshots/reference/Sigma.dat' with 148592 bytes.
[0] Reading file 'output/out/snapshots/reference/vrad.dat' with 150600 bytes.
[0] Reading file 'output/out/snapshots/reference/vazi.dat' with 148592 bytes.
[0] Loading polargrinds at t = 0...
[0] Reading file 'output/out/snapshots/0/Sigma.dat' with 148592 bytes.
[0] Reading file 'output/out/snapshots/0/vrad.dat' with 150600 bytes.
[0] Reading file 'output/out/snapshots/0/vazi.dat' with 148592 bytes.
[0] Restarting planetary system...
[0] Loading planets ...[0] done
[0] Loading rebound ...[0] done
[0] Finished restarting planetary system.
[0] Writing output output/out/snapshots/1, Snapshot Number 1, Time 6.280000.
[0] Writing output output/out/snapshots/2, Snapshot Number 2, Time 12.560000.
[0] Writing output output/out/snapshots/3, Snapshot Number 3, Time 18.840000.
[0] Writing output output/out/snapshots/4, Snapshot Number 4, Time 25.120000.
[0] Writing output output/out/snapshots/5, Snapshot Number 5, Time 31.400000.
[0] Writing output output/out/snapshots/6, Snapshot Number 6, Time 37.680000.
[0] Writing output output/out/snapshots/7, Snapshot Number 7, Time 43.960000.
[0] Writing output output/out/snapshots/8, Snapshot Number 8, Time 50.240000.
[0] Writing output output/out/snapshots/9, Snapshot Number 9, Time 56.520000.
[0] Writing output output/out/snapshots/10, Snapshot Number 10, Time 62.800000.
[0] -- Final: Total Hydrosteps 1728, Time 62.80, Walltime 4.01 seconds, Time per Step: 2.32 milliseconds
0
And finally, let’s see how this system evolved.
Uncomment the first line to get an interactive widget, if your jupyter installation supports it.
%matplotlib inline
Overview("output/out",
vars=["2:Sigma",
"1:Sigma",
"1:vazi:rel",
"1:vrad"]).create()
