# The Code

## pyro2

pyro (from PYthon hydRO) was originally written in 2003-2004 by Michael Zingale (SBU). This version was resurrected beginning in 2012 and rewritten for numpy and brought up to date. pyro2 is a 2D grid-based code mostly written in Python, with a few low-level routines written in Fortran. pyro2 employs a commonly used method to solve numerically the equations in grid-base codes, which is called the finite-volume method. Other methods of solving the equations include the finite-difference method and the finite-element method. Each problem that wants to be modeled is described by different equations. In pyro2 each problem includes a numerical solver according to the equations that we need to solve to describe the evolution of the system. Those solvers are classified in four categories:

• Advection: It models the transport of a substance or conserved property by a fluid at a particular speed. Mathematically can be modeled by the continuity equation. The details of the equations and methods used in this solver can be found on these notes. Below a video showing the advection of a high density region is shown as an example.
• Compressible: It models problems involving fluids in which the fluid density varies significantly in response to a change in pressure. In gases, a change in density is also accompanied by changes in temperature, which complicates considerably the analysis. The details of the equations and methods used in this solver can be found on these notes. Below a video showing the Kevin–Helmholtz instability is shown as an example.
• Diffusion: It models the transport without the needed of a bulk velocity. In simple terms a porperty can spread out in the surrounding medium. The details of the equations and methods used in this solver can be found on these notes. Below a video showing the diffusion of a high density region is shown as an example.
• Incompressible: It models problems involving fluids in whick the fluid density does not change due to change in pressure. As a rough approximation the density can be considered constant, which adds a constraint to the equations $\partial \rho / \partial t = 0$. The details of the equations and methods used in this solver can be found on these notes. Below a video showing the shear of an incompressible fluid is shown as an example.

Shock waves and choked flow are only developed in compressible fluids, which is a key difference between compressible and incompressible flows. Given that, a little more detail of the compressible solver will be given in the next section.

## Solver for compressible problems

Compressible fluids are governed by the ‎Euler equations. These equations are solved using a second order Godunov method. The Godunov’s scheme is a finite-volume method whereby the fluxes through the interfaces are computed by solving the Riemann problem for the system:

$U^{n+1}_i = U^n_i + {\Delta t \over \Delta x} \left( F^{n+1/2}_{i-1/2} - F^{n+1/2}_{i+1/2} \right)$

This says that each of the conserved quantities in $U$ change only due to the flux $F$ of that quantity through the boundary of the cell.

Constructing these interface states requires reconstructing the cell-average data and doing characteristic tracing to see how much of each characteristic quantity comes to the interface over $\Delta t/2$. Since we are comfortable working with the primitive variables, the reconstruction is done on those, and they are then projected into the characteristic variables to determine how much of each characteristic quantity is carried by each of the 3 waves. There are several methods to reconstruct these quantities, among them we can find the piecewise constant, linear, and parabolic polynomial. Below we give more details about the piecewise linear method which is the one is used in the compressible solver.

### Piecewise Linear Reconstruction

In this method the cell average data is approximated by a line with non-zero slope within each cell. (see Fig. 2) But the reconstruction is done using the primitive quantity $q$ associated to the conserved quantity $U$ in order to use a simpler structure.

For example consider the left state at the interface $i + 1/2$ (see Fig. 1). We do a Taylor expansion through $\Delta x/2$ to bring us to the interface, and $\Delta t/2$ to bring us to the midpoint in time. We then obtain:

$q^{n+1/2}_{i+1/2,L} = q^n_i+{1\over 2}(I-A{\Delta t \over \Delta x})\Delta q_i$

where $q_i$ is the cell-centered primitive variable and $\Delta q_i$ is the reconstructed slope of the primitive variable in that cell. The dropped terms in the Taylor expansion are $O(\Delta x^2)$ and $O(\Delta t^2)$ resulting in a second-order accurate method in space and time.
Analogously for the right state at the interface $i + 1/2$, the formula is given by:

$q^{n+1/2}_{i+1/2,R} = q^n_{i+1}-{1\over 2}(I+A{\Delta t \over \Delta x})\Delta q_{i+1}$

As a final note, it is important to note that this solver uses the Rankine-Hugoniot jump conditions to fully resolve shocks.