.. pyrodeo documentation master file, created by sphinx-quickstart on Fri Jun 1 21:12:37 2018. You can adapt this file completely to your liking, but it should at least contain the root `toctree` directive. Welcome to pyrodeo! =================================== .. toctree:: :maxdepth: 2 :caption: Contents: Pyrodeo is a Python implementation of RODEO (ROe solver for Disc-Embedded Objects), a Roe solver implementation aimed at hydrodynamic simulations of astrophysical discs. Installation ------------------------------------- If Python is not installed, download from `here `_ and install. The latest versions of Python come with package manager pip included. Then Pyrodeo can be installed from the command line simply by entering:: pip install pyrodeo Quick start ------------------------------------- Within Python, first import the simulation module:: >>> import pyrodeo Create a simulation in Cartesian geometry with standard parameters:: >>> sim = pyrodeo.Simulation.from_geom('cart') Run the simulation up to t=0.25:: >>> sim.evolve([0.25]) Since the standard initial conditions consist of constant density and pressure and zero velocity, no visible evolution takes place. For more interesting examples, see :ref:`Examples`. Equations solved ------------------------------------- The current version supports inviscid isothermal hydrodynamics in three spatial dimensions. Isothermal means the pressure :math:`p` is related to the density :math:`\rho` simply through :math:`p=c^2\rho`, where the sound speed :math:`c` is either a constant (fully isothermal) or a prescribed function of position (locally isothermal). Four geometries are available: Cartesian coordinates, the shearing sheet, cylindrical coordinates and spherical coordinates. Cartesian coordinates ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In Cartesian coordinates, we have the continuity equation: .. math:: \frac{\partial\rho}{\partial t} + \frac{\partial}{\partial x}(\rho v_x) + \frac{\partial}{\partial y}(\rho v_y) + \frac{\partial}{\partial z}(\rho v_z)=0 Momentum in x-direction: .. math:: \frac{\partial}{\partial t}(\rho v_x) + \frac{\partial}{\partial x}(\rho v_x^2 + c^2\rho) + \frac{\partial}{\partial y}(\rho v_x v_y) + \frac{\partial}{\partial z}(\rho v_x v_z)=0 Momentum in y-direction: .. math:: \frac{\partial}{\partial t}(\rho v_y) + \frac{\partial}{\partial x}(\rho v_x v_y) + \frac{\partial}{\partial y}(\rho v_y^2 + c^2\rho) + \frac{\partial}{\partial z}(\rho v_y v_z)=0 Momentum in z-direction: .. math:: \frac{\partial}{\partial t}(\rho v_z) + \frac{\partial}{\partial x}(\rho v_x v_z) + \frac{\partial}{\partial y}(\rho v_y v_z)+ \frac{\partial}{\partial z}(\rho v_z^2 + c^2\rho) =0 Shearing Sheet ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Shearing Sheet is essentially a Cartesian model of a small patch in an astrophysical disc. This patch is rotating at the local Keplerian angular velocity :math:`\Omega`, which means that Coriolis and centrifugal-type forces need to be included on the right-hand side of the equations. On the other hand, the patch is assumed to be small enough so that a local Cartesian frame can be used in stead of cylindrical coordinates. Usually the computational domain is taken to be periodic in y and shear-periodic in x (periodic but corrected for the shear). We therefore still have the continuity equation: .. math:: \frac{\partial\rho}{\partial t} + \frac{\partial}{\partial x}(\rho v_x) + \frac{\partial}{\partial y}(\rho v_y) + \frac{\partial}{\partial z}(\rho v_z)=0 The x-momentum equation now includes source terms on the right-hand side: .. math:: \frac{\partial}{\partial t}(\rho v_x) + \frac{\partial}{\partial x}(\rho v_x^2 + c^2\rho) + \frac{\partial}{\partial y}(\rho v_x v_y) + \frac{\partial}{\partial z}(\rho v_x v_z)=2\Omega\rho v_y + 3\rho\Omega^2 x Same for the momentum equation in y-direction: .. math:: \frac{\partial}{\partial t}(\rho v_y) + \frac{\partial}{\partial x}(\rho v_x v_y) + \frac{\partial}{\partial y}(\rho v_y^2 + c^2\rho) + \frac{\partial}{\partial z}(\rho v_y v_z)=-2\Omega\rho v_x In the z-direction we get a source term due to the vertical component of the stellar gravity: .. math:: \frac{\partial}{\partial t}(\rho v_z) + \frac{\partial}{\partial x}(\rho v_x v_z) + \frac{\partial}{\partial y}(\rho v_y v_z)+ \frac{\partial}{\partial z}(\rho v_z^2 + c^2\rho) =-\rho \Omega^2 z .. NOTE:: In the shearing sheet the sound speed should really be constant (no locally isothermal shearing sheet). Together, sound speed and angular velocity define a length scale :math:`c/\Omega`, which is a measure of the scale height of the disc. Typically one chooses :math:`c=\Omega=1`, so that distances are measured in scale heights and time in inverse orbital frequency. Cylindrical coordinates ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For a full disc in two dimensions cylindrical coordinates :math:`(R,\varphi,z)` are preferred. This time we have geometrical source terms and gravity from the central object to worry about. The continuity equation now reads: .. math:: \frac{\partial}{\partial t}(R\rho) + \frac{\partial}{\partial R}(R\rho v_R) + \frac{\partial}{\partial \varphi}(\rho v_\varphi) + \frac{\partial}{\partial z}(R\rho v_z)=0 The radial momentum equation now includes source terms representing centrifugal and gravitational forces, in addition to a geometrical pressure source term: .. math:: \frac{\partial}{\partial t}(R\rho v_R) + \frac{\partial}{\partial R}(R\rho v_R^2 + c^2R\rho) + \frac{\partial}{\partial \varphi}(\rho v_R v_\varphi)+\frac{\partial}{\partial z}(R\rho v_R v_z)= \rho v_\varphi^2 - R^2\rho\frac{GM_*}{(R^2+z^2)^{3/2}} + c^2\rho In the :math:`\varphi` direction we get a Coriolis source term: .. math:: \frac{\partial}{\partial t}(\rho v_\varphi) + \frac{\partial}{\partial R}(\rho v_R v_\varphi) + \frac{1}{R}\frac{\partial}{\partial \varphi}(\rho v_\varphi^2 + c^2\rho) +\frac{\partial}{\partial z}(\rho v_\varphi v_z)=-2\rho v_R v_\varphi/R In the vertical direction we again have the vertical component of the stellar gravity: .. math:: \frac{\partial}{\partial t}(R\rho v_z) + \frac{\partial}{\partial R}(R\rho v_R v_z) + \frac{\partial}{\partial \varphi}(\rho v_\varphi v_z)+ \frac{\partial}{\partial z}(R\rho v_z^2 + c^2R\rho) = -R\rho z\frac{GM_*}{(R^2+z^2)^{3/2}} .. NOTE:: The unit of mass is taken to be the mass of the central object. The unit of distance is some reference radius. The unit of time is the inverse angular velocity at the reference radius. In this system of units, the gravitational constant is unity, and one orbit equals :math:`2\pi` time units. Spherical coordinates ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For a full disc in three dimensions spherical coordinates :math:`(r,\theta,\varphi)` are often preferred. The continuity and momentum equations now read: .. math:: \frac{\partial}{\partial t}(r^2\sin\theta\rho) + \frac{\partial}{\partial r}(r^2\sin\theta\rho v_r) + \frac{\partial}{\partial \varphi}(r\rho v_\varphi) + \frac{\partial}{\partial \theta}(r\sin\theta\rho v_\theta)=0 .. math:: \frac{\partial}{\partial t} (r^2\sin\theta\rho v_r) +\frac{\partial}{\partial r}(r^2\sin\theta\rho(v_r^2+c^2)) + \frac{\partial}{\partial\theta}(r\rho v_r v_\theta\sin\theta) + \frac{\partial}{\partial\varphi} (r\rho v_r v_\varphi) = \\ r^2\sin\theta\rho\frac{v_\theta^2 + v_\varphi^2}{r} -r^2\sin\theta\rho\frac{\partial \Phi}{\partial r}+2r\sin\theta c^2\rho .. math:: \frac{\partial}{\partial t} (r^2\sin\theta\rho v_\theta) + \frac{\partial}{\partial r}(r^2\sin\theta\rho v_rv_\theta) + \frac{\partial}{\partial\theta}(r\sin\theta(\rho v_\theta^2 + p)) + \frac{\partial}{\partial\varphi}(r\rho v_\varphi v_\theta) = \\ r\rho v_\varphi^2\cos\theta - r\sin\theta\rho v_rv_\theta + r\cos\theta p .. math:: \frac{\partial}{\partial t} (r^2\sin\theta\rho v_\varphi) + \frac{\partial}{\partial r}(r^2\sin\theta\rho v_r v_\varphi) + \frac{\partial}{\partial\theta}(r\sin\theta\rho v_\theta v_\varphi) + \frac{\partial}{\partial\varphi}(r\rho (v_\varphi^2+c^2)) = \\ -r\sin\theta\rho v_\varphi v_r-r\rho v_\varphi v_\theta \cos\theta Extra source terms ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Pyrodeo solves inviscid isothermal hydrodynamics, and in the shearing sheet and cylindrical and spherical geometries only gravity from the central object is considered. Extra physics, as far as it concerns extra source terms, can be added by a user-defined source integration function. See the :ref:`Examples` section. This function is called once per time step and can also be used for monitoring various quantities (mass, torque on planet, etc.). Numerical method ------------------------------- Dimensional splitting ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Pyrodeo uses dimensional splitting to integrate the equations. First direction (x, R, r) """"""""""""""""""""""""""" For the x direction (therefore neglecting y- and z-derivatives), we can cast the equations into the following form: .. math:: \frac{\partial\bar\rho}{\partial t} + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x) = 0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_x) + \frac{\partial}{\partial \bar x}(\bar\rho \bar v_x^2 + \bar c^2\bar \rho) = S_x .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_y) + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x \bar v_y) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_z) + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x \bar v_z) =0 For the Cartesian setup, we simply have .. math:: \bar\rho =\rho, \bar v_x = v_x, \bar v_y = v_y, \bar v_z = v_z, \bar c = c, \bar x = x, \bar y = y, \bar z = z, S_x = 0 For the shearing sheet, we need .. math:: \bar\rho =\rho, \bar v_x = v_x, \bar v_y = v_y - \Omega x/2, \bar v_z = v_z, \bar c = c, \bar x = x, \bar y = y, \bar z = z , S_x = 2\Omega \rho (v_y +3\Omega x/2) In cylindrical coordinates we need .. math:: \bar\rho =R\rho, \bar v_x = v_R, \bar v_y = Rv_\varphi, \bar v_z = v_z, \bar c = c, \bar x = R, \bar y = \varphi, \bar z = z , S_x = \rho v_\varphi^2 - \frac{R^2\rho GM_*}{(R^2+z^2)^{3/2}} + c^2\rho Finally, for spherical coordinates we need: .. math:: \bar\rho =r^2\sin\theta \rho, \bar v_x =v_r, \bar v_y =rv_\varphi , \bar v_z =rv_\theta, \bar c = c, \bar x = r, \bar y = \varphi, \bar z = \theta, S_x = \bar\rho\frac{\bar v_z^2 + \bar v_y^2}{r^3}-\bar\rho\frac{GM_*}{r^2}+\frac{2 c^2\bar\rho}{r} Second direction (y, :math:`\varphi`) """"""""""""""""""""""""""""""""""""""""""""""" For the y-integration (neglecting x- and z-derivatives) we can cast the equations in the form: .. math:: \frac{\partial\bar\rho}{\partial t} + \frac{\partial}{\partial \bar y}(\bar\rho\bar v_y) = 0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_x) + \frac{\partial}{\partial \bar y}(\bar\rho\bar v_x \bar v_y) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_y) + \frac{\partial}{\partial \bar y}(\bar\rho \bar v_y^2 + \bar c^2\bar \rho) = S_y .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_z) + \frac{\partial}{\partial \bar y}(\bar\rho\bar v_y \bar v_z) =0 For both the Cartesian setup and the shearing sheet, we simply have .. math:: \bar\rho =\rho, \bar v_x = v_x, \bar v_y = v_y, \bar v_z = v_z, \bar c = c, \bar x = x, \bar y = y, \bar z = z, S_y = 0 In cylindrical coordinates we need .. math:: \bar\rho =\rho, \bar v_x = v_R, \bar v_y = v_\varphi/R, \bar v_z = v_z, \bar c = c/R, \bar x = R, \bar y = \varphi, \bar z = z, S_y = 0 Finally, spherical coordinates: .. math:: \bar\rho =\rho, \bar v_x =v_r, \bar v_y =v_\varphi/(r\sin\theta), \bar v_z =v_\theta, \bar c =c/(r\sin\theta), \bar x =r, \bar y =\varphi, \bar z =\theta, S_y =0 Third direction (z, :math:`\theta`) """"""""""""""""""""""""""""""""""""""""""""""" Finally, for the z integration we can cast the equations in the form: .. math:: \frac{\partial\bar\rho}{\partial t} + \frac{\partial}{\partial \bar z}(\bar\rho\bar v_z) = 0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_x) + \frac{\partial}{\partial \bar y}(\bar\rho\bar v_x \bar v_z) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_y) + \frac{\partial}{\partial \bar z}(\bar\rho\bar v_y \bar v_z) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_z) + \frac{\partial}{\partial \bar z}(\bar\rho \bar v_z^2 + \bar c^2\bar \rho) = S_z For both the Cartesian setup, we simply have .. math:: \bar\rho =\rho, \bar v_x = v_x, \bar v_y = v_y, \bar v_z = v_z, \bar c = c, \bar x = x, \bar y = y, \bar z = z, S_z = 0, with the shearing sheet being exactly the same but with a non-zero source :math:`S_z=-\rho\Omega^2 z`. In cylindrical coordinates we need .. math:: \bar\rho =\rho, \bar v_x = v_R, \bar v_y = v_\varphi/R, \bar v_z = v_z, \bar c = c, \bar x = R, \bar y = \varphi, \bar z = z, S_z =-\rho z\frac{GM_*}{(R^2+z^2)^{3/2}} Finally, spherical coordinates: .. math:: \bar\rho =\sin\theta \rho, \bar v_x =v_r, \bar v_y =\sin\theta v_\varphi/r, \bar v_z =v_\theta/r, \bar c =c/r, \bar x =r, \bar y =\varphi, \bar z =\theta, S_z =\bar\rho \cot\theta (\bar v_y^2/\sin^2\theta+\bar c^2) Unified approach """"""""""""""""""""""""""" Note that the resulting equations are very symmetric in x, y and z: if we swap x and y in the y integration the equations have exactly the same form as for the x integration. Similar for the z integration when swapping z and x. Therefore, if we prepare all quantities appropriately, we only need a single hydrodynamic solver that is able to advance the system .. math:: \frac{\partial\bar\rho}{\partial t} + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x) = 0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_x) + \frac{\partial}{\partial \bar x}(\bar\rho \bar v_x^2 + \bar c^2\bar \rho) = S_x .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_y) + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x \bar v_y) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_z) + \frac{\partial}{\partial \bar x}(\bar\rho\bar v_x \bar v_z) =0 This is what is done in the :class:`.Roe` class. The necessary preparation is done in the :class:`.Hydro` class. Orbital advection ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In the case of the shearing sheet and the cylindrical disc there is a large :math:`\bar v_y` that severely limits the time step. This limit can be overcome by splitting the y integration one more by writing :math:`\bar v_y = \bar u_y + v_0` where :math:`v_0` is independent of y .. math:: \frac{\partial\bar\rho}{\partial t} + v_0\frac{\partial\bar\rho}{\partial\bar y} + \frac{\partial}{\partial \bar y}(\bar\rho\bar u_y) = 0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar v_x) + v_0\frac{\partial}{\partial\bar y} (\bar\rho \bar v_x)+ \frac{\partial}{\partial \bar y}(\bar\rho\bar v_x \bar u_y) =0 .. math:: \frac{\partial}{\partial t}(\bar\rho \bar u_y) + v_0\frac{\partial}{\partial\bar y} (\bar\rho \bar u_y)+ \frac{\partial}{\partial \bar y}(\bar\rho \bar u_y^2 + \bar c^2\bar \rho) = 0 The terms involving :math:`v_0` make up a linear advection problem that can be solved straightforwardly for any time step. This is done in the :class:`.LinearAdvection` class. The remaining terms are integrated in the :class:`.Roe` class, but at a much larger time step because presumably :math:`u_y \ll v_0`. Algorithm overview ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A single time step in :func:`Hydro.evolve ` consists of the following steps: 1. Calculate time step using :func:`Hydro.calc_time_step `. 2. Set shear periodic boundary conditions if necessary. 3. Preprocessing step to cast the equations in the same form for all geometries and directions, while at the same time calculating the source term using :func:`Hydro.preprocess `. 4. Use the Roe solver to advance the hydrodynamic equations using :func:`Roe.step `. 5. Do orbital advection if necessary using :func:`Hydro.orbital_advection `. 6. Do the inverse of step 3, getting all quantities back to their original form in :func:`Hydro.postprocess `. 7. Integrate any extra source terms. Boundary conditions ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The available boundary conditions are * 'closed': closed boundary, i.e. no mass flow through the boundary. Waves will be reflected off the boundary. * 'periodic': periodic boundary. * 'nonreflecting': allow waves to pass through unhindered. * 'symmetric': assume boundary is a symmetry plane. Very much like a closed boundary, but less general in that it needs the boundary to be a plane of symmetry. By default, all boundaries are set to be closed; this can be changed by changing the `boundaries` attribute of the :class:`Param ` class. Output ------------------------------- Once the integration routine :func:`Simulation.evolve ` has finished the final state is available through :func:`simulation.state`. In addition, an output file `rodeo.h5` is created containing the state at all specified checkpoints. This is an `HDF5 `_ file created with `h5py `_. It contains the following groups: * param: Simulation parameters as specified in the :class:`Param ` class. * coords: Coordinates from the :class:`Coordinates ` class. * checkpoint#: State at checkpoint, where # stands for an integer. .. NOTE:: Both state and coordinate arrays include two ghost zones on each side in all directions. This is in order to be able to restore a simulation from a checkpoint. .. NOTE:: The value stored in `state.vely` is the y-velocity with the orbital advection velocity removed! In other words, the equilibrium solution in a constant pressure shearing sheet or cylindrical disc has vanishing `state.vely`. An example of reading the file and plotting using `matplotlib `_: .. literalinclude:: ../examples/example_plot.py .. _Examples: Examples ------------------------------- Shock tube ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A simple example in Cartesian geometry is a one-dimensional isothermal shock tube: .. literalinclude:: ../examples/example_cartesian.py The standard grid dimensions are `(100,1)`, which means 100 cells in x and 1 in y. Try a higher resolution by explicitly specifying the dimensions in :func:`Simulation.from_geom `. Instability in shearing sheet ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A more demanding two-dimensional calculation involves the instability of a sharp density ridge in the shearing sheet: .. literalinclude:: ../examples/example_sheet.py The checkpoints have been chosen close enough together to allow for the results to be animated: .. literalinclude:: ../examples/example_animation.py Disc-planet interaction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ As the final, most complex example, consider a planet embedded in a disc in cylindrical coordinates. Since pyrodeo only includes gravity from the central star, we need to provide an extra source term to account for the gravitational force due to the planet. In addition, we define wave-killing zones on the radial edges of the domain to avoid wave reflection. It will take some time to run this simulation, so have a cup of tea and come back to see a Jupiter-like planet carve out a gap in the disc. .. literalinclude:: ../examples/example_planet.py Class reference ========================= .. automodule:: pyrodeo Coordinates ------------------------------- .. automodule:: pyrodeo.coords :members: State ------------------------------- .. automodule:: pyrodeo.state :members: Param ------------------------------- .. automodule:: pyrodeo.param :members: Conservation law solver ------------------------------- .. automodule:: pyrodeo.claw_solver :members: Linear advection solver ------------------------------- .. automodule:: pyrodeo.linear_advection :members: Roe solver ------------------------------- .. automodule:: pyrodeo.roe :members: Hydro ------------------------------- .. automodule:: pyrodeo.hydro :members: Simulation ------------------------------- .. automodule:: pyrodeo.simulation :members: