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}>

Update Formulas

To solve the first source term ODE,

\frac{d}{dt}(\rho u) = - \rho \triangledown \phi

replace with discrete derivatives,

\frac{(\rho u)_i^{n+1} - (\rho u)_i ^{n+1/2}}{dt} = -\frac{\rho_i^{n+1/2}(\phi_{i+1}^{n+1/2}- \phi_i^{n+1/2})}{dx}

which becomes,

\frac{\rho_i^{n+1}(u_i^{n+1} - u_i^{n+1/2})}{dt} = -\frac{\rho_i^{n+1/2}(\phi_{i+1}^{n+1/2}- \phi_i^{n+1/2})}{dx}

since

\rho_i^{n+1} = \rho_i^{n+1/2}

we can just cancel the rho's and solve for the updated u,

u_i^{n+1}=-\frac{dt}{dx}[\phi_{i+1}^{n+1/2}-\phi_i^{n+1/2} ] + u_i^{n+1/2}

This requires a dt (which is the same for the hydro step, solved using the CFL condition), and a gravitational potential phi. Phi is solved using the Jacobi solver, with the source function given by

f(i)=4 \pi G (\rho_i^{n+1/2} - \bar{\rho})

and boundary conditions

\phi(0) = \phi(mx), ~\phi(mx+1) = \phi(1)

To solve the 2nd source term ODE,

\frac{d}{dt}(E) = -\rho u \cdot \triangledown \phi

replace with discrete derivatives,

\frac{E_i^{n+1} - E_i^{n+1/2}}{dt} = -\frac{\rho_i^{n+1/2}u_i^{n+1/2}(\phi_{i+1}^{n+1/2} - \phi_i^{n+1/2})}{dx}

I am not exactly sure which time level to use for u (n+1 or n+½).

Rearranging for the updated version of the energy, we have,

E_i^{n+1} = -\frac{dt}{dx}[\rho_i^{ n+1/2} u_i^{n+1/2}(\phi_{i+1}^{n+1/2}-\phi_i^{n+1/2})] + E_i^{n+1/2}