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 with the Jeans instability test, 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

1. Sets up problem scales
2. Initializes the grid with the prescribed perturbations for Jeans instability
3. Sets boundary conditions on box
4. Finds max speed and uses CFL condition to get dt
5. Takes a hydro step; i.e. solves the local Riemann problem, gets fluxes, updates cells
6. Takes a source step, finds the mean of the density distribution, sets up the source term vector f=4 Pi G (rho-rhoavg), sends to the Jacobi solver, finds the solution for the gravitational potential phi, using phi, updates the velocity and energy of the cells
7. Updates time and checks whether tfinal has been reached
8. Sets boundary conditions on box
9. Gets dt
10. …

Update Formulas