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 = -k/\Gamma 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.

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 then compared this output to that astrobear

gave me.