Initialization

To test the operator splitting, I initialized the grid with the same conditions as the Jeans test I did with astrobear (see that wiki page under the self-gravity section of my user page). Namely, I computed the Jeans length (in cm) for a gas of density rho = 3.34d-21 g/cc, temperature T = 100 K, and gamma = 1.0001. I then seeded the grid (that was 50x Jeans length long) with the prescribed density, velocity, and pressure perturbations as derived before. The density form was straightforward,

\rho = \rho_0 / rscale + Cos(kx)

where rscale = rho_0, and k is the wavenumber of the Jeans mode.

The background pressure (P0) and pressure perturbation (P1) was given by the ideal gas law,

P = P_0 + P_1 = \frac{K_B T}{m_H}(\rho_0 + \rho_1) = 1 + \frac{K_B T}{m_H}(\rho_1) / pscale

The velocity formula was a little more involved, looking like

u = -\frac{\Gamma}{k \rho_0} Sin(kx)

where the growthrate is given by,

\Gamma = 4 \pi G \rho_0[1-(\lambda_J/ \lambda)^2]

In order to keep the velocity perturbation dimensionless, the factor out in front of the Sin(kx) needed to be scaled to computational units. This involved a factor of lscale for k, rhoscale for rho_0 in the growthrate, and a term for scaled gravity in the growthrate.

Scaling gravitational constant

To get a form for 'scalegrav', I noted I needed a scale factor for time, given the units of G are

[G] = \frac{cm^3}{g s^2} = \frac{1}{\rho s^2}

To derive a scale for time, I needed to use the 3 scales I had already defined: lscale, rhoscale, pscale. I did this with the following equations,

P=F/A=(\frac{m L}{T^2})(\frac{1}{L^2})=\frac{M}{T^2 L} = \frac{\rho L^2}{T^2}

which gives an expression for the relationship between the scales. Rearranging and solving for timescale, we have

timescale = (\frac{rhoscale}{pscale})^{1/2} Lscale

Thus, the scaled version of the gravitational constant is given by,

scalegrav = 6.674^{-8}\frac{cm^3}{g s^2}rhoscale (timescale)^2

With all quantities properly scaled, the initial condition (in computational units) looks like,

density:

velocity:

and pressure:

With these initial conditions for the primitive variables of the hydro solver, I then began time-stepping through the solution, using periodic boundaries on both my box and poisson solver (​see).

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. …

Results

To get the final time (in computational units), I had the code print out tfinal=5/growthrate, and used this in my problem.data for the simulation time. I defined a new array to send to the poisson solver for the source function, that was the density array with the mean density subtracted off.