Splitting Method

The non-homogenous 1D Euler equations with self-gravity are given by

\frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho u \\ E\end{pmatrix} + \frac{\partial }{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E+p) \end{pmatrix} = \begin{pmatrix} 0 \\ - \rho \triangledown \phi_g \\ -\rho u \cdot \triangledown \phi_g \end{pmatrix}

which in short hand notation is

\vec{u}_t + \vec{F}(\vec{u})_x = \vec{S}(\vec{u})

where \vec{u} is vector of fluid variables, \vec{F} is their fluxes, and \vec{S} is the source-term vector.

To solve this system of equations, we can employ operator splitting, which is analagous to the procedure for splitting the higher dimensional Euler equations (as I did for 2D — see wiki page).

It is basically as follows, first solve the homogenous equations,

\vec{u}_t + \vec{F}(\vec{u})_x = 0

with initial condition for the grid,

\vec{u}(x,0) = \vec{u}_i ^0

over a time-step dt (found by usual CFL condition for upwind Godunov scheme). This gives the solution,

\vec{u}^{n+1/2}

Next, solve the equation for the source-term (again over time interval dt),

\frac{d}{dt}\vec{u} = \vec{S}(\vec{u})

with initial condition

\vec{u}(x,0) = \vec{u}^{n+1/2}

This gives the solution for the complete time-step,

\vec{u}^{n+1}

Given the ODE's for the 'source' step only involve equations for momentum and energy, we conclude that only u and E change over this step (density does not).

Thus, a schematic is as follows,

Data time level n \rightarrow Hydro step \rightarrow

\vec{u}_i ^{n+1/2} = <\rho^{n+1/2}, u^{n+1/2}, E^{n+1/2}>\rightarrow

Source step \rightarrow \vec{u}_i ^{n+1} = <\rho^{n+1/2}=\rho^{n+1}, u^{n+1}, E^{n+1}>

Boundary Conditions

To test my code, I used periodic boundary conditions on the box and the poisson solver. For periodic BCs on the elliptic solver, this required subtracting off the mean density from the vector for rho, and using this as the source function.

Code Outline