Hydrodynamic Simulation
The hydrodynamic evolution of a jet with a given angular profile can be solved by creating an object jetsimpy.Jet. Below is
simulation
- class jetsimpy.Jet(profile, nwind, nism, tmin=10, tmax=3200000000.0, grid=jetsimpy.Uniform(257), tail=True, spread=True, cal_level=1, rtol=1e-06, cfl=0.9)
Create an object containing the simulation results.
- Parameters:
profile (tuple) – a tuple of the tabulated data of the jet angular profile taking the form of (theta, Eiso, lf). theta is the polar angle array starting from 0 and end with pi. Eiso is the isotropic equivalent energy. lf is the Lorentz factor. All three elements must be 1D np.array of the same length.
nwind (float) – wind density scale.
nism (float) – ism density scale.
tmin (float) – start time of the simulation.
tmax (float) – end time of the simulation.
grid (np.array) – the polar angles of the “cell edges” of the simulation. This parameter specifies the resolution of the simulation.
tail (bool) – Add an low-energy isotropic tail or not. This tail ensures that the energy is positive and the Lorentz factor is above unity.
spread (bool) – turning on the spreading or not.
cal_level (int) – calibration level. “cal_level=0”: no calibration. “cal_level=1”: calibrate with Blandford-McKee solutions throughout the simulation. “cal_level=2”: calibrate with Blandford-McKee solution in the relativistic limit and Sedov-Taylor solution in the Newtonian limit. Setting cal_level=2 might be aggressive since the thin shell approximation is no longer valid in the Newtonian limit.
rtol (float) – relative tolerance in the primitive variable solver. Never change this parameter unless you know what is going on.
cfl (float) – the CFL number. Never change this parameter unless you know what is going on.
jet profile
The jet profile can be arbitrary but must be physically reasonable. For general purpose, users should always provide tabulated data even if the profile is analytic. To save time, the code provides some built-in profiles.
- jetsimpy.TopHat(theta_c, Eiso, lf0=1e+100)
A top-hat jet profile.
- jetsimpy.Gaussian(theta_c, Eiso, lf0=1e+100)
A Gaussian jet profile defined by
\[ \begin{align}\begin{aligned}E(\theta) &= E_{\rm iso} \exp\left[{-\frac{1}{2}\left(\frac{\theta}{\theta_{\rm c}}\right)^2}\right]\\\Gamma(\theta) &= (\Gamma_0 - 1)\exp\left[{-\frac{1}{2}\left(\frac{\theta}{\theta_{\rm c}}\right)^2}\right] + 1\end{aligned}\end{align} \]
- jetsimpy.PowerLaw(theta_c, Eiso, lf0=1e+100, s=4.0)
A power-law jet profile defined by
\[ \begin{align}\begin{aligned}E(\theta) &= E_{\rm iso} \left[1 + \left(\frac{\theta}{\theta_{\rm c}}\right)^2\right]^{-s/2}\\\Gamma(\theta) &=(\Gamma_0 - 1)\left[1 + \left(\frac{\theta}{\theta_{\rm c}}\right)^2\right]^{-s/2} + 1\end{aligned}\end{align} \]
coasting phase
To turn off the coasting phase, just set a very large (but finite) initial Lorentz factor such as \(10^{100}\).
external density
The external density is assumed to follow the profile below.
start and end time
The start time tmin is non-zero because the initial radius must be positive. The blastwave will coast at the initial speed until tmin, and then the simulation starts. Physically tmin can be set to a value around the photospherical radius (timescale). A moderate value such as tmin=10 is preferred.
The simulation ends at tmax. This value should be sufficiently large to cover the time period of afterglow modeling.
resolution
In a finite volume method, the simulation domain is discretized into a number of cells. In this code the discretization is determined by the grid parameter. The grid parameter is a 1D numpy array of the “cell edges”. Namely, the cell centers locate at the middle of adjacent elements.
The resolution directly determines the speed of the code. The optimal way to choose the resolution is to set a high resolution inside the jet core and low resolution outside the core. However, to maintain numerical stability, the size of adjacent cells shouldn’t change dramatically.
The code provides some built-in resolution setup for jet profiles with well-defined opening angles.
- jetsimpy.Uniform(theta_c, npoints)
Uniform cell distribution. A conservative choice for all cases.
- jetsimpy.ForwardJetRes(theta_c, npoints)
Resolution optimized for forward-jet only.
- jetsimpy.CounterJetRes(theta_c, npoints)
Resolution optimized for counter-jet only.
simulation results
One can access the simulation data from the jetsimpy.Jet object.
- property .t_pde
time data (1d numpy.array)
- property .y_pde
primitive variable data with the shape (5, ntheta, ntime). The 5 variables are [\(M_{\rm sw}\), \(M_{\rm ej}\), \(\beta^2\gamma^2\), \(\beta_{\theta}\), \(R\)].
- property .theta_pde
array of cell centers
- .dE0_dOmega(t, theta)
The interpolated energy profile (rest mass excluded).
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle
- .dMsw_dOmega(t, theta)
The interpolated swept-up mass profile.
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle
- .dMej_dOmega(t, theta)
The interpolated ejecta mass profile.
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle
- .beta_gamma(t, theta)
The interpolated 4-velocity profile.
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle
- .beta_theta(t, theta)
The interpolated tangent velocity profile.
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle
- .R(t, theta)
The interpolated Radius angular profile.
- Parameters:
t (numpy.array) – time
theta (numpy.array) – polar angle