## C  Exact solution details

### C.1  Errors

The error norms calculated when exact solutions are plotted are as follows: The error for each particle is given by

 ei = fi − fexact,     (9)

where the exact solution fexact(x) is the solution returned from the exact solution subroutines (with resolution adjustable in the exact solution options menu option) interpolated to the position of the current particle xi via a simple linear interpolation. The absolute L1 error norm is simply the average of the errors across the domain, calculated according to

∥ e ∥L1 =
 1 N fmax

 N ∑ i=1
| ei |,     (10)

where fmax is the maximum value of the exact solution in the region in which the particles lie (also only particles in the current plot are used) which is used to normalise the error estimate. A better error norm is the L2 or Root Mean Square (RMS) norm given by

∥ e ∥L2 =

 1 N

 1 fmax2

 N ∑ i=1
| ei |2

 1/2
.     (11)

Finally the maximum error, or L norm is calculated according to

∥ e ∥L =
 1 fmax
maxi | ei |.     (12)

which is the most stringent error norm.

The inset plot of the individual particle errors shows the fractional deviation for each particle given by

 ei,frac = (fi − fexact) / fexact.     (13)

### C.2  Shock tubes (Riemann problem)

The subroutine `exact_shock` plots the exact solution for a one-dimensional shock tube (Riemann problem). The difficult bit of the problem is to determine the jump in pressure and velocity across the shock front given the initial left and right states. This is performed in a separate subroutine (riemannsolver) as there are many different methods by which this can be done (see e.g. ???). The actual subroutine exact_shock reconstructs the shock profile (consisting of a rarefaction fan, contact discontinuity and shock, summarised in Figure 19), given the post-shock values of pressure and velocity.

The speed at which the shock travels into the ‘right’ fluid can be computed from the post shock velocity using the relation

vshock = vpost
 (ρpost/ρR) (ρpost/ρR)− 1
,     (14)

where the jump conditions imply

 ρpost ρR
=
 (Ppost/PR) + β 1 + β (Ppost/PR)
(15)

with

β =
 γ − 1 γ + 1
.     (16)

#### C.2.1   Riemann solver

The algorithm for determining the post-shock velocity and pressure is taken from ???.

### C.3  Polytrope

The subroutine `exact_polytrope` computes the exact solution for a static polytrope with arbitrary γ. From Poisson’s equation

 ∇2 φ = 4π G ρ,     (17)

assuming only radial dependence this is given by

 1 r2

 d dr

r2
 dφ dr

= 4π G ρ(r).     (18)

The momentum equation assuming an equilibrium state (v = 0) and a polytropic equation of state P = Kργ gives

 dφ dr
= −
 γ K γ−1
 d dr

ρ(γ −1)
(19)

Combining (18) and (19) we obtain an equation for the density profile

 γ K 4π G (γ − 1)

 1 r2

 d dr

r2
 d dr

ργ−1

+ ρ(r) = 0.     (20)

This equation can be rearranged to give

 γ K 4π G (γ − 1)

 d2 dr2

rργ−1
rρ = 0.     (21)

The program solves this equation numerically by defining a variable

 E = r ργ−1     (22)

and finite differencing the equation according to

 Ei+1 − Ei + Ei−1 (Δ r)2
=
 4π G (γ − 1) γ K
r

 E r

 1/(γ−1)
.     (23)

### C.4  Linear wave

The subroutine `exact_wave` simply plots a sine function on a given graph. The function is of the form

 y = sin(k x − ω t)     (24)

where k is the wavenumber and ω is the angular frequency. These parameters are set via the input values of wavelength λ = 2π/k and wave period P = 2π/ω.

 λ wavelength P period
 Table 4: Input parameters for the linear wave exact solution

### C.5  Sedov blast wave

The subroutine `exact_sedov` computes the self-similar Sedov solution for a blast wave.

### C.6  Toy stars

The subroutine `exact_toystar1D` computes the exact solutions for the ‘Toy Stars’ described in ???. The system is one dimensional with velocity v, density ρ, and pressure P. The acceleration equation is

 dv dt
= −
 1 ρ

 ∂ P ∂ x
− Ω2 x,     (25)

We assume the equation of state is

 P = K ργ,     (26)

The exact solutions provided assume the equations are scaled such that Ω2 = 1.

#### C.6.1   Static structure

The static structure is given by

 ρ= 1− x2,     (27)

#### C.6.2   Linear solutions

The linear solution for the velocity is given by

 v = 0.05 Cs Gn(x) cosω t )     (28)

density is

 ρ = ρ + η     (29)

where

 η = 0.1 Cs ω Pn+1(x) sin(ω t))     (30)

#### C.6.3   Non-linear solution

In this case the velocity is given by

 v = A(t) x,     (31)

whilst the density solution is

 ργ −1 = H(t) − C(t) x2.     (32)

where the parameters A, H and C are determined by solving the ordinary differential equations

= −AH(γ −1),     (33)
Ȧ =

 2K γ γ −1
C − 1 − A2
(34)
Ċ = −AC(1+ γ),     (35)

The relation

A2 = −1 −
 2 σ C γ −1
+ kC
 2 γ +1

,     (36)

is used to check the quality of the solution of the differential equations by evaluating the constant k (which should remain close to its initial value).

### C.7  MHD shock tubes

These are some tabulated solutions for specific MHD shock tube problems at a given time taken from the tables given in ??? and ???.

### C.8  h vs ρ

The subroutine exact_hrho simply plots the relation between smoothing length and density, i.e.,

h = hfact

 m ρ

 1/ν
(37)

where ν is the number of spatial dimensions. The parameter hfact is output by the code into the header of each timestep. For particles of different masses, a different curve is plotted for each different mass value.

SPLASH: A visualisation tool for SPH data ©2004–2012 Daniel Price.
http://users.monash.edu.au/~dprice/splash/