Commit d200065c authored by James Willis's avatar James Willis
Browse files

Merge branch 'master' into dopair-vectorisation

Conflicts:
	src/Makefile.am
	src/cache.h
parents 06c0d64f 2781ef21
......@@ -22,9 +22,6 @@ doc/Doxyfile
examples/swift
examples/swift_mpi
examples/*.xmf
examples/used_parameters.yml
examples/energy.txt
examples/*/*.xmf
examples/*/*.hdf5
examples/*/*.png
......
......@@ -15,24 +15,26 @@ Usage: swift [OPTION]... PARAMFILE
swift_mpi [OPTION]... PARAMFILE
Valid options are:
-a Pin runners using processor affinity
-c Run with cosmological time integration
-C Run with cooling
-a Pin runners using processor affinity.
-c Run with cosmological time integration.
-C Run with cooling.
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
Allows user to check validy of parameter and IC files as well as memory limits.
-D Always drift all particles even the ones far from active particles.
-e Enable floating-point exceptions (debugging mode)
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements
-g Run with an external gravitational potential
-G Run with self-gravity
-D Always drift all particles even the ones far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-e Enable floating-point exceptions (debugging mode).
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential.
-G Run with self-gravity.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with SPH
-s Run with hydrodynamics.
-S Run with stars.
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-v [12] Increase the level of verbosity
1: MPI-rank 0 writes
2: All MPI-ranks write
-y {int} Time-step frequency at which task graphs are dumped
-h Print this help message and exit
-y {int} Time-step frequency at which task graphs are dumped.
-h Print this help message and exit.
See the file examples/parameter_example.yml for an example of parameter file.
......@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
AC_INIT([SWIFT],[0.4.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
AC_INIT([SWIFT],[0.5.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
......
......@@ -762,6 +762,7 @@ WARN_LOGFILE =
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
INPUT += @top_srcdir@/src/hydro/Minimal
INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/stars/Default
INPUT += @top_srcdir@/src/riemann
INPUT += @top_srcdir@/src/potential/point_mass
INPUT += @top_srcdir@/src/cooling/const_du
......
......@@ -91,6 +91,9 @@ for i in range(402):
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
print "Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0]
print "Ending energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1]
# Plot energy evolution
figure()
plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
......
......@@ -31,6 +31,7 @@ P0 = 0. # Constant additional pressure (should have no impact on the d
import matplotlib
matplotlib.use("Agg")
from pylab import *
from scipy import stats
import h5py
# Plot parameters
......@@ -104,6 +105,26 @@ u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
# Bin te data
r_bin_edge = np.arange(0., 1., 0.02)
r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
v_bin,_,_ = stats.binned_statistic(r, v_phi, statistic='mean', bins=r_bin_edge)
P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
v2_bin,_,_ = stats.binned_statistic(r, v_phi**2, statistic='mean', bins=r_bin_edge)
P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
# Plot the interesting quantities
figure()
......@@ -113,6 +134,7 @@ subplot(231)
plot(r, v_phi, '.', color='r', ms=0.5)
plot(solution_r, solution_v_phi, '--', color='k', alpha=0.8, lw=1.2)
errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
......@@ -125,6 +147,7 @@ subplot(232)
plot(r, rho, '.', color='r', ms=0.5)
plot(solution_r, solution_rho, '--', color='k', alpha=0.8, lw=1.2)
errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
......@@ -138,6 +161,7 @@ subplot(233)
plot(r, P, '.', color='r', ms=0.5)
plot(solution_r, solution_P, '--', color='k', alpha=0.8, lw=1.2)
errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
......@@ -150,6 +174,7 @@ subplot(234)
plot(r, u, '.', color='r', ms=0.5)
plot(solution_r, solution_u, '--', color='k', alpha=0.8, lw=1.2)
errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
......@@ -163,6 +188,7 @@ subplot(235)
plot(r, S, '.', color='r', ms=0.5)
plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
......
......@@ -4,7 +4,7 @@
if [ ! -e Isothermal.hdf5 ]
then
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 1000 1
python makeIC.py 1000 0
fi
rm -rf Isothermal_*.hdf5
......
......@@ -36,110 +36,154 @@ eta = 1.2349 # 48 ngbs with cubic spline kernel
rhoDM = 1.
Ldm = int(sys.argv[2]) # Number of particles along one axis
fileName = "multiTypes.hdf5"
massStars = 0.1
Lstars = int(sys.argv[3]) # Number of particles along one axis
fileBaseName = "multiTypes"
num_files = int(sys.argv[4])
#---------------------------------------------------
numGas = Lgas**3
massGas = boxSize**3 * rhoGas / numGas
numGas_tot = Lgas**3
massGas = boxSize**3 * rhoGas / numGas_tot
internalEnergy = P / ((gamma - 1.)*rhoGas)
numDM = Ldm**3
massDM = boxSize**3 * rhoDM / numDM
numDM_tot = Ldm**3
massDM = boxSize**3 * rhoDM / numDM_tot
numStars_tot = Lstars**3
massStars = massDM * massStars
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numGas, numDM, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
grp.attrs["Dimension"] = 3
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.
grp.attrs["Unit time in cgs (U_t)"] = 1.
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
# Gas Particle group
grp = file.create_group("/PartType0")
v = zeros((numGas, 3))
ds = grp.create_dataset('Velocities', (numGas, 3), 'f')
ds[()] = v
v = zeros(1)
m = full((numGas, 1), massGas)
ds = grp.create_dataset('Masses', (numGas,1), 'f')
ds[()] = m
m = zeros(1)
h = full((numGas, 1), eta * boxSize / Lgas)
ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
ds[()] = h
h = zeros(1)
u = full((numGas, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numGas,1), 'f')
ds[()] = u
u = zeros(1)
ids = linspace(0, numGas, numGas, endpoint=False).reshape((numGas,1))
ds = grp.create_dataset('ParticleIDs', (numGas, 1), 'L')
ds[()] = ids + 1
x = ids % Lgas;
y = ((ids - x) / Lgas) % Lgas;
z = (ids - x - Lgas * y) / Lgas**2;
coords = zeros((numGas, 3))
coords[:,0] = z[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
coords[:,1] = y[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
coords[:,2] = x[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
ds = grp.create_dataset('Coordinates', (numGas, 3), 'd')
ds[()] = coords
# DM Particle group
grp = file.create_group("/PartType1")
v = zeros((numDM, 3))
ds = grp.create_dataset('Velocities', (numDM, 3), 'f')
ds[()] = v
v = zeros(1)
m = full((numDM, 1), massDM)
ds = grp.create_dataset('Masses', (numDM,1), 'f')
ds[()] = m
m = zeros(1)
ids = linspace(0, numDM, numDM, endpoint=False).reshape((numDM,1))
ds = grp.create_dataset('ParticleIDs', (numDM, 1), 'L')
ds[()] = ids + Lgas**3 + 1
x = ids % Ldm;
y = ((ids - x) / Ldm) % Ldm;
z = (ids - x - Ldm * y) / Ldm**2;
coords = zeros((numDM, 3))
coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
ds = grp.create_dataset('Coordinates', (numDM, 3), 'd')
ds[()] = coords
file.close()
offsetGas = 0
offsetDM = 0
offsetStars = 0
for n in range(num_files):
# File name
if num_files == 1:
fileName = fileBaseName + ".hdf5"
else:
fileName = fileBaseName + ".%d.hdf5"%n
# File
file = h5py.File(fileName, 'w')
# Number of particles
numGas = numGas_tot / num_files
numDM = numDM_tot / num_files
numStars = numStars_tot / num_files
if n == num_files - 1:
numGas += numGas_tot % num_files
numDM += numDM_tot % num_files
numStars += numStars_tot % num_files
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numGas_tot, numDM_tot, 0, 0, numStars_tot, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, numStars, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = num_files
grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
grp.attrs["Dimension"] = 3
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.
grp.attrs["Unit time in cgs (U_t)"] = 1.
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
# Gas Particle group
grp = file.create_group("/PartType0")
v = zeros((numGas, 3))
ds = grp.create_dataset('Velocities', (numGas, 3), 'f', data=v)
m = full((numGas, 1), massGas)
ds = grp.create_dataset('Masses', (numGas,1), 'f', data=m)
h = full((numGas, 1), eta * boxSize / Lgas)
ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f', data=h)
u = full((numGas, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numGas,1), 'f', data=u)
ids = linspace(offsetGas, offsetGas+numGas, numGas, endpoint=False).reshape((numGas,1))
ds = grp.create_dataset('ParticleIDs', (numGas, 1), 'L', data=ids+1)
x = ids % Lgas;
y = ((ids - x) / Lgas) % Lgas;
z = (ids - x - Lgas * y) / Lgas**2;
coords = zeros((numGas, 3))
coords[:,0] = z[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
coords[:,1] = y[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
coords[:,2] = x[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
ds = grp.create_dataset('Coordinates', (numGas, 3), 'd', data=coords)
# DM Particle group
grp = file.create_group("/PartType1")
v = zeros((numDM, 3))
ds = grp.create_dataset('Velocities', (numDM, 3), 'f', data=v)
m = full((numDM, 1), massDM)
ds = grp.create_dataset('Masses', (numDM,1), 'f', data=m)
ids = linspace(offsetDM, offsetDM+numDM, numDM, endpoint=False).reshape((numDM,1))
ds = grp.create_dataset('ParticleIDs', (numDM, 1), 'L', data=ids + numGas_tot + 1)
ds[()] = ids + Lgas**3 + 1
x = ids % Ldm;
y = ((ids - x) / Ldm) % Ldm;
z = (ids - x - Ldm * y) / Ldm**2;
coords = zeros((numDM, 3))
coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
ds = grp.create_dataset('Coordinates', (numDM, 3), 'd', data=coords)
# Star Particle group
grp = file.create_group("/PartType4")
v = zeros((numStars, 3))
ds = grp.create_dataset('Velocities', (numStars, 3), 'f', data=v)
m = full((numStars, 1), massStars)
ds = grp.create_dataset('Masses', (numStars,1), 'f', data=m)
ids = linspace(0, numStars, numStars, endpoint=False).reshape((numStars,1))
ds = grp.create_dataset('ParticleIDs', (numStars, 1), 'L', data=ids + numGas_tot + numDM_tot + 1)
x = ids % Ldm;
y = ((ids - x) / Ldm) % Ldm;
z = (ids - x - Ldm * y) / Ldm**2;
coords = zeros((numStars, 3))
coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
ds = grp.create_dataset('Coordinates', (numStars, 3), 'd', data=coords)
# Shift stuff
offsetGas += numGas
offsetDM += numDM
offsetStars += numStars
file.close()
......@@ -32,6 +32,7 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./multiTypes.hdf5 # The file to read
replicate: 2 # Replicate all particles twice along each axis
# External potential parameters
PointMassPotential:
......
......@@ -4,7 +4,7 @@
if [ ! -e multiTypes.hdf5 ]
then
echo "Generating initial conditions for the multitype box example..."
python makeIC.py 50 60
python makeIC.py 9 13 7 1
fi
../swift -s -g -t 16 multiTypes.yml 2>&1 | tee output.log
../swift -s -g -S -t 1 multiTypes.yml 2>&1 | tee output.log
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import h5py
from numpy import *
# Generates a swift IC file for the 1D Noh problem test in a periodic box
# Parameters
numPart = 1001
gamma = 5./3. # Gas adiabatic index
rho0 = 1. # Background density
P0 = 1.e-6 # Background pressure
fileName = "noh.hdf5"
#---------------------------------------------------
coords = zeros((numPart, 3))
h = zeros(numPart)
vol = 2.
for i in range(numPart):
coords[i,0] = i * vol/numPart + vol/(2.*numPart)
h[i] = 1.2348 * vol / numPart
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m = zeros(numPart)
u = zeros(numPart)
m[:] = rho0 * vol / numPart
u[:] = P0 / (rho0 * (gamma - 1))
v[coords[:,0]<vol/2. ,0] = 1
v[coords[:,0]>vol/2. ,0] = -1
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = [vol, vol, vol]
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
grp.attrs["Dimension"] = 1
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 1
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.
grp.attrs["Unit time in cgs (U_t)"] = 1.
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
UnitLength_in_cgs: 1 # Centimeters
UnitVelocity_in_cgs: 1 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.6 # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: noh # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 5e-2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-5 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions
InitialConditions:
file_name: ./noh.hdf5 # The file to read
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
# Computes the analytical solution of the Noh problem and plots the SPH answer
# Parameters
gas_gamma = 5./3. # Polytropic index
rho0 = 1. # Background density
P0 = 1.e-6 # Background pressure